tensor() function and sets
Jonathan Rougier
Jonathan Rougier <J.C.Rougier@durham.ac.uk>
Tue, 20 Jul 1999 15:26:35 +0100 (BST)
Hi Everyone,
To complete the outer() and kronecker() functions in the base, may I
suggest the following tensor() function, which allows the multiplication
of arrays through sets of conformable dimensions. I am happy to write a
help page if required.
The code also needs a setdiff() function which prompts me to ask: what
about simple set functions? I expect many of us have written our own
(Brian has a setdiff() in drop1.lm(), for example), which seems like a
good reason for putting versions in the base. I would be happy to provide
mine for general scrutiny.
Cheers, Jonathan.
Jonathan Rougier Science Laboratories
Department of Mathematical Sciences South Road
University of Durham Durham DH1 3LE
"setdiff" <-
function (x, y) x[!(x %in% y)]
"tensor" <-
function (A, B, da, db)
{
# tensor product of A and B through da and db
no.na <- is.null(na <- dimnames(A <- as.array(A)))
dima <- dim(A)
no.nb <- is.null(nb <- dimnames(B <- as.array(B)))
dimb <- dim(B)
if (any(dima[da] != dimb[db]))
stop("Mismatched dimensions")
kpa <- setdiff(seq(along = dima), da)
kpb <- setdiff(seq(along = dimb), db)
# fix up the dimnames (see outer)
if (no.na && no.nb)
nms <- NULL
else {
if (no.na)
na <- vector("list", length(dima))
else if (no.nb)
nb <- vector("list", length(dimb))
nms <- c(na[kpa], nb[kpb])
}
A <- matrix(aperm(A, c(kpa, da)), ncol = prod(dima[da]))
B <- matrix(aperm(B, c(db, kpb)), nrow = prod(dimb[db]))
array(A %*% B, c(dima[kpa], dimb[kpb]), dimnames = nms)
}
# examples
A <- array(1:6, 2:3, list(LETTERS[1:2], LETTERS[3:5]))
B <- 1:3
tensor(A, B, 2, 1) # same as drop(A %*% B)
A <- outer(A, B) # A now 2 by 3 by 3
tensor(A, B, 2, 1); tensor(A, B, 3, 1) # both 2 by 3 (nb dimnames)
B <- outer(1:2, B) # B now 2 by 3
tensor(A, B, c(1, 3), 1:2) # must be 3 vector
B <- outer(1:3, B) # B now 3 by 2 by 3
tensor(A, B, 1:3, c(2, 1, 3)) # must be a scalar
sum(A * aperm(B, c(2, 1, 3))) # same scalar value
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._