# Efficient evaluation of quadratic form for a Kronecker-decomposable covariance matrix?

Well, it's in the title. I'd like to evaluate $\mathbf{y}'K^{-1}\mathbf{y}$ in the case that $K$ is decomposable into a Kronecker product (e.g. $K = A\otimes B\otimes C$). (EDIT: $K$ is also symmetric and positive definite - it's a covariance matrix in the space of $\mathbf{y}$). The full result I want to calculate is the likelihood:

$$\mathcal{L(y}) := (2\pi)^{-\frac{N}{2}}\lvert K\rvert^{-\frac{1}{2}}\exp\left(-\frac{1}{2}\mathbf{y}'K^{-1}\mathbf{y}\right)$$

What I've tried so far

One nice property of the Kronecker product is that the Cholesky decomposition of the product is the product of Cholesky decompositions (saw this in Schäcke 2013):

$$A \otimes B = (L_A L_A') \otimes (L_BL_B') = (L_A \otimes L_B)(L_A \otimes L_B)'$$ So using this to get the Cholesky decomposition of the overall matrix, I follow the routine suggested in Rasmussen 2006 (Gaussian Processes in Machine Learning): \begin{align} L&:=\text{cholesky}(K)\\ \alpha&:=L' \backslash (L\backslash \mathbf{y}) \\ \log \mathcal{L}(y)&:=-\frac{1}{2}\mathbf{y}'\mathbf{\alpha}-\sum_i L_{ii} - \frac{n}{2}\log(2\pi) \end{align} Where $A\backslash y$ is the solution to the equation $Ax=y$ (which is fast to calculate since L is triangular). I'm finding that this is slow and unstable when $K$ is huge. Could you suggest a simpler way to make the calculation, further taking advantage of the structure of $K$?

EDIT: It would be especially beneficial if some operations can be computed ahead of introducing $\mathbf{y}$. That's because I need to evaluate the function for many values of $\mathbf{y}$ but fixed $K$. Thanks!

• So, your $K$ must be symmetric, positive definite, right? – RTJ Dec 14 '15 at 18:50
• @CTNT yes! In fact it's a covariance matrix, which has those requirements. I will add that now. – cgreen Dec 14 '15 at 19:03