tensor() function and sets

[Jonathan Rougier]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._