Khatri-Rao Matrix Product

Computes Khatri-Rao products for any kind of matrices.

The Khatri-Rao product is a column-wise Kronecker product. Originally introduced by Khatri and Rao (1968), it has many different applications, see Liu and Trenkler (2008) for a survey. Notably, it is used in higher-dimensional tensor decompositions, see Bader and Kolda (2008).

Usage

KhatriRao(X, Y = X, FUN = "*", sparseY = TRUE, make.dimnames = FALSE)

Arguments

X,Y: matrices of with the same number of columns.
FUN: the (name of the) function to be used for the column-wise Kronecker products, see kronecker, defaulting to the usual multiplication.
sparseY: logical specifying if Y should be coerced and treated as sparseMatrix. Set this to FALSE, e.g., to distinguish structural zeros from zero entries.
make.dimnames: logical indicating if the result should inherit dimnames from X and Y in a simple way.

Value

a "CsparseMatrix", say R, the Khatri-Rao product of X (\(n \times k\)) and Y (\(m \times k\)), is of dimension \((n\cdot m) \times k\), where the j-th column, R[,j] is the kronecker product kronecker(X[,j], Y[,j]).

Note

The current implementation is efficient for large sparse matrices.

References

Khatri, C. G., and Rao, C. Radhakrishna (1968) Solutions to Some Functional Equations and Their Applications to Characterization of Probability Distributions. Sankhya: Indian J. Statistics, Series A 30, 167–180.

Bader, Brett W, and Tamara G Kolda (2008) Efficient MATLAB Computations with Sparse and Factored Tensors. SIAM J. Scientific Computing 30, 205–231.

Author

Original by Michael Cysouw, Univ. Marburg; minor tweaks, bug fixes etc, by Martin Maechler.

Examples

## Example with very small matrices:
m <- matrix(1:12,3,4)
d <- diag(1:4)
KhatriRao(m,d)
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>                 
#>  [1,] 1  .  .  .
#>  [2,] .  8  .  .
#>  [3,] .  . 21  .
#>  [4,] .  .  . 40
#>  [5,] 2  .  .  .
#>  [6,] . 10  .  .
#>  [7,] .  . 24  .
#>  [8,] .  .  . 44
#>  [9,] 3  .  .  .
#> [10,] . 12  .  .
#> [11,] .  . 27  .
#> [12,] .  .  . 48
KhatriRao(d,m)
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>                 
#>  [1,] 1  .  .  .
#>  [2,] 2  .  .  .
#>  [3,] 3  .  .  .
#>  [4,] .  8  .  .
#>  [5,] . 10  .  .
#>  [6,] . 12  .  .
#>  [7,] .  . 21  .
#>  [8,] .  . 24  .
#>  [9,] .  . 27  .
#> [10,] .  .  . 40
#> [11,] .  .  . 44
#> [12,] .  .  . 48
dimnames(m) <- list(LETTERS[1:3], letters[1:4])
KhatriRao(m,d, make.dimnames=TRUE)
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>       a  b  c  d
#>  [1,] 1  .  .  .
#>  [2,] .  8  .  .
#>  [3,] .  . 21  .
#>  [4,] .  .  . 40
#>  [5,] 2  .  .  .
#>  [6,] . 10  .  .
#>  [7,] .  . 24  .
#>  [8,] .  .  . 44
#>  [9,] 3  .  .  .
#> [10,] . 12  .  .
#> [11,] .  . 27  .
#> [12,] .  .  . 48
KhatriRao(d,m, make.dimnames=TRUE)
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>                 
#>  [1,] 1  .  .  .
#>  [2,] 2  .  .  .
#>  [3,] 3  .  .  .
#>  [4,] .  8  .  .
#>  [5,] . 10  .  .
#>  [6,] . 12  .  .
#>  [7,] .  . 21  .
#>  [8,] .  . 24  .
#>  [9,] .  . 27  .
#> [10,] .  .  . 40
#> [11,] .  .  . 44
#> [12,] .  .  . 48
dimnames(d) <- list(NULL, paste0("D", 1:4))
KhatriRao(m,d, make.dimnames=TRUE)
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>       a  b  c  d
#>  [1,] 1  .  .  .
#>  [2,] .  8  .  .
#>  [3,] .  . 21  .
#>  [4,] .  .  . 40
#>  [5,] 2  .  .  .
#>  [6,] . 10  .  .
#>  [7,] .  . 24  .
#>  [8,] .  .  . 44
#>  [9,] 3  .  .  .
#> [10,] . 12  .  .
#> [11,] .  . 27  .
#> [12,] .  .  . 48
KhatriRao(d,m, make.dimnames=TRUE)
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>       D1 D2 D3 D4
#>  [1,]  1  .  .  .
#>  [2,]  2  .  .  .
#>  [3,]  3  .  .  .
#>  [4,]  .  8  .  .
#>  [5,]  . 10  .  .
#>  [6,]  . 12  .  .
#>  [7,]  .  . 21  .
#>  [8,]  .  . 24  .
#>  [9,]  .  . 27  .
#> [10,]  .  .  . 40
#> [11,]  .  .  . 44
#> [12,]  .  .  . 48
dimnames(d) <- list(paste0("d", 10*1:4), paste0("D", 1:4))
(Kmd <- KhatriRao(m,d, make.dimnames=TRUE))
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>       a  b  c  d
#> d10:A 1  .  .  .
#> d20:A .  8  .  .
#> d30:A .  . 21  .
#> d40:A .  .  . 40
#> d10:B 2  .  .  .
#> d20:B . 10  .  .
#> d30:B .  . 24  .
#> d40:B .  .  . 44
#> d10:C 3  .  .  .
#> d20:C . 12  .  .
#> d30:C .  . 27  .
#> d40:C .  .  . 48
(Kdm <- KhatriRao(d,m, make.dimnames=TRUE))
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>       D1 D2 D3 D4
#> A:d10  1  .  .  .
#> B:d10  2  .  .  .
#> C:d10  3  .  .  .
#> A:d20  .  8  .  .
#> B:d20  . 10  .  .
#> C:d20  . 12  .  .
#> A:d30  .  . 21  .
#> B:d30  .  . 24  .
#> C:d30  .  . 27  .
#> A:d40  .  .  . 40
#> B:d40  .  .  . 44
#> C:d40  .  .  . 48

nm <- as(m, "nsparseMatrix")
nd <- as(d, "nsparseMatrix")
KhatriRao(nm,nd, make.dimnames=TRUE)
#> 12 x 4 sparse Matrix of class "ngCMatrix"
#>       a b c d
#> d10:A | . . .
#> d20:A . | . .
#> d30:A . . | .
#> d40:A . . . |
#> d10:B | . . .
#> d20:B . | . .
#> d30:B . . | .
#> d40:B . . . |
#> d10:C | . . .
#> d20:C . | . .
#> d30:C . . | .
#> d40:C . . . |
KhatriRao(nd,nm, make.dimnames=TRUE)
#> 12 x 4 sparse Matrix of class "ngCMatrix"
#>       D1 D2 D3 D4
#> A:d10  |  .  .  .
#> B:d10  |  .  .  .
#> C:d10  |  .  .  .
#> A:d20  .  |  .  .
#> B:d20  .  |  .  .
#> C:d20  .  |  .  .
#> A:d30  .  .  |  .
#> B:d30  .  .  |  .
#> C:d30  .  .  |  .
#> A:d40  .  .  .  |
#> B:d40  .  .  .  |
#> C:d40  .  .  .  |

stopifnot(dim(KhatriRao(m,d)) == c(nrow(m)*nrow(d), ncol(d)))
## border cases / checks:
zm <- nm; zm[] <- FALSE # all FALSE matrix
stopifnot(all(K1 <- KhatriRao(nd, zm) == 0), identical(dim(K1), c(12L, 4L)),
          all(K2 <- KhatriRao(zm, nd) == 0), identical(dim(K2), c(12L, 4L)))

d0 <- d; d0[] <- 0; m0 <- Matrix(d0[-1,])
stopifnot(all(K3 <- KhatriRao(d0, m) == 0), identical(dim(K3), dim(Kdm)),
    all(K4 <- KhatriRao(m, d0) == 0), identical(dim(K4), dim(Kmd)),
    all(KhatriRao(d0, d0) == 0), all(KhatriRao(m0, d0) == 0),
    all(KhatriRao(d0, m0) == 0), all(KhatriRao(m0, m0) == 0),
    identical(dimnames(KhatriRao(m, d0, make.dimnames=TRUE)), dimnames(Kmd)))

## a matrix with "structural" and non-structural zeros:
m01 <- new("dgCMatrix", i = c(0L, 2L, 0L, 1L), p = c(0L, 0L, 0L, 2L, 4L),
           Dim = 3:4, x = c(1, 0, 1, 0))
D4 <- Diagonal(4, x=1:4) # "as" d
DU <- Diagonal(4)# unit-diagonal: uplo="U"
(K5  <- KhatriRao( d, m01))
#> 12 x 4 sparse Matrix of class "dgCMatrix"
#>              
#>  [1,] . . . .
#>  [2,] . . . .
#>  [3,] . . . .
#>  [4,] . . . .
#>  [5,] . . . .
#>  [6,] . . . .
#>  [7,] . . 3 .
#>  [8,] . . . .
#>  [9,] . . 0 .
#> [10,] . . . 4
#> [11,] . . . 0
#> [12,] . . . .
K5d  <- KhatriRao( d, m01, sparseY=FALSE)
K5Dd <- KhatriRao(D4, m01, sparseY=FALSE)
K5Ud <- KhatriRao(DU, m01, sparseY=FALSE)
(K6  <- KhatriRao(diag(3),     t(m01)))
#> 12 x 3 sparse Matrix of class "dgCMatrix"
#>            
#>  [1,] . . .
#>  [2,] . . .
#>  [3,] 1 . .
#>  [4,] 1 . .
#>  [5,] . . .
#>  [6,] . . .
#>  [7,] . . .
#>  [8,] . 0 .
#>  [9,] . . .
#> [10,] . . .
#> [11,] . . 0
#> [12,] . . .
K6D  <- KhatriRao(Diagonal(3), t(m01))
K6d  <- KhatriRao(diag(3),     t(m01), sparseY=FALSE)
K6Dd <- KhatriRao(Diagonal(3), t(m01), sparseY=FALSE)
stopifnot(exprs = {
    all(K5 == K5d)
    identical(cbind(c(7L, 10L), c(3L, 4L)),
              which(K5 != 0, arr.ind = TRUE, useNames=FALSE))
    identical(K5d, K5Dd)
    identical(K6, K6D)
    all(K6 == K6d)
    identical(cbind(3:4, 1L),
              which(K6 != 0, arr.ind = TRUE, useNames=FALSE))
    identical(K6d, K6Dd)
})