R help - Aug 2018 - A slightly unorthodox matrix product.

> Date: Sat, 4 Aug 2018 17:52:37 +1200
>
> From: Rolf Turner <r.turner at auckland.ac.nz>
> To: r-help <r-help at r-project.org>
> Subject: [R] A slightly unorthodox matrix product.
> Message-ID: <e0073e07-31b3-20df-c20d-6c565c857554 at 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
Thomas Jagger

	[[alternative HTML version deleted]]

On 05/08/18 16:41, Thomas Jagger wrote:
>  > Date: Sat, 4 Aug 2018 17:52:37 +1200
>  >
>  > From: Rolf Turner <r.turner at auckland.ac.nz 
> <mailto:r.turner at auckland.ac.nz>>
>  > To: r-help <r-help at r-project.org <mailto:r-help at
r-project.org>>
>  > Subject: [R] A slightly unorthodox matrix product.
>  > Message-ID: <e0073e07-31b3-20df-c20d-6c565c857554 at
auckland.ac.nz
> <mailto:e0073e07-31b3-20df-c20d-6c565c857554 at 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