Kronecker products on Array objects
Array-kronecker-methods.RdThe S4Arrays package implements kronecker() methods
for Array objects that work out-of-the-box on Array
derivatives that support [ and *.
Note that kronecker is a generic function defined in the
methods package but documented in the base package.
See ?base::kronecker.
Usage
# S4 method for class 'Array,ANY'
kronecker(X, Y, FUN="*", make.dimnames=FALSE, ...)
# S4 method for class 'ANY,Array'
kronecker(X, Y, FUN="*", make.dimnames=FALSE, ...)
# S4 method for class 'Array,Array'
kronecker(X, Y, FUN="*", make.dimnames=FALSE, ...)
## The workhorse behind the three above methods.
kronecker2(X, Y, FUN="*", make.dimnames=FALSE, ...)Arguments
- X, Y
Array-like objects. Alternatively
Xand/orYcan be vectors, in which case they are converted to 1D-arrays withas.array().Note that
XandYare expected to have the same number of dimensions. However, when they don't, ineffective dimensions (i.e. dimensions with an extent of 1) are added to the object with less dimensions.For the S4 methods, at least one of
XorYmust be an Array derivative.- FUN, make.dimnames, ...
See
?base::kroneckerfor a description of these arguments.
Details
The three kronecker() methods listed above
delegate to kronecker2(), a re-implementation of
base::kronecker().
kronecker2() is semantically equivalent to
base::kronecker. However, unlike the latter which
calls as.array() on X and Y internally,
the former operates natively on the input objects regardless
of their internal representations, as long as they are array-like
objects that support [ (single-bracket subsetting) and *.
In particular, when X and Y use the same internal
representations, the returned object will also use that representation.
In other words, the output object will have the same class as the
input objects (endomorphism).
The endomorphism property is particularly important when the input
objects are sparse (e.g. SVT_SparseArray objects
from the SparseArray package) or when they use an on-disk
representation (e.g. DelayedArray objects from
the DelayedArray package). For example, if X and Y
are DelayedArray objects, the returned object is
another DelayedArray object. Also in that case,
calling kronecker2() is virtually instantaneous because all the
operations that the function performs internally on X and
Y by the are delayed!
See also
base::kroneckerin the base package.The "Kronecker product" page on Wikipedia: https://en.wikipedia.org/wiki/Kronecker_product
The Array class.
SparseArray objects implemented in the SparseArray package.
DelayedArray objects implemented in the DelayedArray package.
TENxMatrix objects implemented in the HDF5Array package.
Examples
## ---------------------------------------------------------------------
## SIMPLE kronecker2() EXAMPLES
## ---------------------------------------------------------------------
m1 <- matrix(1:10, nrow=2) # 2 x 5 matrix
m2 <- matrix(101:106, nrow=3) # 3 x 2 matrix
kronecker2(m1, m2) # 6 x 10 matrix
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 101 104 303 312 505 520 707 728 909 936
#> [2,] 102 105 306 315 510 525 714 735 918 945
#> [3,] 103 106 309 318 515 530 721 742 927 954
#> [4,] 202 208 404 416 606 624 808 832 1010 1040
#> [5,] 204 210 408 420 612 630 816 840 1020 1050
#> [6,] 206 212 412 424 618 636 824 848 1030 1060
a1 <- array(1:16, dim=c(4, 2, 2))
a2 <- array(1:30, dim=c(3, 5, 2))
kronecker2(a1, a2) # 12 x 10 x 4 array
#> , , 1
#>
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 1 4 7 10 13 5 20 35 50 65
#> [2,] 2 5 8 11 14 10 25 40 55 70
#> [3,] 3 6 9 12 15 15 30 45 60 75
#> [4,] 2 8 14 20 26 6 24 42 60 78
#> [5,] 4 10 16 22 28 12 30 48 66 84
#> [6,] 6 12 18 24 30 18 36 54 72 90
#> [7,] 3 12 21 30 39 7 28 49 70 91
#> [8,] 6 15 24 33 42 14 35 56 77 98
#> [9,] 9 18 27 36 45 21 42 63 84 105
#> [10,] 4 16 28 40 52 8 32 56 80 104
#> [11,] 8 20 32 44 56 16 40 64 88 112
#> [12,] 12 24 36 48 60 24 48 72 96 120
#>
#> , , 2
#>
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 16 19 22 25 28 80 95 110 125 140
#> [2,] 17 20 23 26 29 85 100 115 130 145
#> [3,] 18 21 24 27 30 90 105 120 135 150
#> [4,] 32 38 44 50 56 96 114 132 150 168
#> [5,] 34 40 46 52 58 102 120 138 156 174
#> [6,] 36 42 48 54 60 108 126 144 162 180
#> [7,] 48 57 66 75 84 112 133 154 175 196
#> [8,] 51 60 69 78 87 119 140 161 182 203
#> [9,] 54 63 72 81 90 126 147 168 189 210
#> [10,] 64 76 88 100 112 128 152 176 200 224
#> [11,] 68 80 92 104 116 136 160 184 208 232
#> [12,] 72 84 96 108 120 144 168 192 216 240
#>
#> , , 3
#>
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 9 36 63 90 117 13 52 91 130 169
#> [2,] 18 45 72 99 126 26 65 104 143 182
#> [3,] 27 54 81 108 135 39 78 117 156 195
#> [4,] 10 40 70 100 130 14 56 98 140 182
#> [5,] 20 50 80 110 140 28 70 112 154 196
#> [6,] 30 60 90 120 150 42 84 126 168 210
#> [7,] 11 44 77 110 143 15 60 105 150 195
#> [8,] 22 55 88 121 154 30 75 120 165 210
#> [9,] 33 66 99 132 165 45 90 135 180 225
#> [10,] 12 48 84 120 156 16 64 112 160 208
#> [11,] 24 60 96 132 168 32 80 128 176 224
#> [12,] 36 72 108 144 180 48 96 144 192 240
#>
#> , , 4
#>
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 144 171 198 225 252 208 247 286 325 364
#> [2,] 153 180 207 234 261 221 260 299 338 377
#> [3,] 162 189 216 243 270 234 273 312 351 390
#> [4,] 160 190 220 250 280 224 266 308 350 392
#> [5,] 170 200 230 260 290 238 280 322 364 406
#> [6,] 180 210 240 270 300 252 294 336 378 420
#> [7,] 176 209 242 275 308 240 285 330 375 420
#> [8,] 187 220 253 286 319 255 300 345 390 435
#> [9,] 198 231 264 297 330 270 315 360 405 450
#> [10,] 192 228 264 300 336 256 304 352 400 448
#> [11,] 204 240 276 312 348 272 320 368 416 464
#> [12,] 216 252 288 324 360 288 336 384 432 480
#>
## The Kronecker product is **not** commutative:
m1 <- matrix(LETTERS[1:10], nrow=2)
m2 <- matrix(letters[1:6], nrow=3)
kronecker2(m1, m2, paste, sep="*")
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] "A*a" "A*d" "C*a" "C*d" "E*a" "E*d" "G*a" "G*d" "I*a" "I*d"
#> [2,] "A*b" "A*e" "C*b" "C*e" "E*b" "E*e" "G*b" "G*e" "I*b" "I*e"
#> [3,] "A*c" "A*f" "C*c" "C*f" "E*c" "E*f" "G*c" "G*f" "I*c" "I*f"
#> [4,] "B*a" "B*d" "D*a" "D*d" "F*a" "F*d" "H*a" "H*d" "J*a" "J*d"
#> [5,] "B*b" "B*e" "D*b" "D*e" "F*b" "F*e" "H*b" "H*e" "J*b" "J*e"
#> [6,] "B*c" "B*f" "D*c" "D*f" "F*c" "F*f" "H*c" "H*f" "J*c" "J*f"
kronecker2(m2, m1, paste, sep="*")
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] "a*A" "a*C" "a*E" "a*G" "a*I" "d*A" "d*C" "d*E" "d*G" "d*I"
#> [2,] "a*B" "a*D" "a*F" "a*H" "a*J" "d*B" "d*D" "d*F" "d*H" "d*J"
#> [3,] "b*A" "b*C" "b*E" "b*G" "b*I" "e*A" "e*C" "e*E" "e*G" "e*I"
#> [4,] "b*B" "b*D" "b*F" "b*H" "b*J" "e*B" "e*D" "e*F" "e*H" "e*J"
#> [5,] "c*A" "c*C" "c*E" "c*G" "c*I" "f*A" "f*C" "f*E" "f*G" "f*I"
#> [6,] "c*B" "c*D" "c*F" "c*H" "c*J" "f*B" "f*D" "f*F" "f*H" "f*J"
## Sanity checks:
stopifnot(
identical(kronecker2(m1, m2, paste0), kronecker(m1, m2, paste0)),
identical(kronecker2(m2, m1, paste0), kronecker(m2, m1, paste0))
)
## ---------------------------------------------------------------------
## USING kronecker() ON Array DERIVATIVES
## ---------------------------------------------------------------------
## The user should typically avoid direct calls to kronecker2() and
## stick to kronecker(). Because this is a generic function, it will
## dispatch to the appropriate method based on the classes of the input
## objects. If one of them is an Array derivative, kronecker2() will
## be called thanks to the methods defined in the S4Arrays package and
## listed in the Usage section above.
## With SparseMatrix objects (Array derivatives):
library(SparseArray)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#>
#> Attaching package: ‘MatrixGenerics’
#> The following objects are masked from ‘package:matrixStats’:
#>
#> colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#> colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#> colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#> colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#> colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#> colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#> colWeightedMeans, colWeightedMedians, colWeightedSds,
#> colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#> rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#> rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#> rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#> rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#> rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#> rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#> rowWeightedSds, rowWeightedVars
sm1 <- poissonSparseMatrix(300, 15, density=0.25)
sm2 <- poissonSparseMatrix(80, 500, density=0.15)
kronecker2(sm1, sm2) # 24000 x 7500 SparseMatrix object
#> <24000 x 7500 SparseMatrix> of type "integer" [nzcount=6576528 (3.7%)]:
#> [,1] [,2] [,3] [,4] ... [,7497] [,7498] [,7499] [,7500]
#> [1,] 0 0 0 0 . 0 0 0 1
#> [2,] 0 0 0 0 . 0 0 0 0
#> [3,] 0 0 0 0 . 0 0 0 0
#> [4,] 0 0 0 0 . 0 0 0 0
#> [5,] 0 0 0 0 . 0 0 0 0
#> ... . . . . . . . . .
#> [23996,] 0 0 0 0 . 0 0 0 0
#> [23997,] 0 0 0 0 . 0 0 0 0
#> [23998,] 0 0 0 0 . 0 0 0 0
#> [23999,] 0 0 0 0 . 0 0 0 0
#> [24000,] 0 0 0 0 . 0 0 0 0
## With TENxMatrix objects (DelayedArray derivatives, therefore also
## Array derivatives):
library(HDF5Array)
#> Loading required package: DelayedArray
#>
#> Attaching package: ‘DelayedArray’
#> The following objects are masked from ‘package:base’:
#>
#> apply, scale, sweep
#> Loading required package: h5mread
#> Loading required package: rhdf5
#>
#> Attaching package: ‘h5mread’
#> The following object is masked from ‘package:rhdf5’:
#>
#> h5ls
M1 <- writeTENxMatrix(sm1) # 300 x 15 TENxMatrix object
M2 <- writeTENxMatrix(sm2) # 80 x 500 TENxMatrix object
K <- kronecker2(M1, M2) # instantaneous! (all operations are delayed)
showtree(K) # show delayed operations details
#> 24000x7500 integer: DelayedMatrix object
#> └─ 24000x7500 integer: N-ary iso op
#> ├─ 24000x7500 integer: Subset
#> │ └─ 300x15 integer, sparse: [seed] TENxMatrixSeed object
#> └─ 24000x7500 integer: Subset
#> └─ 80x500 integer, sparse: [seed] TENxMatrixSeed object
## ---------------------------------------------------------------------
## VERIFYING THE MIXED-PRODUCT PROPERTY (JUST FOR FUN!)
## ---------------------------------------------------------------------
## See https://en.wikipedia.org/wiki/Kronecker_product for details
## about "The mixed-product property".
## We verify the property on 4 random matrices:
A <- matrix(runif(1000), ncol=40)
B <- matrix(runif(800), ncol=100)
C <- matrix(runif(600), nrow=40)
D <- matrix(runif(5000), nrow=100)
kAB <- kronecker2(A, B)
kCD <- kronecker2(C, D)
kAB_x_kCD <- kAB %*% kCD
A_x_C <- A %*% C
B_x_D <- B %*% D
stopifnot(all.equal(kAB_x_kCD, kronecker2(A_x_C, B_x_D)))
## The mixed-product property also for the element-wise product
## (Hadamard product). We verify this on 4 random arrays:
A <- array(1:60, dim=5:3)
B <- array(101:180, dim=c(2,10,4))
C <- array(runif(60), dim=5:3)
D <- array(runif(80), dim=c(2,10,4))
kAB <- kronecker2(A, B)
kCD <- kronecker2(C, D)
kAB_o_kCD <- kAB * kCD
A_o_C <- A * C
B_o_D <- B * D
stopifnot(all.equal(kAB_o_kCD, kronecker2(A_o_C, B_o_D)))