For any two matrices \(\textbf{A} = \{a_{ij}\}\) of dimensions \(m\times n\) and \(\textbf{B} = \{b_{ij}\}\) of dimensions \(p\times q\), the direct Kronecker product between them is the matrix of dimensions \(mp\times nq\) defined as the block matrix
\[ \textbf{A}\otimes\textbf{B} = \{a_{ij}\textbf{B}\} \]
This can be computed using the kronecker()
function from
the ‘base’ R-package. Alternatively, the Kronecker()
function from the ‘tensorEVD’ R-package can be used.
Simple scalar multiplication. Let \(a\) be a scalar and \(\textbf{B}\) any matrix. Then, computing their Kronecker product is the same as multiplying \(\textbf{B}\) by the scalar:
a <- 10
( B <- matrix(1:6, ncol=2) )
## [,1] [,2]
## [1,] 1 4
## [2,] 2 5
## [3,] 3 6
Kronecker(a, B)
## [,1] [,2]
## [1,] 10 40
## [2,] 20 50
## [3,] 30 60
# In this case, should be equal to Kronecker(B, a)
Kronecker(B, a)
## [,1] [,2]
## [1,] 10 40
## [2,] 20 50
## [3,] 30 60
Block diagonal matrix. A block diagonal matrix can be formed using the Kronecker product with a diagonal matrix
D <- diag(1, 2)
Kronecker(D, B)
## [,1] [,2] [,3] [,4]
## [1,] 1 4 0 0
## [2,] 2 5 0 0
## [3,] 3 6 0 0
## [4,] 0 0 1 4
## [5,] 0 0 2 5
## [6,] 0 0 3 6
# Is not equal to
Kronecker(B, D) # this a 'striped' matrix
## [,1] [,2] [,3] [,4]
## [1,] 1 0 4 0
## [2,] 0 1 0 4
## [3,] 2 0 5 0
## [4,] 0 2 0 5
## [5,] 3 0 6 0
## [6,] 0 3 0 6
Outer product. Let \(\textbf{a}\) and \(\textbf{b}\) be any two vectors. Then, the Kronecker product between \(\textbf{a}\) and the transpose of \(\textbf{b}\) form an outer product:
a <- c(1,2,3)
b <- c(4,5)
Kronecker(a, t(b))
## [,1] [,2]
## [1,] 4 5
## [2,] 8 10
## [3,] 12 15
# Should be equal to the product a b'
tcrossprod(a, b)
## [,1] [,2]
## [1,] 4 5
## [2,] 8 10
## [3,] 12 15
Here we compare tensorEVD’s Kronecker()
function with
the kronecker()
function from the ‘base’ R-package in the
computation of a Kronecker product of two general matrices \(\textbf{A}\) and \(\textbf{B}\).
# Simulating matrices A and B
m = 20; n = 20
p = 40; q = 30
A <- matrix(rnorm(m*n), ncol=n)
B <- matrix(rnorm(p*q), ncol=q)
# Making the Kronecker product
K1 <- kronecker(A, B)
K2 <- Kronecker(A, B)
# Should be equal
all.equal(K1,K2)
## [1] TRUE
Here we compare these methods in terms of computational speed in
scenarios with small and large matrices. We included in the benchmark
the kronecker.prod()
function from the ‘fastmatrix’
R-package. The following benchmark was performed using the code provided
in this script
run on a Linux environment based on the following system settings:
Selecting specific rows and columns from the Kronecker can be done by pre- and post- multiplication with incidence matrices, for instance
\[ \textbf{R}(\textbf{A}\otimes\textbf{B})\textbf{C}' \]
where \(\textbf{R}\) is an incidence matrix for rows and \(\textbf{C}\) is an incidence matrix for columns.
This sub-matrix can be obtained by indexing rows or columns using,
for instance, integer vectors rows
and cols
as
Kronecker(A, B)[rows,cols]
. However, this approach computes
first the Kronecker product then makes the sub-setting.
This approach can be inefficient if a relatively small number of rows
and columns are to be selected. The Kronecker()
function
can derive this sub-matrix directly from \(\textbf{A}\) and \(\textbf{B}\) on the fly without forming the
whole Kronecker product using rows
and cols
as
arguments. For example,
dm <- c(nrow(A)*nrow(B), ncol(A)*ncol(B)) # dimension of the Kronecker
# Subsetting a matrix with 30% of rows/columns
rows <- sample(seq(dm[1]), 0.3*dm[1])
cols <- sample(seq(dm[2]), 0.3*dm[2])
K1 <- Kronecker(A, B)[rows,cols]
K2 <- Kronecker(A, B, rows=rows, cols=cols)
dim(K1) # small size
## [1] 240 180
all.equal(K1, K2)
## [1] TRUE
Here we show some benchmark results on the time performance of
subsetting a Kronecker matrix using
Kronecker(A, B, rows, cols)
and
Kronecker(A, B)[rows,cols]
, for small (30% of the
rows/columns) and large (300% of the rows/columns) sub-matrices
scenarios, from a Kronecker product matrix of dimension \(5000\times 5000\). The code to perform this
benchmark is provided in this script.
If \(\textbf{A}\) is a matrix and \(\textbf{B}\) is a scalar, the resulting Kronecker will be of the same dimensions as \(\textbf{A}\); therefore, we could overwrite the result on the memory occupied by \(\textbf{A}\).
Usually, assigning an output to an object will occupy a different memory address than inputs:
B <- rnorm(1) # B is a scalar
K <- Kronecker(A, B)
c(K=pryr::address(K), A=pryr::address(A))
## K A
## "0x7fd4a1e12200" "0x7fd4a1277c00"
The parameter inplace
can be used to overwrite the
output at the same address as the input:
A <- Kronecker(A, B, inplace=TRUE)
c(A=pryr::address(A)) # output address remain unchanged
## A
## "0x7fd4a1277c00"
all.equal(K, A) # contains the desired result
## [1] TRUE
Row and column names for the Kronecker product can be retrieved using
the make.dimnames
argument. Attribute dimnames
of the Kronecker will be produced by crossing rownames
and
colnames
of input \(\textbf{A}\) with those of input \(\textbf{B}\). For instance,
A <- matrix(1:9, ncol=3)
B <- matrix(10*(1:4), ncol=2)
dimnames(A) <- list(paste("week",1:3), month.abb[1:3])
dimnames(B) <- list(c("office","home"), LETTERS[1:2])
Kronecker(A, B, make.dimnames=TRUE)
## Jan:A Jan:B Feb:A Feb:B Mar:A Mar:B
## week 1:office 10 30 40 120 70 210
## week 1:home 20 40 80 160 140 280
## week 2:office 20 60 50 150 80 240
## week 2:home 40 80 100 200 160 320
## week 3:office 30 90 60 180 90 270
## week 3:home 60 120 120 240 180 360