ENH: QR updating for penalized L2 estimation

[josef-pkt]

Just an idea while looking at the code for TheilGLS

(x'x + a S) beta = x'y
can be computed using data augmentation S_half = S^{-1/2}, vstack([x, S_half])
We use this in PIRLS in GAM

linalg:

we need `(x'x + a S)^{-1} = (x'x + a S_half' S_half)^{-1}

If S_half is low rank, i.e. has few rows, then this is a low rank updating of the inverse matrix, which is cheaper to compute than the full matrix inverse #6265.
Similar would apply if we use the QR decomposition of x and update the decomposition to vstack([x, S_half]) (but there we wouldn't have the inverse yet.

The advantage would be that we can do the updating for different pen_weights which are included in S_half

aside:
we have code for S_half in singular case with reduced number of rows.
tools.linalg.matrix_sqrt (adding with GAM PR)