# Re-using computations in several least squares problems

I have $$K$$ least squares problems of the form $$Y_k = X_k\beta_k$$ for $$k = 1, \dots, n$$. If the matrix $$X_k$$ is the same for each index $$k$$, we can rewrite the problem as $$Y_k = X \beta_k$$. How can I speed up least squares by knowing this additional information? The values of $$Y_k$$ change.

Of course if $$K=1$$, it's just a a single linear regression problem. But how can I reuse my work? Is this possible? Thanks

Reusing your work is a good idea. To start, let's review how the lm function in R calculates least squares estimates. We'll assume throughout that $$X \in \mathbb{R}^{n \times p}$$ has full rank and that $$n \geq p$$.

The function calculates that $$X=QR$$, where $$Q \in \mathbb{R}^{n \times p}$$ has orthonormal columns and $$R \in \mathbb{R}^{p \times p}$$ is upper triangular. Then, it uses this decomposition to quickly calculate $$(X^TX)^{-1}X^T Y_k = R^{-1} Q^T Y_k$$ by taking the matrix-vector product $$Q^T Y_k$$ and then backsolving. Notice that the decomposition of $$X$$ does not rely on the response $$Y_k$$. This suggests a general strategy: we can calculate a something involving $$X$$ once and then reuse it for each response vector.

If you have a very large number of response vectors, it will be advantageous to actually compute $$(X^TX)^{-1}X^T \in \mathbb{R}^{p \times n}$$. Then each calculation of a least squares coefficient will only involve a matrix-vector product and avoid backsolving.

Otherwise, simply storing the $$QR$$ decomposition of $$X$$ will greatly speed up your least squares run time.

• This is a great answer--I wonder why you deleted it? Am I overlooking some subtlety or limitation?
– whuber
Commented Apr 15, 2022 at 13:12
• Thanks @whuber . I wanted to add some flop counts and fix the error saying the R uses LU.
– Ben
Commented Apr 15, 2022 at 21:15