[R] A slightly unorthodox matrix product.

[Rolf Turner]

On 05/08/18 16:41, Thomas Jagger wrote:
>  > Date: Sat, 4 Aug 2018 17:52:37 +1200
>  >
>  > From: Rolf Turner <r.turner using auckland.ac.nz 
> <mailto:r.turner using auckland.ac.nz>>
>  > To: r-help <r-help using r-project.org <mailto:r-help using r-project.org>>
>  > Subject: [R] A slightly unorthodox matrix product.
>  > Message-ID: <e0073e07-31b3-20df-c20d-6c565c857554 using auckland.ac.nz 
> <mailto:e0073e07-31b3-20df-c20d-6c565c857554 using auckland.ac.nz>>
>  > Content-Type: text/plain; charset="utf-8"; Format="flowed"
>  >
>  >
>  > Can anyone think of a sexy way of forming following "product"?
>  >
>  > Given matrices A and B, both with m rows, form a 3 dimensional array C
>  > such that:
>  >
>  >      C[i,j,k] = A[i,j]*B[i,k]
>  >
>  > I *think* that the following does what I want.  (I keep confusing
>  > myself, so I'm not sure!)
>  >
>  > library(abind)
>  > xxx <- lapply(1:nrow(a),function(i,a,b){a[i,]%o%b[i,]},a=A,b=B)
>  > do.call(abind,c(xxx,list(along=3)))
>  >
>  > Is there a cleverer way?
>  >
>  > cheers,
>  >
>  > Rolf Turner
>  >
>  > --
>  > Technical Editor ANZJS
>  > Department of Statistics
>  > University of Auckland
>  > Phone: +64-9-373-7599 ext. 88276
> 
> Dear Rolf,
> Try the following:
> 
> B<-matrix(1:12,3,4)
> 
> C<-as.vector(A[,rep(seq(ncol(A)),ncol(B))])*as.vector(B[,rep(seq(ncol(B)),each=ncol(A))])
> dim(C) <- c(nrow(A),ncol(A),ncol(B))
> 
> #test it on column 2 should return true
> all(C[,,2]==A*B[,rep(2,ncol(A))])
> #on all columns (sapply returns 9 rows with 3 columns all values are TRUE)
> 
> all( sapply(seq(ncol(C)),function(i) (C[,,i]==A*B[,rep(i,ncol(A))]) ) )
> 
> Note that it creates the final array by taking advantage of the 
> column-major ordering in R.
> Initially, we create a vector by multiplying elementwise the 2 vectors 
> internally associated with each matrix,
> finally,  we generate our  3D array by adding the dimensions attribute, 
> a vector of 3  elements.
> 
> This method should be fairly fast since we are using internal R matrix 
> addressing, and not multiple function calls required by lapply ().
> I hope this helps

Neat, and much better than anything I could have thought of.  However
microbenchmark() indicates that the Chuck Berry/Jeff Newmiller solution 
is about twice as fast.  Not that this speed difference is Any Big Deal, 
but.

Thanks for taking an interest in this obscure query of mine.

cheers,

Rolf
-- 
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276