A slightly unorthodox matrix product.

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

A slightly unorthodox matrix product.

Thomas Jagger
> Date: Sat, 4 Aug 2018 17:52:37 +1200
>
> From: Rolf Turner <[hidden email]>
> To: r-help <[hidden email]>
> Subject: [R] A slightly unorthodox matrix product.
> Message-ID: <[hidden email]>
> 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]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[hidden email]
> <mailto:[hidden email]>>
>  > To: r-help <[hidden email] <mailto:[hidden email]>>
>  > Subject: [R] A slightly unorthodox matrix product.
>  > Message-ID: <[hidden email]
> <mailto:[hidden email]>>
>  > 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

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.