Extension to outer()

Jonathan Rougier J.C.Rougier@durham.ac.uk
Mon, 22 Mar 1999 12:40:03 +0000 (GMT)


  This message is in MIME format.  The first part should be readable text,
  while the remaining parts are likely unreadable without MIME-aware tools.
  Send mail to mime@docserver.cac.washington.edu for more info.

---559023410-1804928587-922106403=:691
Content-Type: TEXT/PLAIN; charset=US-ASCII

A couple of weeks ago I asked if anyone had generalised outer to include
bound extents -- effectively taking the diagonal of an outer product.
Attached is my take on the solution.  There is an extra argument to outer,
BIND, which if null results in a simple outer product.  I've used a
slightly different algorithm which I think is more transparent than the
original code, and also generalises better.

If BIND is not null, then outer(A, B, BIND = c(i, j)) computes an array
with extent j of B bound along extent i of A.  For example, if A and B are
both matrices and V is the result with BIND=c(1,1), then V[i, j, k] =
FUN(A[i, j], B[i, k], ...).  One use for the BIND argument is in computing
Euclidean distances:

"gdist" <- function(X, Y)
{
	X <- as.matrix(X)
	Y <- as.matrix(Y)
	sqrt(apply(outer(X, Y, "-", BIND = c(2,2))^2, c(1,3), sum))
}

which, unlike dist() from `mva' can provide submatrices of the full
distance matrix -- very useful for covariances.

This BIND argument to outer turns out to be very useful in lots of
different areas of our work.  If anyone has any feedback on the code I
would be very grateful for it.  Many thanks, Jonathan. 

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

"[B]egin upon the precept ... that the things we see are to be 
 weighed in the scale with what we know"  (Meredith, 1879, The Egoist)

---559023410-1804928587-922106403=:691
Content-Type: TEXT/PLAIN; charset=US-ASCII; name="outer.R"
Content-Transfer-Encoding: BASE64
Content-ID: <Pine.GSO.3.96.990322124003.691B@laplace>
Content-Description: 

Im91dGVyIiA8LQ0KZnVuY3Rpb24gKEEsIEIsIEZVTiA9ICIqIiwgQklORCA9
IE5VTEwsIC4uLikgDQp7DQogICAgRlVOIDwtIG1hdGNoLmZ1bihGVU4pDQog
ICAgQSA8LSBhcy5hcnJheShBKQ0KICAgIEIgPC0gYXMuYXJyYXkoQikNCiAg
ICBkaW1zIDwtIGxpc3QoZGltKEEpLCBkaW0oQikpDQogICAgbmEgPC0gbGVu
Z3RoKGRpbXNbWzFdXSkNCiAgICBuYiA8LSBsZW5ndGgoZGltc1tbMl1dKQ0K
ICAgIGlmIChpcy5udWxsKGRuYSA8LSBkaW1uYW1lcyhBKSkpIA0KICAgICAg
ICBkbmEgPC0gdmVjdG9yKCJsaXN0IiwgbmEpDQogICAgaWYgKGlzLm51bGwo
ZG5iIDwtIGRpbW5hbWVzKEIpKSkgDQogICAgICAgIGRuYiA8LSB2ZWN0b3Io
Imxpc3QiLCBuYikNCiAgICBpZiAoaXMubnVsbChCSU5EKSkgew0KICAgICAg
ICAjIHNpbXBsZSBvdXRlcg0KICAgICAgICBBIDwtIGFycmF5KEEsIGMoZGlt
c1tbMV1dLCBkaW1zW1syXV0pKQ0KICAgICAgICBCIDwtIGFycmF5KEIsIGMo
ZGltc1tbMl1dLCBkaW1zW1sxXV0pKQ0KICAgICAgICBCIDwtIGFwZXJtKEIs
IGMobmIgKyAxOm5hLCAxOm5iKSkNCiAgICAgICAgcm9iaiA8LSBhcnJheShG
VU4oQSwgQiwgLi4uKSwgZGltKEEpKQ0KICAgICAgICBkaW1uYW1lcyhyb2Jq
KSA8LSBjKGRuYSwgZG5iKQ0KICAgICAgICByb2JqDQogICAgfQ0KICAgIGVs
c2Ugew0KICAgICAgICAjIG1vcmUgY29tcGxpY2F0ZWQgYm91bmQgb3V0ZXIN
CiAgICAgICAgaWYgKCEoQklORFsxXSAlaW4lIDE6bmEgJiYgQklORFsyXSAl
aW4lIDE6bmIpKSANCiAgICAgICAgICAgIHN0b3AoImludmFsaWQgZXh0ZW50
IGZvciBcIkJJTkRcIiIpDQogICAgICAgIGVsc2UgaWYgKGRpbXNbWzFdXVtC
SU5EWzFdXSAhPSBkaW1zW1syXV1bQklORFsyXV0pIA0KICAgICAgICAgICAg
c3RvcCgiRXh0ZW50cyB0byBiZSBib3VuZCBkbyBub3QgbWF0Y2giKQ0KICAg
ICAgICBpZiAobmIgPT0gMSkgDQogICAgICAgICAgICByb2JqIDwtIHN3ZWVw
KEEsIEJJTkRbMV0sIEIsIEZVTiA9IEZVTiwgLi4uKQ0KICAgICAgICBlbHNl
IHsNCiAgICAgICAgICAgICMgc3dpdGNoIGJvdW5kIGV4dGVudCB0byB0aGUg
ZnJvbnQNCiAgICAgICAgICAgIG9kYSA8LSBjKEJJTkRbMV0sICgxOm5hKVst
QklORFsxXV0pDQogICAgICAgICAgICBkYSA8LSBkaW0oQSA8LSBhcGVybShB
LCBvZGEpKQ0KICAgICAgICAgICAgZGIgPC0gZGltKEIgPC0gYXBlcm0oQiwg
YyhCSU5EWzJdLCAoMTpuYilbLUJJTkRbMl1dKSkpDQogICAgICAgICAgICAj
IHJlcGxpY2F0ZSBBIGFsb25nIEIgKGV4Y2x1ZGluZyBib3VuZCBleHRleHQp
IGFuZCB2aWNlIHZlcnNhDQogICAgICAgICAgICBBIDwtIGFycmF5KEEsIGMo
ZGEsIGRiWy0xXSkpDQogICAgICAgICAgICBCIDwtIGFycmF5KEIsIGMoZGIs
IGRhWy0xXSkpDQogICAgICAgICAgICAjIGB0cmFuc3Bvc2UnIEIsIGNvbXB1
dGUgdGhlbiB1bi1wZXJtDQogICAgICAgICAgICBpZiAobmEgPiAxKSANCiAg
ICAgICAgICAgICAgICBCIDwtIGFwZXJtKEIsIGMoMSwgbmIgKyAxOihuYSAt
IDEpLCAyOm5iKSkNCiAgICAgICAgICAgIHJvYmogPC0gRlVOKEEsIEIsIC4u
LikNCiAgICAgICAgICAgIGRpbShyb2JqKSA8LSBjKGRpbXNbWzFdXVtvZGFd
LCBkaW1zW1syXV1bLUJJTkRbMl1dKQ0KICAgICAgICAgICAgcm9iaiA8LSBh
cGVybShyb2JqLCBjKG9yZGVyKG9kYSksIG5hICsgMToobmIgLSAxKSkpDQog
ICAgICAgIH0NCiAgICAgICAgZGltbmFtZXMocm9iaikgPC0gYyhkbmEsIGRu
YlstQklORFsyXV0pDQogICAgICAgIHJvYmoNCiAgICB9DQp9DQo=
---559023410-1804928587-922106403=:691--
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._