# tensor() function and sets

Martin Maechler Martin Maechler <maechler@stat.math.ethz.ch>
Wed, 21 Jul 1999 18:44:18 +0200

```>>>>> On Tue, 20 Jul 1999 15:26:35 +0100 (BST), Jonathan Rougier <J.C.Rougier@durham.ac.uk> said:

JonR> Hi Everyone,

JonR> To complete the outer() and kronecker() functions in the base, may I
JonR> suggest the following tensor() function, which allows the multiplication
JonR> of arrays through sets of conformable dimensions.  I am happy to write a
JonR> help page if required.

JonR> The code also needs a setdiff() function which prompts me to ask:
JonR> what about simple set functions?  I expect many of us have written our own
JonR> (Brian has a setdiff() in drop1.lm(), for example), which seems like a
JonR> good reason for putting versions in the base.  I would be happy to provide
JonR> mine for general scrutiny.

Actually your version of "setdiff" is in
example("%in%")

"%w/o%" <- function(x,y) x[!x %in% y] #--  x without y

Maybe we should make this part of core R, since it allows the nicer syntax
(1:10) %w/o% c(3,7,12).

One problem with all these might be that you forget to think about efficiency
{ %w/o% -> %in% -> match(.) }
but that might be ok.

union
and
intersect
are two other set functions, the first of which is also contained in
example(match) { == example("%in%") }:
intersect <- function(x, y) y[match(x, y, nomatch = 0)]

What do others think?  Should (some of) these be put into "core R" ?

JonR> Jonathan Rougier                       Science Laboratories
JonR> Department of Mathematical Sciences    South Road
JonR> University of Durham                   Durham DH1 3LE

JonR> "setdiff" <-
JonR> function (x, y) x[!(x %in% y)]

JonR> "tensor" <-
JonR> function (A, B, da, db)
JonR> {
JonR> # tensor product of A and B through da and db
JonR> no.na <- is.null(na <- dimnames(A <- as.array(A)))
JonR> dima <- dim(A)
JonR> no.nb <- is.null(nb <- dimnames(B <- as.array(B)))
JonR> dimb <- dim(B)
JonR> if (any(dima[da] != dimb[db]))
JonR> stop("Mismatched dimensions")
JonR> kpa <- setdiff(seq(along = dima), da)
JonR> kpb <- setdiff(seq(along = dimb), db)
JonR> # fix up the dimnames (see outer)
JonR> if (no.na && no.nb)
JonR> nms <- NULL
JonR> else {
JonR> if (no.na)
JonR> na <- vector("list", length(dima))
JonR> else if (no.nb)
JonR> nb <- vector("list", length(dimb))
JonR> nms <- c(na[kpa], nb[kpb])
JonR> }
JonR> A <- matrix(aperm(A, c(kpa, da)), ncol = prod(dima[da]))
JonR> B <- matrix(aperm(B, c(db, kpb)), nrow = prod(dimb[db]))
JonR> array(A %*% B, c(dima[kpa], dimb[kpb]), dimnames = nms)
JonR> }

JonR> # examples

JonR> A <- array(1:6, 2:3, list(LETTERS[1:2], LETTERS[3:5]))
JonR> B <- 1:3
JonR> tensor(A, B, 2, 1)			# same as drop(A %*% B)
JonR> A <- outer(A, B)			# A now 2 by 3 by 3
JonR> tensor(A, B, 2, 1); tensor(A, B, 3, 1)	# both 2 by 3 (nb dimnames)
JonR> B <- outer(1:2, B)			# B now 2 by 3
JonR> tensor(A, B, c(1, 3), 1:2)		# must be 3 vector
JonR> B <- outer(1:3, B)			# B now 3 by 2 by 3
JonR> tensor(A, B, 1:3, c(2, 1, 3))		# must be a scalar
JonR> sum(A * aperm(B, c(2, 1, 3)))		# same scalar value

Looks interesting, however is it really useful
{for doing statistics, graphics, numerical computing, I mean} ?

Yes, if you'd provide a *.Rd help file  (use  ``prompt(tensor)''),
we'd be even more inclined to include the tensor function....
Also, I'd like to see a "real" use for it.

Is
tensor(A, B, 1:3, c(2, 1, 3))
really easier than
sum(A * aperm(B, c(2, 1, 3)))

[probably yes, you have to think less, but still:
couldn't  "da" and "db" (3rd and 4th arg) get reasonable default values ?]

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

```