tensor() function and sets

About this list Date view Thread view Subject view Author view Other groups

Subject: tensor() function and sets
From: Jonathan Rougier (J.C.Rougier@durham.ac.uk)
Date: Wed 21 Jul 1999 - 00:26:35 EST


Message-ID: <Pine.GSO.3.96.990720143125.9359C-100000@laplace>

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


About this list Date view Thread view Subject view Author view Other groups

This archive was generated by hypermail 2b25 : Tue 04 Jan 2000 - 14:16:06 EST