# A slightly unorthodox matrix product.

8 messages
Open this post in threaded view
|

## A slightly unorthodox matrix product.

 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 ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: A slightly unorthodox matrix product.

 Hi Rolf, A few edits because (i) nrow(a) should be nrow(A) and (ii) you have calculated C[j,k,i] = A[i,j]*B[i,k], (iii) minor style change on lapply. library(abind) xxx <- lapply(1:nrow(A),function(i){A[i,]%o%B[i,]}) yyy <- do.call(abind,c(xxx,list(along=3))) zzz <- aperm(yyy,c(3,1,2)) HTH, Eric On Sat, Aug 4, 2018 at 8:52 AM, Rolf Turner <[hidden email]> wrote: > > 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 > > ______________________________________________ > [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/posti> ng-guide.html > and provide commented, minimal, self-contained, reproducible code. >         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: A slightly unorthodox matrix product.

 > On Aug 4, 2018, at 10:01 AM, Eric Berger <[hidden email]> wrote: > > Hi Rolf, > A few edits because (i) nrow(a) should be nrow(A) and (ii) you have > calculated C[j,k,i] = A[i,j]*B[i,k], (iii) minor style change on lapply. > > library(abind) > xxx <- lapply(1:nrow(A),function(i){A[i,]%o%B[i,]}) > yyy <- do.call(abind,c(xxx,list(along=3))) Or use the simplify="array" gambit in sapply: yyy <- sapply(1:nrow(A), function(i) A[i,] %o% B[i,], simplify="array") > zzz <- aperm(yyy,c(3,1,2)) > HTH, Chuck ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: A slightly unorthodox matrix product.

 Sometimes a good old for loop performs best, even if it doesn't look sexy: ######## A <- matrix( 1:12, nrow = 3 ) B <- matrix( 1:15, nrow = 3 ) library(abind) # Eric ans1 <- function( a, b ) {    xxx <- lapply( seq.int( nrow( A ) )                 , function( i ) {                      A[ i, ] %o% B[ i, ]                   }                 )    yyy <- do.call( abind, c( xxx, list( along = 3 ) ) )    zzz <- aperm( yyy, c( 3, 1, 2 ) )    zzz } # Charles ans1b <- function( a, b ) {    xxx <- lapply( seq.int( nrow( A ) )                 , function( i ) {                      A[ i, ] %o% B[ i, ]                   }                 )    yyy <- sapply( seq.int( nrow( a ) )                 , function( i ) a[ i, ] %o% b[ i, ]                 , simplify = "array"                 )    zzz <- aperm( yyy, c( 3, 1, 2 ) )    zzz } # Jeff #1 ans2 <- function( a, b ) {    zzz <- array( rep( NA, nrow( a ) * ncol( a ) * ncol( b ) )                , dim = c( nrow( a ), ncol( a ), ncol( b ) )                )    jseq <- seq.int( ncol( a ) )    kseq <- seq.int( ncol( b ) )    for ( i in seq.int( nrow( a ) ) ) {      zzz[ i, jseq, kseq ] <- outer( a[ i, ], b[ i, ] )    }    zzz } # Jeff #2 ans3 <- function( a, b ) {    idxs <- expand.grid( i = seq.int( nrow( a ) )                       , j = seq.int( ncol( a ) )                       , k = seq.int( ncol( b ) )                       )    ij <- as.matrix( idxs[ , c( "i", "j" ) ] )    ik <- as.matrix( idxs[ , c( "i", "k" ) ] )    array( a[ ij ] * b[ ik ]         , dim = c( nrow( a ), ncol( a ), ncol( b ) )         ) } library(microbenchmark) microbenchmark( res1 <- ans1( A, B )                , res1b <- ans1b( A, B )                , res2 <- ans2( A, B )                , res3 <- ans3( A, B )                ) #> Unit: microseconds #>                  expr     min       lq      mean   median       uq #>    res1 <- ans1(A, B) 660.489 688.3460 4199.5385 742.5505 805.1860 #>  res1b <- ans1b(A, B) 224.436 236.2250  427.4806 246.3240 269.6425 #>    res2 <- ans2(A, B)  91.538  96.9075  287.9596 102.0335 110.8825 #>    res3 <- ans3(A, B) 508.642 528.9700  860.6295 563.5470 619.5285 #>        max neval #>  344769.27   100 #>   17062.63   100 #>   18212.11   100 #>   23041.89   100 all( res1 == res2 ) #> [1] TRUE all( res1 == res1b ) #> [1] TRUE all( res1 == res3 ) #> [1] TRUE res3 #> , , 1 #> #>      [,1] [,2] [,3] [,4] #> [1,]    1    4    7   10 #> [2,]    4   10   16   22 #> [3,]    9   18   27   36 #> #> , , 2 #> #>      [,1] [,2] [,3] [,4] #> [1,]    4   16   28   40 #> [2,]   10   25   40   55 #> [3,]   18   36   54   72 #> #> , , 3 #> #>      [,1] [,2] [,3] [,4] #> [1,]    7   28   49   70 #> [2,]   16   40   64   88 #> [3,]   27   54   81  108 #> #> , , 4 #> #>      [,1] [,2] [,3] [,4] #> [1,]   10   40   70  100 #> [2,]   22   55   88  121 #> [3,]   36   72  108  144 #> #> , , 5 #> #>      [,1] [,2] [,3] [,4] #> [1,]   13   52   91  130 #> [2,]   28   70  112  154 #> [3,]   45   90  135  180 #' Created on 2018-08-04 by the [reprex package](http://reprex.tidyverse.org) (v0.2.0). ######## On Sat, 4 Aug 2018, Berry, Charles wrote: > > >> On Aug 4, 2018, at 10:01 AM, Eric Berger <[hidden email]> wrote: >> >> Hi Rolf, >> A few edits because (i) nrow(a) should be nrow(A) and (ii) you have >> calculated C[j,k,i] = A[i,j]*B[i,k], (iii) minor style change on lapply. >> >> library(abind) >> xxx <- lapply(1:nrow(A),function(i){A[i,]%o%B[i,]}) >> yyy <- do.call(abind,c(xxx,list(along=3))) > > Or use the simplify="array" gambit in sapply: > > yyy <- sapply(1:nrow(A), function(i) A[i,] %o% B[i,], simplify="array") > >> zzz <- aperm(yyy,c(3,1,2)) >> > > HTH, > > Chuck > > ______________________________________________ > [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. > --------------------------------------------------------------------------- Jeff Newmiller                        The     .....       .....  Go Live... DCN:<[hidden email]>        Basics: ##.#.       ##.#.  Live Go...                                        Live:   OO#.. Dead: OO#..  Playing Research Engineer (Solar/Batteries            O.O#.       #.O#.  with /Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: A slightly unorthodox matrix product.

 > On Aug 4, 2018, at 11:43 AM, Jeff Newmiller <[hidden email]> wrote: > > Sometimes a good old for loop performs best, even if it doesn't look sexy: > > Fair enough, but a vectorized solution beats them all (see below). Also, [SNIP] > # Charles > ans1b <- function( a, b ) > { The lapply you put here was from Eric's solution: >  xxx <- lapply( seq.int( nrow( A ) ) >               , function( i ) { >                    A[ i, ] %o% B[ i, ] >                 } This is what I had in mind: ans1b.corrected <- function( a, b ) {   yyy <- sapply( seq.int( nrow( a ) )                  , function( i ) a[ i, ] %o% b[ i, ]                  , simplify = "array"   )   zzz <- aperm( yyy, c( 3, 1, 2 ) )   zzz } On my system it is slower than a for loop but a lot faster than your benchmark showed with the superfluous code from Eric's solution. For speed, a vectorized solution is faster than a for loop by a factor of 3 on my laptop: ans0 <- function(A,B){   nca <- ncol(A)   ncb <- ncol(B)   j.index <- rep(1:nca, times = ncb)   k.index <- rep(1:nca, each = ncb)   res <- array(A[, j.index] * B[, k.index], c(nrow(A), nca, ncb))   res   } > microbenchmark( +   res0 <- ans0(A, B), +   res1 <- ans1(A, B), +   res1b <- ans1b.corrected(A, B), +   res2 <- ans2(A, B), +   res3 <- ans3(A, B) + ) Unit: microseconds                            expr     min       lq      mean   median       uq     max neval   cld              res0 <- ans0(A, B)  13.281  18.4960  21.52723  19.9905  23.4750  61.556   100 a                  res1 <- ans1(A, B) 353.121 369.8635 409.77788 381.5840 444.3290 701.256   100     e  res1b <- ans1b.corrected(A, B)  82.816  89.4185 101.52321  95.4275 107.1700 217.357   100   c                res2 <- ans2(A, B)  49.674  54.4825  61.78278  58.7540  65.5265 172.625   100  b                res3 <- ans3(A, B) 317.772 342.4220 392.25065 360.4675 436.2125 602.346   100    d > FWIW, if there are dimnames on A and B, sapply( row names(A), ..., simplify="array") preserves them without further ado. Chuck ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: A slightly unorthodox matrix product.

 Sorry I missed your intent on the sapply. Slick work on the vectorizing, but for the future reference it was slightly buggy: ####### A <- matrix( 1:12, nrow = 3 ) B <- matrix( 1:15, nrow = 3 ) # for loop ans2 <- function( a, b ) {    zzz <- array( rep( NA, nrow( a ) * ncol( a ) * ncol( b ) )                , dim = c( nrow( a ), ncol( a ), ncol( b ) )                )    jseq <- seq.int( ncol( a ) )    kseq <- seq.int( ncol( b ) )    for ( i in seq.int( nrow( a ) ) ) {      zzz[ i, jseq, kseq ] <- outer( a[ i, ], b[ i, ] )    }    zzz } # fast but buggy ans0 <- function( A, B ) {    nca <- ncol( A )    ncb <- ncol( B )    j.index <- rep( seq.int( nca ), times = ncb)    k.index <- rep( seq.int( nca ), each = ncb)    res <- array( A[ , j.index ] * B[ , k.index ]                , c( nrow( A ), nca, ncb )                )    res    } # bugfixed ans0b <- function( A, B ) {    nca <- ncol( A )    ncb <- ncol( B )    j.index <- rep( seq.int( nca ), times = ncb )    k.index <- rep( seq.int( ncb ), each = nca )    res <- array( A[ , j.index ] * B[ , k.index ]                , c( nrow( A ), nca, ncb )                )    res    } library(microbenchmark) microbenchmark( res2 <- ans2( A, B )                , res0b <- ans0b( A, B )                , res0 <- ans0( A, B )                ) #> Unit: microseconds #>                  expr    min      lq     mean  median      uq      max #>    res2 <- ans2(A, B) 84.987 87.8185 270.2153 96.4315 99.4175 17531.77 #>  res0b <- ans0b(A, B) 17.940 19.2055 126.8974 20.8800 22.2865 10616.36 #>    res0 <- ans0(A, B) 18.041 19.1670 126.1183 20.5530 21.9545 10532.44 #>  neval #>    100 #>    100 #>    100 all( res2 == res0 ) #> [1] FALSE all( res2 == res0b ) #> [1] TRUE #' Created on 2018-08-04 by the [reprex package](http://reprex.tidyverse.org) (v0.2.0). ####### On Sat, 4 Aug 2018, Berry, Charles wrote: > > >> On Aug 4, 2018, at 11:43 AM, Jeff Newmiller <[hidden email]> wrote: >> >> Sometimes a good old for loop performs best, even if it doesn't look sexy: >> >> > > Fair enough, but a vectorized solution beats them all (see below). > > Also, > > [SNIP] > > >> # Charles >> ans1b <- function( a, b ) >> { > > The lapply you put here was from Eric's solution: > >>  xxx <- lapply( seq.int( nrow( A ) ) >>               , function( i ) { >>                    A[ i, ] %o% B[ i, ] >>                 } > > > This is what I had in mind: > > ans1b.corrected <- function( a, b ) { >  yyy <- sapply( seq.int( nrow( a ) ) >                 , function( i ) a[ i, ] %o% b[ i, ] >                 , simplify = "array" >  ) >  zzz <- aperm( yyy, c( 3, 1, 2 ) ) >  zzz > } > > On my system it is slower than a for loop but a lot faster than your benchmark showed with the superfluous code from Eric's solution. > > For speed, a vectorized solution is faster than a for loop by a factor of 3 on my laptop: > > ans0 <- function(A,B){ >  nca <- ncol(A) >  ncb <- ncol(B) >  j.index <- rep(1:nca, times = ncb) >  k.index <- rep(1:nca, each = ncb) >  res <- array(A[, j.index] * B[, k.index], c(nrow(A), nca, ncb)) >  res >  } > > >> microbenchmark( > +   res0 <- ans0(A, B), > +   res1 <- ans1(A, B), > +   res1b <- ans1b.corrected(A, B), > +   res2 <- ans2(A, B), > +   res3 <- ans3(A, B) > + ) > Unit: microseconds >                           expr     min       lq      mean   median       uq     max neval   cld >             res0 <- ans0(A, B)  13.281  18.4960  21.52723  19.9905  23.4750  61.556   100 a >             res1 <- ans1(A, B) 353.121 369.8635 409.77788 381.5840 444.3290 701.256   100     e > res1b <- ans1b.corrected(A, B)  82.816  89.4185 101.52321  95.4275 107.1700 217.357   100   c >             res2 <- ans2(A, B)  49.674  54.4825  61.78278  58.7540  65.5265 172.625   100  b >             res3 <- ans3(A, B) 317.772 342.4220 392.25065 360.4675 436.2125 602.346   100    d >> > > > FWIW, if there are dimnames on A and B, sapply( row names(A), ..., simplify="array") preserves them without further ado. > > Chuck > --------------------------------------------------------------------------- Jeff Newmiller                        The     .....       .....  Go Live... DCN:<[hidden email]>        Basics: ##.#.       ##.#.  Live Go...                                        Live:   OO#.. Dead: OO#..  Playing Research Engineer (Solar/Batteries            O.O#.       #.O#.  with /Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.