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:

  1. 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$.
  2. 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.