Repeated linear least squares problem with large dimensions

[stevengj]

@JADekker

As such, I’d like to save some computational effort and recycle some work if possible. The thing that comes to mind is to compute pinv(X) and to carry out the regression as pinv(X)*y^k. I was wondering if this is the way to go here, or if there is a better approach that I’m missing. I look forward to hearing your thoughts!

The pseudo-inverse (computed by SVD in pinv) is expensive, and I recently read that multiplying by the explicit pinv matrix is not a backwards-stable way of using the pseudo-inverse (Liu & Barnett, 2016, and references therein) . The default method in Julia for least-squares problems via β = X \ y is column-pivoted QR, as explained in this thread: Efficient way of doing linear regression - #33 by stevengj

Because of this, I would tend to use QR = qr(X, ColumnNorm()) and then do β = QR \ y on each step. This re-uses the efficient QR factorization and should be robust and backwards stable.