I have a matrix $K$ which is the sum of a identity and a Kronecker product of two symmetric matrices as following and I want to calculate the inverse of it $K^{-1}$. \begin{eqnarray} K=\mathbf{I}_{mn}+XX^\top\otimes YY^\top \end{eqnarray} where $X$ and $Y$ are $m\times p$ and $n\times q$ matrices respectively and $m>p,n>q$.

Since the size of Kronecker product of two matrices is the product of their sizes, the matrix to be inverted is very large. Can this calculation be simplified?

Using Woodbury matrix identity, the matrix to be inverted can be reduced from $mn\times mn$ to $pq\times pq$ as \begin{eqnarray} \left(\mathbf{I}+XX^\top\otimes YY^\top\right)^{-1} &=&\left(\mathbf{I}+(X\otimes Y)(X\otimes Y)^\top\right)^{-1}\\\ &=&\mathbf{I}-(X\otimes Y)\left(\mathbf{I}+(X\otimes Y)^\top(X\otimes Y)\right)^{-1}(X\otimes Y)^\top \end{eqnarray} However, in my application, $p$ and $q$ are not much less than $m$ and $n$. Thus inverting a $pq\times pq$ matrix is still a consuming work and I want to find a more simple calculation of $K^{-1}$.

Consider the following problem. If there is not the identity matrix $\mathbf{I}_{m\times n}$, we can invert $XX^\top$ and $YY^\top$ respectively and then calculate the Kronecker product as $$ (XX^\top\otimes YY^\top)^{-1}=(XX^\top)^{-1}\otimes(YY^\top)^{-1} $$ using the property of Kronecker product. However, if the identity item exists, this property cannot be utilized.

Then my question is that if there is any way to decompose the calculation of $K^{-1}$ into inverses of some smaller matrices whose size is linear with $m$ and $n$ but not their products? Or is it possible to extract the Kronecker product $\otimes$ out of the inverse?

If you have any suggestion or idea, please let me know. Thank you very much for your help!

  • 1
    $\begingroup$ If $X$ is $m \times p$ with $m>p$, $(XX^T)^{-1}$ does not exist. $\endgroup$ Commented Jun 11, 2012 at 16:16
  • $\begingroup$ Do you really need the inverse, or do you just want to solve linear systems with it? The latter task is much easier. For the inverse, I don't think there's anything better than what you already suggested with the Woodbury identity. $\endgroup$ Commented Jun 11, 2012 at 16:19
  • $\begingroup$ @Federico Poloni, yes, you are right! $(XX^\top)^{-1}$ really does not exist. Thank you! $\endgroup$
    – ppyang
    Commented Jun 12, 2012 at 2:07
  • $\begingroup$ @Federico Poloni, thank you for your suggestion. What is more, the purpose that I took that example is to show that I want to simplify the calculation utilizing the property of Kronecker product, however, there is some mistake in it. Thank you! $\endgroup$
    – ppyang
    Commented Jun 12, 2012 at 2:12
  • $\begingroup$ @ppyang: So is the thing off or not? Because the original question still makes sense, I think. $\endgroup$ Commented Jun 12, 2012 at 9:08

2 Answers 2


I assume spectral decomposition of a symmetric matrix does not cost too much. Let $PXX^TP=\Lambda_1$, $QYY^TQ=\Lambda_2$ be the spectral decomposition of $XX^T, YY^T$, respectively. Here $\Lambda_1, \Lambda_2$ are diagonal and $P, Q$ are orthogonal. Then $$K^{-1}=(P\otimes Q)^T(I+\Lambda_1\otimes \Lambda_2)^{-1}(P\otimes Q).$$

$(I+\Lambda_1\otimes \Lambda_2)^{-1}$ can be read directly.

  • $\begingroup$ You provide an excellent idea that reduce the most complex calculation of an $mn\times mn$ matrix into a $m\times m$ matrix and a $n\times n$ matrix, where the inverse is taken on only a diagonal. I think it solve my problem. Thank you very much! $\endgroup$
    – ppyang
    Commented Jun 17, 2012 at 11:38

Excellent answer by Betrand. But it's easy to verify that the final form should be $K^{-1} = (P\otimes Q)(I + \Lambda_1 \otimes \Lambda_2)^{-1} (P \otimes Q)^T$.

Here is some R code:

A = matrix(c(3, 2, 2, 3), ncol = 2)
B = matrix(c(2, 2, 3, 4), ncol = 2)
AA = A %*% t(A)
BB = B %*% t(B)
AA_decomp = eigen(AA)
P = AA_decomp$vectors
LAMB_A = diag(AA_decomp$values)
P %*% LAMB_A %*% t(P)
BB_decomp = eigen(BB)
Q = BB_decomp$vectors
LAMB_B = diag(BB_decomp$values)
Q %*% LAMB_B %*% t(Q)
K = diag(dim(A)[1] * dim(B)[1]) + kronecker(AA, BB)
kronecker(P, Q) %*% solve(diag(dim(A)[1] * dim(B)[1]) + kronecker(LAMB_A, LAMB_B)) %*% t(kronecker(P, Q))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.