Solve the corrected semi-normal equations R'Rx=A'b.

x, r = csne(R, A, b) solves the least-squares problem

minimize ||r||_2, where r := b - A*x

using the corrected semi-normal equation approach described by Bjork (1987). Assumes that R is upper triangular.


Add a column to a QR factorization without using Q.

R = qraddcol(A,R,v) adds the m-vector v to the QR factorization of the m-by-n matrix A without using Q. On entry, R is a dense n-by-n upper triangular matrix such that R'*R = A'*A.

R = qraddcol(A,R,v,β) is similar to the above, except that the routine updates the QR factorization of

[A; β I], and R'R = (A'A + β^2I) = R'R.

A should have fewer columns than rows. A and v may be sparse or dense. On exit, R is the (n+1)-by-(n+1) factor corresponding to

Anew = [A V ], Rnew = [R u ], Rnew'Rnew = Anew'Anew. [beta*I ] [ gamma] [ beta]

The new column v is assumed to be nonzero. If A has no columns yet, input A = [], R = [].

15 Jun 2007: First version of QRaddcol.m (without beta).
             Michael Friedlander ( and
             Michael Saunders (
             Where necessary, Ake Bjorck's CSNE
             (Corrected Semi-Normal Equations) method
             is used to improve the accuracy of u and gamma.
             See p143 of Ake Bjork's Least Squares book.
18 Jun 2007: R is now the exact size on entry and exit.
19 Oct 2007: Sparse A, a makes c and u sparse.
             Force them to be dense.
             For dense R we probably should use linsolve,
             which requires c and u to be dense anyway.
04 Aug 2008: QRaddcol2 updates u using du, rather than
             u = R*z as in Ake's book.
             We guess that it might be slightly more accurate,
             but it's hard to tell.  No R*z makes it a little cheaper.
03 Sep 2008: Generalize A to be [A; beta*I] for some scalar beta.
             Update u using du, but keep Ake's version in comments.
29 Dec 2015: Converted to Julia.

Add a row and update a Q-less QR factorization.

qraddrow!(R, a) returns the triangular part of a QR factorization of [A; a], where A = QR for some Q. The argument 'a' should be a row vector.


Delete the k-th column and update a Q-less QR factorization.

R = qrdelcol(R,k) deletes the k-th column of the upper-triangular R and restores it to upper-triangular form. On input, R is an n x n upper-triangular matrix. On output, R is an n-1 x n-1 upper triangle.

18 Jun 2007: First version of QRdelcol.m.
             Michael Friedlander ( and
             Michael Saunders (
             To allow for R being sparse,
             we eliminate the k-th row of the usual
             Hessenberg matrix rather than its subdiagonals,
             as in Reid's Bartel-Golub LU update and also
             the Forrest-Tomlin update.
             (But Matlab will still be pretty inefficient.)
18 Jun 2007: R is now the exact size on entry and exit.
30 Dec 2015: Translate to Julia.