Efficient recalculation of weighted least squares regression when weights change

Question

I am performing weighted least squares regression as described on wiki: WLS

I need to solve this equation: $B= (t(X)WX)^{-1}*t(X)Wy$

I use SVD to find: $(t(X)WX)^{-1}$ and store it in a matrix. In addition I store the matrix $H= (t(X)WX)^{-1}*t(X)W$ and simply do the following for any new value of y: B= Hy. This way I am able to save the cost of repeating SVD and matrix multiplications as y's change.

W is a diagonal matrix and generally does not change. However sometimes I change one or two elements on the diagonal in the W matrix. In that case I need to do SVD again and recalc the H matrix. This is clearly slow and time consuming.

My question is: If I know what changed in W and nothing changes in X is there a more efficient method to recalculate (t(X)WX)^-1?

Or put differently is there an efficient analytic method to find B given that only diagonal elements in W can change by a known amount?

Answer 1

There is such a method, in the case that the inverse you compute is a true inverse and not a generalised inverse (ie none of the singular values are 0). However some caution in using this is recommended. If you were doing your sums in infinite precision, all would be well. With finite precision, and particularly with nearly singular problems -- if some of the singular values are very large -- these formulae could result in loss of precision.

I'll call inverse you store C. If you add d (which can be positive or negative) to the m'th weight, then the modified C matrix, C~ say, and the modified H, H~, can be computed like this:

(' denotes transpose, and e_m is row the vector that's all 0, except the m'th slot is 1)

Let

c = the m'th column of H, divided by the original m'th weight    
a = m'th row of the data matrix X    
f = e_m - a*H    
gamma = 1/d + a*c

(so c is a column vector, while a and f are row vectors)

Then

C~ = C - c*c'/gamma
H~ = H + c*f/gamma

If you want to find the new B, B~ say, for a given y, it can be calculated via:

r = y[m] - a*B
B~ = B + (r/gamma) * c

The derivation of these formulae is straightforward, but tedious, matrix algebra. The matrix inversion lemma comes in handy.