# Tensor EVD

### Definition

Let the \(n\times n\) matrix \(\textbf{K}\) to be a Hadamard product involving two smaller symmetric positive semi-definite (e.g., covariance structures) matrices \(\textbf{K}_1\) and \(\textbf{K}_2\) of dimensions \(n_1\times n_1\) and \(n_2\times n_2\), respectively,

\[ \textbf{K} = (\textbf{Z}_1 \textbf{K}_1 \textbf{Z}'_1) \odot (\textbf{Z}_2 \textbf{K}_2 \textbf{Z}'_2) \]

where \(\textbf{Z}_1\) and \(\textbf{Z}_2\) are incidence matrices mapping from rows (and columns) of the resulting Hadamard \(\textbf{K}\) to rows (and columns) of \(\textbf{K}_1\) and \(\textbf{K}_2\), respectively.

Let the eigenvalue decomposition (EVD) of \(\textbf{K}_1\) and \(\textbf{K}_2\) to be \(\textbf{K}_1 = \textbf{U}_1 \textbf{D}_1 \textbf{U}'_1\) and \(\textbf{K}_2 = \textbf{U}_2 \textbf{D}_2 \textbf{U}'_2\). Using properties of the Hadamard and Kronecker (\(\otimes\)) products, an EVD of the Hadamard product \(\textbf{K}\) can be derived from the EVD of the corresponding Kronecker product between \(\textbf{K}_1\) and \(\textbf{K}_2\) as

\[ \textbf{K} = (\textbf{Z}_1\star\textbf{Z}_2)(\textbf{U}_1\otimes\textbf{U}_2)(\textbf{D}_1\otimes\textbf{D}_2)(\textbf{U}_1\otimes\textbf{U}_2)'(\textbf{Z}_1\star\textbf{Z}_2)' \] This can be simplified as

\[ \textbf{K} = \tilde{\textbf{U}} \textbf{D} \tilde{\textbf{U}}' \]

where \(\textbf{D} = \textbf{D}_1\otimes\textbf{D}_2 = diag(d_1,\dots,d_N)\) is a diagonal matrix containing \(N = n_1\times n_2\) tensor eigenvalues \(d_1 \ge \dots \ge d_N \ge 0\), and \(\tilde{\textbf{U}} = (\textbf{Z}_1\star\textbf{Z}_2)(\textbf{U}_1\otimes\textbf{U}_2) = [\boldsymbol{u}_1,\dots,\boldsymbol{u}_N]\) is matrix containing \(N\) tensor eigenvectors \(\boldsymbol{u}_k\). The term \(\textbf{Z}_1\star\textbf{Z}_2\) is an \(n\times N\) “face-splitting product” matrix (aka “transposed Khatri–Rao product”) of matrices \(\textbf{Z}_1\) and \(\textbf{Z}_2\), obtained as a row-by-row Kronecker product

\[ \textbf{Z}_1\star\textbf{Z}_2 = \begin{pmatrix} \boldsymbol{z}_{11}\otimes\boldsymbol{z}_{12} \\ \boldsymbol{z}_{21}\otimes\boldsymbol{z}_{22}\\ \vdots \\ \boldsymbol{z}_{n1}\otimes\boldsymbol{z}_{n2} \end{pmatrix} \]

with \(\boldsymbol{z}_{i1}\) and \(\boldsymbol{z}_{i2}\) being the \(i^{th}\) row of \(\textbf{Z}_1\) and \(\textbf{Z}_2\), respectively.

### Implementation

The `tensorEVD()` function derives the decomposition of the Hadamard product \(\textbf{K} = \tilde{\textbf{U}} \textbf{D} \tilde{\textbf{U}}'\) formed from two matrices \(\textbf{K}_1\) and \(\textbf{K}_2\). The matrix \(\tilde{\textbf{U}} = [\boldsymbol{u}_1,\dots,\boldsymbol{u}_N]\) is efficiently computed by deriving each column as a Hadamard product (‘\(\odot\)’) using the corresponding \(i_k^{th}\) and \(j_k^{th}\) eigenvectors \(\boldsymbol{u}_{1i_k}\) and \(\boldsymbol{u}_{2j_k}\) of \(\textbf{U}_1\) and \(\textbf{U}_2\), respectively, that form the \(k^{th}\) eigenvector of \(\textbf{U}_1\otimes\textbf{U}_2\), this is

\[ \boldsymbol{u}_k = (\textbf{Z}_1 \boldsymbol{u}_{1i_k}) \odot (\textbf{Z}_2 \boldsymbol{u}_{2j_k}) \]

The terms \(\textbf{Z}_1 \boldsymbol{u}_{1i_k}\) and \(\textbf{Z}_2 \boldsymbol{u}_{2j_k}\) can be obtained by matrix indexing using integer vectors `ID1` and `ID2`, for instance, as `u1ik[ID1]` and `u2jk[ID2]`. Objects `ID1` and `ID2` are vectors of length \(n\) whose entries of are the row (and column) number of \(\textbf{K}_1\) and \(\textbf{K}_2\) that are mapped at each row of \(\textbf{Z}_1\) and \(\textbf{Z}_2\), respectively. Therefore, the `tensorEVD()` can be implemented using as inputs the covariance structure matrices and IDs, e.g.,

    tensorEVD(K1, K2, ID1, ID2)

### Examples

#### 1. Balanced design

When \(\textbf{Z}_1\) and \(\textbf{Z}_2\) are such that the Hadamard \(\textbf{K}\) represent a Kronecker product between \(\textbf{K}_1\) and \(\textbf{K}_2\) (i.e, \(n = n_1\times n_2\) combinations, each element of \(\textbf{K}_1\) crossed with one and only one element of \(\textbf{K}_2\)), then the results of the `tensorEVD()` are exactly the same as those obtained by performing the EVD directly on \(\textbf{K}\) using the `eigen()` function from the ‘base’ R-package.

    # Simulating covariance matrices K1 and K2
    n1 = 10; n2 = 15
    K1 <- crossprod(matrix(rnorm(n1*(n1+10),sd=sqrt(1/n1)), ncol=n1))
    K2 <- crossprod(matrix(rnorm(n2*(n2+10),sd=sqrt(1/n2)), ncol=n2))

    ID1 <- rep(seq(n1), each=n2)
    ID2 <- rep(seq(n2), times=n1)

    # Direct EVD of the Hadamard product
    K <- Hadamard(K1, K2, ID1, ID2, ID1, ID2)   # same as K = K1[ID1,ID1]*K2[ID2,ID2]
    EVD0 <- eigen(K)

    # Tensor EVD using K1 and K2
    EVD <- tensorEVD(K1, K2, ID1, ID2)

    # Eigenvalues and (absolute) eigenvectors are numerically equal
    all.equal(EVD0$values, EVD$values)

    ## [1] TRUE

    all.equal(abs(EVD0$vectors), abs(EVD$vectors))

    ## [1] TRUE

#### 2. Unbalanced design

For unbalanced (i.e., combinations appearing at different frequencies) designs, eigenvectors and eigenvalues are no longer equivalent:

    n <- n1*n2   # size of the Hadamard
    ID1 <- sample(seq(n1), n, replace=TRUE) # random sample of size n of ID1
    ID2 <- sample(seq(n2), n, replace=TRUE) # random sample of size n of ID2

    K <- Hadamard(K1, K2, ID1, ID2, ID1, ID2)
    EVD0 <- eigen(K)
    EVD <- tensorEVD(K1, K2, ID1, ID2)

    all.equal(EVD0$values, EVD$values)

    ## [1] "Mean relative difference: 0.256392"

    all.equal(abs(EVD0$vectors), abs(EVD$vectors))

    ## [1] "Mean relative difference: 1.084041"

However, `tensorEVD()` will still produce a total sum of eigenvalues always equal to the \(trace(\textbf{K})\) and will provide the same approximation \(\textbf{K} = \tilde{\textbf{U}} \textbf{D} \tilde{\textbf{U}}'\)

    # Sum of eigenvalues
    c(sum(EVD0$values), sum(EVD$values), sum(diag(K)))

    ## [1] 477.4273 477.4273 477.4273

    # Approximation for K
    K01 <- EVD0$vectors%*%diag(EVD0$values)%*%t(EVD0$vectors)
    K02 <- EVD$vectors%*%diag(EVD$values)%*%t(EVD$vectors)
    c(all.equal(K,K01), all.equal(K,K02))

    ## [1] TRUE TRUE

### Dimension and rank

The set of eigenvectors with positive eigenvalue are the only ones needed to span the Hadamard matrix \(\textbf{K} = (\textbf{Z}_1 \textbf{K}_1 \textbf{Z}'_1) \odot (\textbf{Z}_2 \textbf{K}_2 \textbf{Z}'_2)\). The size of this set (i.e., the rank of \(\textbf{K}\)) is at most the minimum between \(n\) and \(n_1\times n_2\). The `tensorEVD()` algorithm produces the complete basis containing \(n_1\times n_2\) eigenvectors for the Kronecker matrix product \(\textbf{K}_1\otimes \textbf{K}_2\).

As consequence, `tensorEVD()` can provide more vectors than the ones needed to span \(\textbf{K}\) if the size of the Hadamard product (\(n\)) is considerably smaller than the corresponding Kronecker product (\(n_1\times n_2\)). For example,

    n = n1*n2/2    # size of the Hadamard is half of n1 x n2
    ID1 <- sample(seq(n1), n, replace=TRUE)
    ID2 <- sample(seq(n2), n, replace=TRUE)

    K <- Hadamard(K1, K2, ID1, ID2, ID1, ID2)
    EVD0 <- eigen(K)
    EVD <- tensorEVD(K1, K2, ID1, ID2)

    # Number of eigenvectors with positive eigenvalue
    c(eigen=sum(EVD0$values>1E-10), tensorEVD=sum(EVD$values>1E-10))

    ##     eigen tensorEVD
    ##        63       150

However, when \(n\) is larger than \(n_1\times n_2\), both approaches provide similar number of eigenvectors. For the balanced replicated case, the number of eigenvectors will be the same:

    # Size of the Hadamard is three times n1 x n2
    # Balanced and replicated case
    ID1 <- rep(rep(seq(n1), each=n2), 3)
    ID2 <- rep(rep(seq(n2), times=n1), 3)

    K <- Hadamard(K1, K2, ID1, ID2, ID1, ID2)
    EVD0 <- eigen(K)
    EVD <- tensorEVD(K1, K2, ID1, ID2)

    c(eigen=sum(EVD0$values>1E-10), tensorEVD=sum(EVD$values>1E-10))

    ##     eigen tensorEVD
    ##       150       150

### Proportion of variance explained

Instead of forming all possible eigenvectors, `tensorEVD()` allows to specify a proportion of variance explained (\(0\lt\alpha\leq 1\)) and build only the eigenvectors needed to achieve such proportion of variance. For example,

    alpha <- 0.95
    EVD <- tensorEVD(K1, K2, ID1, ID2, alpha=alpha)
    ncol(EVD$vectors)

    ## [1] 101

    # For the direct EVD
    varexp = cumsum(EVD0$values/sum(EVD0$values))
    index = 1:which.min(abs(varexp-alpha))
    ncol(EVD0$vectors[,index])

    ## [1] 101

### Extras

#### Making dimension names

Row and column names for the eigenvectors of the Hadamard product can be retrieved using the `make.dimnames` argument. Attribute `rownames` of the eigenvectors will be produced by crossing `rownames` of \(\textbf{K}_1\) with those of \(\textbf{K}_2\). Attribute `colnames` will contain the cross between the eigenvector position of the EVD of \(\textbf{K}_1\) with those of \(\textbf{K}_1\) forming each eigenvector of the Hadamard product. For instance,

    dimnames(K1) <- list(paste0("i",seq(n1)), paste0("i",seq(n1)))
    dimnames(K2) <- list(paste0("j",seq(n2)), paste0("j",seq(n2)))

    EVD <- tensorEVD(K1, K2, ID1, ID2, make.dimnames=TRUE)
    EVD$vectors[1:6,1:5]

    ##                1:1           1:2          1:3         2:1          1:4
    ## i1:j1 -0.010991300 -0.0058710318 -0.001176204 -0.04564372  0.005121892
    ## i1:j2 -0.003279555  0.0002818522  0.009558760 -0.01361905 -0.004069084
    ## i1:j3  0.020825652 -0.0153317996  0.011277039  0.08648297 -0.030535547
    ## i1:j4  0.002747478  0.0181337532  0.017851513  0.01140949 -0.023367477
    ## i1:j5  0.007694427 -0.0092708171  0.001547025  0.03195275 -0.010503801
    ## i1:j6  0.017989476 -0.0176471652  0.004568838  0.07470514  0.025347504

#### Pre-calculated EVD

Pre-calculated EVD can be also provided to the function, so, I will not be calculated again.

    EVD2 <- eigen(K2)
    EVD0 <- tensorEVD(K1=K1, EVD2=EVD2, ID1=ID1, ID2=ID2)
    EVD <- tensorEVD(K1=K1, K2=K2, ID1=ID1, ID2=ID2)

    all.equal(EVD0$values, EVD$values)

    ## [1] TRUE

    all.equal(abs(EVD0$vectors), abs(EVD$vectors))

    ## [1] TRUE