A slightly unorthodox matrix product.

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

A slightly unorthodox matrix product.

Rolf Turner

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/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.

Eric Berger
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-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.

Berry, Charles


> 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.
Reply | Threaded
Open this post in threaded view
|

Re: A slightly unorthodox matrix product.

Jeff Newmiller
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-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.

Berry, Charles


> 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-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.

Jeff Newmiller
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-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.

Berry, Charles


> On Aug 4, 2018, at 12:59 PM, Jeff Newmiller <[hidden email]> wrote:
>
...

> Slick work on the vectorizing, but for the future reference it was slightly buggy:
>

Thanks for catching that!

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.
Reply | Threaded
Open this post in threaded view
|

Re: [FORGED] Re: A slightly unorthodox matrix product.

Rolf Turner


Thanks to Eric Berger, Jeff Newmiller and Chuck Berry for their help
with this.  Thanks especially to Eric for catching the bugs in my
proposed solution.  I can *never* keep the indexing straight in
multidimensional arrays!  It's tough being a <expletive deleted>-wit!

I wouldn't have figured out the "ans0b" solution in a million years.

Thanks again.

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.