Update the "Q-less" QR factorization of a matrix. Routines are provided for adding and deleting columns, adding rows, and solving associated linear least-squares problems.
The least-squares solver uses Björck's corrected semi-normal equation (CSNE) approach [1] with one step of iterative refinement. Using double precision, the method should be stable for matrices A
with condition number up to 10^8
.
import Pkg
Pkg.add("QRupdate")
Build the "Q-less" QR factorization of A
one column at a time:
m, n = 100, 50
A = randn(m,0)
R = Array{Float64, 2}(undef, 0, 0)
for i in 1:n
a = randn(m)
R = qraddcol(A, R, a)
A = [A a]
end
Minimize the least-squares residual ||Ax - b||₂
using the computed R
:
b = randn(m)
x, r = csne(R, A, b)
Delete a random column and compute new R
:
n = size(A,2)
k = rand(1:n)
A = A[:, 1:n .!= k]
R = qrdelcol(R, k)
Minimize the least-squares residual ||Ax - b||₂
using the computed R
:
x, r = csne(R, A, b)
Add a row to A
:
n = size(A,2)
a = randn(n)' # must be row vector
A = [A; a]
R = qraddrow(R, a)
Minimize the least-squares residual ||Ax - b||₂
using the computed R
:
b = [b; randn()]
x, r = csne(R, A, b)
These operations consider that you preivously allocate the matrices involved. A depth argument is required when adding columns.
m, n = 100, 4
# Allocate matrices
A = zeros(m,n)
R = zeros(n,n)
# Then add/remove
a = randn(m)
current_R_size = 0
qraddcol!(A,R,a,current_R_size)
current_R_size = 1
a = randn(m)
qraddcol!(A,R,a,current_R_size)
qrdelcol!(A,R,2)
[1] Björck, A. (1996). Numerical methods for least squares problems. SIAM.
β
).
u
and γ
. See p143 of Ake Bjork's Least Squares book.R
is now the exact size on entry and exit.A
, a makes c
and u
sparse. Force them to be dense.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.A
to be [A; β*I]
for some scalar β
. Update u
using du
, but keep Ake's version in comments.