How can split the calculation of a real matrix inverse $(S + D)^{-1}$ when I know that $S$ is symmetric and PSD and $D$ diagonal with only a handful of unique values (=diag$(a,a...a,b...b,c...c,c)$)? Ideally the split would be in two parts, where the one dependent on $D$ would be significantly faster to calculate than the rest.


I have a weighted linear-least-squares system with regularisation: $$W A~x \approx Wc~,$$ with ~$X$0,000 equations and 50 unknowns. There are $X$ (typically 2-5) blocks of equations, roughly 10,000 each, and the weights are the same within a block (i.e. $W$ is diagonal with each block being $W_i=w_i I$).

Now I need to solve this system many times over and over, with different $c$. I can precompute the pseudoinverse $(A^TWA+\lambda I)^{-1}A^TW$ ($\lambda$ is Tikhonov regularisation) and then just multiply each new $c$ with that, which is very fast. BUT when I need to change the weights (and I do), I have to redo the matrix inversion, which is rather slow (>200 ms).

So far I've got (using the PosDef identity [The Matrix Cookbook, Eq. (185)] and some shuffling around) to: $$(A^TWA+\lambda I)^{-1}A^TW = A^T(AA^T+\lambda W^{-1})^{-1}~.$$ Now I'd like to split this to one part which can be slow to calculate but independent of $W$, and the remainder which needs to be calculated much faster.


This is similar to Inverse of sum of two marices, one being diagonal and other unitary. and Inverse of the sum of a symmetric and diagonal matrices but not quite the same. I feel like the fact that $$\lambda W^{-1} = \text{diag}(a',a',...a',b',...b',c',...c',c')$$ should be somehow helpful. Maybe I could use the Woodbury identity as stated in https://en.wikipedia.org/wiki/Woodbury_matrix_identity#Applications - if I could find such $U$ and $V$ as to go from a tiny $X\times X$ matrix $D'=\text{diag}(a',b',c')$ to the full diagonal $D$, that could be a great speed-up with precomputed $AA^T$? And how would it change if I set $\lambda = 0$, i.e. I'd want to decompose only $(A^TWA)^{-1}$?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.