In some algorithm I need to compute $a=(X^tX)^{-1}X^ty$ and $b=(X^tX)^{-1} z$ in each step, where $X$ is a $n \times p$ non-square matrix ($n \geq p$, $p$ is increased in each step) and $y,z$ are some appropriate vectors, but I'd like to do so in an efficient way.
What I came up with so far:
- To compute just $a$ it is worth noting that $a$ is a least squares solution, i.e. $a$ minimizes $\Vert Xa-y\Vert_2^2$, so I would have calculated $a$ using the QR decomposition of $X$ (cost is $O(np^2)$, but now I also need to compute $b$.
- So another naive approach would be computing $W := X^tX$ and $v = X^ty$, and then computing $a= W^{-1} v$, $b= W^{-1} z$. This approach has the advantage that we could use e.g. a Cholesky decomposition of $W$ (cost $O(p^3)$), but the backdraws are that we need to compute $W$ first (cost $O(np^2)$) and that $X^tX$ has a worse condition number than $X$.
But is there a more efficient way to compute $a$ and $b$?
With $X=QR$ in the small variant where $R$ is square, you get $X^tX=R^tR$ so that $a=R^{-1}Q^ty$ and $b=R^{-1}(R^t)^{-1}z$ where you have only triangular systems to solve.