# Inverse of a diagonal matrix plus a Kronecker product?

Given two matrices $X$ and $Y$, it's easy to take the inverse of their Kronecker product:

$(X\otimes Y)^{-1} = X^{-1}\otimes Y^{-1}$

Now, suppose we have some diagonal matrix $\Lambda$ (or more generally an easily inverted matrix, or one for which we already know the inverse). Is there a closed-form expression or efficient algorithm for computing $(\Lambda + (X\otimes Y))^{-1}$?

• +1 What's the motivation for your question? – draks ... Jan 30 '13 at 21:37
• I have a convex optimization problem I'm trying to solve, and the Hessian of the objective takes the form above. If I can compute the inverse Hessian efficiently I can use Newton's method, which is far preferable to gradient methods. – David Pfau Jan 30 '13 at 22:36
• Did you ever figure out a way to do this that avoids eigenvalue decompositions? I also have come across Hessians with this structure and would be interested in efficient algorithms to solve this problem. – Nick Alger May 9 '16 at 16:20

yes there is.

See equation 5 in http://books.nips.cc/papers/files/nips24/NIPS2011_0443.pdf

Stegle et al. Efficient inference in matrix-variate Gaussian models with iid observation noise

• The matrices corresponding to $X$ and $Y$ in the cited paper are positive semidefinite and $\Gamma$ is a multiple of $I$, but I think what the OP asks for is more general. – user1551 Feb 27 '13 at 8:49
• Yes, its not a complete answer, I agree... but I think the eigen value decomposition as given in the paper will allow at least a partial solution. – David Rohde Mar 4 '13 at 4:40
• To elaborate further... the matrix must be diagonalizable.... maybe another matrix decomposition method could be used in the general case The expression for the inverse is given here en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix I haven't checked, but I think the \Gamma is a multiple of I can be dealt with i.e. this constraint can be at least paritaly removed . So to qualify, I haven't given a full answer here (because, I don't know) - but I think this is the right approach... – David Rohde Mar 4 '13 at 4:42

Let $C := \left\{ c_{i,j} \right\}_{i,j=1}^N$ and $A := \left\{ a_{i,j} \right\}_{i,j=1}^n$ be symmetric matrices. The spectral decompositions of the two matrices read $A = O^T D_A O$ and $C = U^T D_C U$ where $D_A := Diag(\lambda_k)_{k=1}^n$ and $D_C := Diag(\mu_k)_{k=1}^N$ and $O \cdot O^T = O^T \cdot O = 1_n$ and $U \cdot U^T = U^T \cdot U = 1_N$. We use equation 5 from the cited paper in order to compute the resolvent $G_{A \otimes C}(z) := \left(z 1_{ n N} - A \otimes C\right)^{-1}$. We have: \begin{eqnarray} G_{A \otimes C}(z) &=& \left( O^T \otimes U^T\right) \cdot \left( z 1_{ n N} - D_A \otimes D_C \right)^{-1} \cdot \left(O \otimes U \right) \\ &=& \left\{ \sum\limits_{k=1}^n O_{k,i} O_{k,j} U^T \cdot \left( z 1_{N} - \lambda_k D_C\right)^{-1} U \right\}_{i,j=1}^n \\ &=& \sum\limits_{p=0}^\infty \frac{1}{z^{1+p}} A^p \otimes C^p \\ &=& \frac{\sum\limits_{p=0}^{d-1} z^{d-1-p} \sum\limits_{l=d-p}^d (-1)^{d-l} {\mathfrak a}_{d-l} \left(A \otimes C\right)^{p-d+l}}{\sum\limits_{p=0}^d z^{d-p} (-1)^p {\mathfrak a}_p} \end{eqnarray} where $\det( z 1_{n N} - A \otimes C) := \sum\limits_{l=0}^{n N} (-1)^l {\mathfrak a}_l z^{n N-l}$. The first two lines from the top are straightforward. In the third line we expanded the inverse matrix in a series and finally in the fourth line we summed up the infinite series using Cayley-Hamilton's theorem.

Use the following trick:

$$(\Lambda+A)^{-1}=A^{-1}(A^{-1}\Lambda+I)^{-1}$$ where $$A=X\bigotimes Y$$. Then we can decompose $$\Lambda=\Lambda_1\bigotimes \Lambda_2$$ by any feasible $$\Lambda_1$$ and $$\Lambda_2$$ because $$\Lambda$$ is diagonal. Now we have $$A^{-1}\Lambda=X^{-1}\Lambda_1 \bigotimes Y^{-1}\Lambda_2$$ and use the method mentioned by David Rohde