# Inverse of matrix sum, one symmetric PSD and one near-constant diagonal

### Question

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.

### Motivation

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.

### Ideas

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}$$?