TITLE: Solve many linear systems of similar structure
QUESTION [2 upvotes]: Let $\Phi \in R^{N \times n}$, $N > n$, with $\Phi$ having linearly independent columns.
I have to solve a lot of linear systems of the form
$$
(\Phi^\top D_i \Phi) x = \Phi^\top b_i,
$$
where $D_i \in R^{N \times N}$ is a different strictly positive definite diagonal matrix and $b_i \in R^N$ a different RHS for every system but the matrix $\Phi$ stays constant.  
Is there any way to precompute some factorization such that the linear system is fast to solve every iteration? I stumbled upon this "solution": Prefactoring to solve many similar linear systems, but it doesn't seem to work:
>> Phi = randn(100,10);
>> D = diag(randn(100,1));
>> [Q, R] = qr(Phi, 0);
>> norm(inv(Phi'*D*Phi) - inv(R) * Q' * inv(D) * Q * inv(R)')

ans =

    5.8901

REPLY [0 votes]: When $D$ is a positive diagonal matrix, is is very easy to find $\sqrt{D}$ with $\sqrt{D}^2=D$ by simply taking the root of all elements. 
That leaves you with 
$\Phi^T \sqrt{D}^T \sqrt{D} \Phi x = \Phi b$. Using $\Psi = \sqrt{D}\Phi$ you find $\Psi^T \Psi x = \Psi \sqrt{D}^{-1} b$. Lastly you set $\sqrt{D}^{-1}b=r$ and arrive at
$$\Psi^T \Psi x =\Psi^Tr,$$
something that is known as Normal Equations and is used to solve least square problems, in the case of a tall and thin $\Psi$, as you describe it. 
A way better solution to this problem (in terms of numerical stabilty and calculation time) is to solve the least squares problem $\min_x(||\Psi x -r||_2)$. 
If all entries of $D$ are about the same value, you can simply use the QR decomposition of $\Phi$ and solve $$Q_\Phi R_\Phi x-\sqrt{D}^{-1} \sqrt{D}^{-1} b.$$ The error, when using this method is scaled by $||\sqrt{D}^{-1}||_2$, ie the root of the smallest entry of $D$. If that is not too large, you could use the solution $x$ as starting point for iterative refinement.
{"set_name": "stack_exchange", "score": 2, "question_id": 2699557}