From 2ed56b72897408aa36edff97a72f2cb2a23b5d58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=A1s=20Barnafi?= Date: Mon, 24 Jun 2024 16:57:25 -0400 Subject: [PATCH 1/6] Added in-place column delete --- src/QRupdate.jl | 42 +++++++++++++++++++++++++++++++++++++++++- test/runtests.jl | 10 ++++++++++ 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/src/QRupdate.jl b/src/QRupdate.jl index ad16218..245801a 100644 --- a/src/QRupdate.jl +++ b/src/QRupdate.jl @@ -2,7 +2,7 @@ module QRupdate using LinearAlgebra -export qraddcol, qraddrow, qrdelcol, csne +export qraddcol, qraddrow, qrdelcol, qrdelcol!, csne """ Add a column to a QR factorization without using Q. @@ -143,6 +143,46 @@ function qrdelcol(R::AbstractMatrix{T}, k::Int) where {T} return R end +""" +This function is identical to the previous one, but instead leaves R +with a column of zeros. This is useful to avoid copying the matrix. +""" +function qrdelcol!(R::AbstractMatrix{T}, k::Int) where {T} + + m = size(R,1) + n = size(R,2) # Should have m=n+1 + + # Shift columns. This is apparently faster than copying views. + for j in (k+1):n, i in 1:m + @inbounds R[i,j-1] = R[i, j] + end + for i in 1:m + @inbounds R[i,n] = 0.0 + end + + + for j in k:(n-1) # Forward sweep to reduce k-th row to zeros + G, y = givens(R[j+1,j], R[k,j], 1, 2) + R[j+1,j] = y + if j Date: Mon, 24 Jun 2024 18:05:03 -0400 Subject: [PATCH 2/6] Added in-place column update and its test --- src/QRupdate.jl | 111 ++++++++++++++++++++++++++++++++++++++++++++++- test/runtests.jl | 23 ++++++++++ 2 files changed, 133 insertions(+), 1 deletion(-) diff --git a/src/QRupdate.jl b/src/QRupdate.jl index 245801a..01f9077 100644 --- a/src/QRupdate.jl +++ b/src/QRupdate.jl @@ -2,7 +2,41 @@ module QRupdate using LinearAlgebra -export qraddcol, qraddrow, qrdelcol, qrdelcol!, csne +export qraddcol, qraddcol!, qraddrow, qrdelcol, qrdelcol!, csne + +""" +Auxiliary function used to solve fully allocated but incomplete R matrices. +See documentation of qraddcol! . +""" +function solveR!(R::Matrix, b::Vector, sol::Vector, realSize::Int64) + # Note: R is upper triangular + @inbounds sol[realSize] = b[realSize] / R[realSize, realSize] + for i in (realSize-1):-1:1 + @inbounds sol[i] = b[i] + for j in realSize:-1:(i+1) + @inbounds sol[i] = sol[i] - R[i,j] * sol[j] + end + @inbounds sol[i] = sol[i] / R[i,i] + end +end + +""" +Auxiliary function used to solve transpose of fully allocated but incomplete R matrices. +See documentation of qraddcol! . +""" +function solveRT!(R::Matrix, b::Vector, sol::Vector, realSize::Int64) + # Note: R is upper triangular + @inbounds sol[1] = b[1] / R[1, 1] + for i in 2:realSize + @inbounds sol[i] = b[i] + for j in 1:(u) + @inbounds sol[i] = sol[i] - R[i,j] * sol[j] + end + @inbounds sol[i] = sol[i] / R[i,i] + end +end + + """ Add a column to a QR factorization without using Q. @@ -83,6 +117,81 @@ function qraddcol(A::AbstractMatrix{T}, Rin::AbstractMatrix{T}, a::Vector{T}, β return Rout end +""" +This function is identical to the previous one, but it does in-place calculations on R. It requires +knowing which is the last and new column of A. The logic is that 'A' is allocated at the beginning +of an iterative procedure, and thus does not require further allocations: + +It 0 --> It 1 --> It 2 +A = [0 0 0] A = [a1 0 0] A = [a1 a2 0] + +and so on. This yields that R is + +It 0 --> It 1 --> It 2 +R = [0 0 0 R = [r11 0 0 R = [r11 r12 0 + 0 0 0 0 0 0 0 r22 0 + 0 0 0] 0 0 0] 0 0 0] +""" +function qraddcol!(A::AbstractMatrix{T}, R::AbstractMatrix{T}, a::Vector{T}, N::Int64, β::T = zero(T)) where {T} + + m, n = size(A) + anorm = norm(a) + anorm2 = anorm^2 + β2 = β^2 + if β != 0 + anorm2 = anorm2 + β2 + anorm = sqrt(anorm2) + end + + if n == 0 + #return reshape([anorm], 1, 1) + R[1,1] = anorm + end + + c = zeros(n) + u = zeros(n) + du = zeros(n) + mul!(c, A', a) #c = A'a + solveRT!(R, c, u, N) #u = R'\c + unorm2 = norm(u)^2 + d2 = anorm2 - unorm2 + + z = zeros(n) + dz = zeros(n) + r = zeros(m) + + if d2 > anorm2 #DISABLE 0.01*anorm2 # Cheap case: γ is not too small + γ = sqrt(d2) + else + solveR!(R, u, z, N) #z = R\u # First approximate solution to min ||Az - a|| + #r = a - A*z + mul!(r, A, z) # r = Az + axpby!(1.0, a, -1, r) + mul!(c, A', r) #c = A'r + if β != 0 + axpy!(-β2, z, c) #c = c - β2*z + end + solveRT!(R, c, du, N) #du = R'\c + solveR!(R, du, dz, N) #dz = R\du + axpy!(1, dz, z) #z += dz # Refine z + # u = R*z # Original: Bjork's version. + axpy!(1, du, u) #u += du # Modification: Refine u + #r = a - A*z + mul!(r, A, z) + axpby!(1, a, -1, r) + γ = norm(r) # Safe computation (we know gamma >= 0). + if !iszero(β) + γ = sqrt(γ^2 + β2*norm(z)^2 + β2) + end + end + + # This seems to be faster than concatenation, ie: + # [ Rin u + # zeros(1,n) γ ] + R[:,N+1] .= u + R[N+1,N+1] = γ +end + """ Add a row and update a Q-less QR factorization. diff --git a/test/runtests.jl b/test/runtests.jl index 1dae50f..04f92dc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -68,3 +68,26 @@ end @test norm(R) - norm(R2) < 1e-14 end + +@testset "qraddcol!" begin + m = 100 + A = randn(m, 0) + b = randn(m) + R = zeros(0, 0) + + # Do two steps + a1 = randn(m) + R = qraddcol(A, R, a1) + A = [A a1] + a2 = randn(m) + R = qraddcol(A, R, a2) + A = [A a2] + + # Now compute two steps of inplace + Rin = zeros(2,2) + Ain = zeros(m, 2) + qraddcol!(Ain, Rin, a1, 0) + qraddcol!(Ain, Rin, a2, 1) + @test norm(R) - norm(Rin) < 1e-10 +end + From fe6fb74bfc3d221685476c9ae7e286d3f2d792b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=A1s=20Barnafi?= Date: Mon, 24 Jun 2024 19:46:33 -0400 Subject: [PATCH 3/6] Improved tests to be similar to existing ones --- src/QRupdate.jl | 7 ++++--- test/runtests.jl | 36 ++++++++++++++++++------------------ 2 files changed, 22 insertions(+), 21 deletions(-) diff --git a/src/QRupdate.jl b/src/QRupdate.jl index 01f9077..b3e6a0a 100644 --- a/src/QRupdate.jl +++ b/src/QRupdate.jl @@ -29,7 +29,7 @@ function solveRT!(R::Matrix, b::Vector, sol::Vector, realSize::Int64) @inbounds sol[1] = b[1] / R[1, 1] for i in 2:realSize @inbounds sol[i] = b[i] - for j in 1:(u) + for j in 1:(i-1) @inbounds sol[i] = sol[i] - R[i,j] * sol[j] end @inbounds sol[i] = sol[i] / R[i,i] @@ -143,9 +143,10 @@ function qraddcol!(A::AbstractMatrix{T}, R::AbstractMatrix{T}, a::Vector{T}, N:: anorm = sqrt(anorm2) end - if n == 0 + if N == 0 #return reshape([anorm], 1, 1) R[1,1] = anorm + return end c = zeros(n) @@ -284,7 +285,7 @@ function qrdelcol!(R::AbstractMatrix{T}, k::Int) where {T} # Shift k-th row. We skipped the removed column. for j in k:(n-1) - for i in k:(j+1) + for i in k:j @inbounds R[i,j] = R[i+1, j] end R[j+1,j] = 0 diff --git a/test/runtests.jl b/test/runtests.jl index 04f92dc..84a7b83 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -63,31 +63,31 @@ end m = 100 A = randn(m,m) Q, R = qr(A) - R2 = qrdelcol(R, 10) - qrdelcol!(R, 10) - @test norm(R) - norm(R2) < 1e-14 + Qin, Rin = qr(A) + for i in 100:-1:1 + k = rand(1:i) + A = A[:,1:i .!= k] + R = qrdelcol(R, k) + qrdelcol!(Rin, k) + @test norm(R) - norm(Rin) < 1e-14 + end end @testset "qraddcol!" begin m = 100 A = randn(m, 0) - b = randn(m) R = zeros(0, 0) - - # Do two steps - a1 = randn(m) - R = qraddcol(A, R, a1) - A = [A a1] - a2 = randn(m) - R = qraddcol(A, R, a2) - A = [A a2] + Rin = zeros(m,m) + Ain = zeros(m, m) + for i in 1:m + a = randn(m) + R = qraddcol(A, R, a) + A = [A a] + qraddcol!(Ain, Rin, a, i-1) + @test norm(R) - norm(Rin) < 1e-10 + end - # Now compute two steps of inplace - Rin = zeros(2,2) - Ain = zeros(m, 2) - qraddcol!(Ain, Rin, a1, 0) - qraddcol!(Ain, Rin, a2, 1) - @test norm(R) - norm(Rin) < 1e-10 end + From f207e4285aaa543036627f4cab1da1e250e6ac22 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=A1s=20Barnafi?= Date: Mon, 24 Jun 2024 19:56:36 -0400 Subject: [PATCH 4/6] Added column deletion to A for consistency and expanded the README --- README.md | 21 +++++++++++++++++++++ src/QRupdate.jl | 8 +++++--- test/runtests.jl | 3 ++- 3 files changed, 28 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 8dc59ab..751a993 100644 --- a/README.md +++ b/README.md @@ -75,6 +75,27 @@ b = [b; randn()] x, r = csne(R, A, b) ``` +### In-place operations + +These operations consider that you preivously allocate the matrices involved. A depth argument is required when adding columns. + +```julia +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) +``` + ## Reference [1] Björck, A. (1996). Numerical methods for least squares problems. SIAM. diff --git a/src/QRupdate.jl b/src/QRupdate.jl index b3e6a0a..70339e5 100644 --- a/src/QRupdate.jl +++ b/src/QRupdate.jl @@ -257,17 +257,19 @@ end This function is identical to the previous one, but instead leaves R with a column of zeros. This is useful to avoid copying the matrix. """ -function qrdelcol!(R::AbstractMatrix{T}, k::Int) where {T} +function qrdelcol!(A::AbstractMatrix{T},R::AbstractMatrix{T}, k::Int) where {T} - m = size(R,1) - n = size(R,2) # Should have m=n+1 + m = size(A,1) + n = size(A,2) # Should have m=n+1 # Shift columns. This is apparently faster than copying views. for j in (k+1):n, i in 1:m @inbounds R[i,j-1] = R[i, j] + @inbounds A[i,j-1] = A[i, j] end for i in 1:m @inbounds R[i,n] = 0.0 + @inbounds A[i,n] = 0.0 end diff --git a/test/runtests.jl b/test/runtests.jl index 84a7b83..cdd8d29 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,13 +62,14 @@ end @testset "qrdelcol!" begin m = 100 A = randn(m,m) + Ain = copy(A) Q, R = qr(A) Qin, Rin = qr(A) for i in 100:-1:1 k = rand(1:i) A = A[:,1:i .!= k] R = qrdelcol(R, k) - qrdelcol!(Rin, k) + qrdelcol!(Ain, Rin, k) @test norm(R) - norm(Rin) < 1e-14 end end From 6117828f9fac34aef37ecba8a11eee9f83e89a85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=A1s=20Barnafi?= Date: Mon, 24 Jun 2024 23:48:39 -0400 Subject: [PATCH 5/6] Removed redundant operations and fixed some nasty bugs --- src/QRupdate.jl | 69 ++++++++++++++++++++++++++++++++++-------------- test/runtests.jl | 30 ++++++++++++++++++--- 2 files changed, 76 insertions(+), 23 deletions(-) diff --git a/src/QRupdate.jl b/src/QRupdate.jl index 70339e5..27b81e8 100644 --- a/src/QRupdate.jl +++ b/src/QRupdate.jl @@ -8,7 +8,7 @@ export qraddcol, qraddcol!, qraddrow, qrdelcol, qrdelcol!, csne Auxiliary function used to solve fully allocated but incomplete R matrices. See documentation of qraddcol! . """ -function solveR!(R::Matrix, b::Vector, sol::Vector, realSize::Int64) +function solveR!(R::AbstractMatrix{T}, b::Vector{T}, sol::Vector{T}, realSize::Int64) where {T} # Note: R is upper triangular @inbounds sol[realSize] = b[realSize] / R[realSize, realSize] for i in (realSize-1):-1:1 @@ -24,13 +24,13 @@ end Auxiliary function used to solve transpose of fully allocated but incomplete R matrices. See documentation of qraddcol! . """ -function solveRT!(R::Matrix, b::Vector, sol::Vector, realSize::Int64) +function solveRT!(R::AbstractMatrix{T}, b::Vector{T}, sol::Vector{T}, realSize::Int64) where {T} # Note: R is upper triangular @inbounds sol[1] = b[1] / R[1, 1] for i in 2:realSize @inbounds sol[i] = b[i] for j in 1:(i-1) - @inbounds sol[i] = sol[i] - R[i,j] * sol[j] + @inbounds sol[i] = sol[i] - R[j,i] * sol[j] end @inbounds sol[i] = sol[i] / R[i,i] end @@ -135,6 +135,12 @@ R = [0 0 0 R = [r11 0 0 R = [r11 r12 0 function qraddcol!(A::AbstractMatrix{T}, R::AbstractMatrix{T}, a::Vector{T}, N::Int64, β::T = zero(T)) where {T} m, n = size(A) + + # First add vector to A + for i in 1:m + @inbounds A[i,N+1] = a[i] + end + anorm = norm(a) anorm2 = anorm^2 β2 = β^2 @@ -149,26 +155,42 @@ function qraddcol!(A::AbstractMatrix{T}, R::AbstractMatrix{T}, a::Vector{T}, N:: return end - c = zeros(n) - u = zeros(n) - du = zeros(n) - mul!(c, A', a) #c = A'a + c = zeros(N) + u = zeros(N) + du = zeros(N) + + for i in 1:N #c = A'a + for j in 1:m + @inbounds c[i] += A[j,i] * a[j] + end + end solveRT!(R, c, u, N) #u = R'\c unorm2 = norm(u)^2 d2 = anorm2 - unorm2 - z = zeros(n) - dz = zeros(n) + z = zeros(N) + dz = zeros(N) r = zeros(m) - if d2 > anorm2 #DISABLE 0.01*anorm2 # Cheap case: γ is not too small + if d2 > anorm2 γ = sqrt(d2) else - solveR!(R, u, z, N) #z = R\u # First approximate solution to min ||Az - a|| - #r = a - A*z - mul!(r, A, z) # r = Az - axpby!(1.0, a, -1, r) - mul!(c, A', r) #c = A'r + solveR!(R, u, z, N) #z = R\u + #mul!(r, A, z, -1, 1) #r = a - A*z + for i in 1:m + @inbounds r[i] = a[i] + for j in 1:N + @inbounds r[i] -= A[i,j] * z[j] + end + end + #mul!(c, A', r) #c = A'r + c[:] .= 0.0 + for i in 1:N + for j in 1:m + @inbounds c[i] += A[j,i] * r[j] + end + end + if β != 0 axpy!(-β2, z, c) #c = c - β2*z end @@ -177,9 +199,14 @@ function qraddcol!(A::AbstractMatrix{T}, R::AbstractMatrix{T}, a::Vector{T}, N:: axpy!(1, dz, z) #z += dz # Refine z # u = R*z # Original: Bjork's version. axpy!(1, du, u) #u += du # Modification: Refine u - #r = a - A*z - mul!(r, A, z) - axpby!(1, a, -1, r) + #r = a - A*z + for i in 1:m + @inbounds r[i] = a[i] + for j in 1:N + @inbounds r[i] -= A[i,j] * z[j] + end + end + γ = norm(r) # Safe computation (we know gamma >= 0). if !iszero(β) γ = sqrt(γ^2 + β2*norm(z)^2 + β2) @@ -189,7 +216,9 @@ function qraddcol!(A::AbstractMatrix{T}, R::AbstractMatrix{T}, a::Vector{T}, N:: # This seems to be faster than concatenation, ie: # [ Rin u # zeros(1,n) γ ] - R[:,N+1] .= u + for i in 1:N + @inbounds R[i, N+1] = u[i] + end R[N+1,N+1] = γ end @@ -277,7 +306,7 @@ function qrdelcol!(A::AbstractMatrix{T},R::AbstractMatrix{T}, k::Int) where {T} G, y = givens(R[j+1,j], R[k,j], 1, 2) R[j+1,j] = y if j Date: Tue, 25 Jun 2024 00:17:48 -0400 Subject: [PATCH 6/6] . --- src/QRupdate.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/QRupdate.jl b/src/QRupdate.jl index 27b81e8..d0e627c 100644 --- a/src/QRupdate.jl +++ b/src/QRupdate.jl @@ -231,7 +231,7 @@ vector. function qraddrow(R::AbstractMatrix{T}, a::AbstractMatrix{T}) where {T} n = size(R,1) - @inbounds @simd for k in 1:n + @inbounds for k in 1:n G, r = givens( R[k,k], a[k], 1, 2 ) B = G * [ reshape(R[k,k:n], 1, n-k+1) reshape(a[:,k:n], 1, n-k+1) ] @@ -271,7 +271,7 @@ function qrdelcol(R::AbstractMatrix{T}, k::Int) where {T} G, y = givens(R[j+1,j], R[k,j], 1, 2) R[j+1,j] = y if j