Add Golub Riley iteration method (JuliaSmoothOptimizers#123)

Co-authored-by: Dominique <dominique.orban@gmail.com>

Author: MaxenceGollier

src/QRMumps.jl

@@ -36,7 +36,7 @@ include("wrapper/qr_mumps_api.jl")
 
 include("utils.jl")
 
-export qrm_spmat, qrm_spfct,
+export qrm_spmat, qrm_shifted_spmat, qrm_spfct,
     qrm_init, qrm_finalize,
     qrm_spmat_init!, qrm_spmat_init, qrm_spmat_destroy!,
     qrm_spfct_init!, qrm_spfct_init, qrm_spfct_destroy!,
@@ -55,7 +55,10 @@ export qrm_spmat, qrm_spfct,
     qrm_min_norm_semi_normal!, qrm_min_norm_semi_normal,
     qrm_residual_norm!, qrm_residual_norm,
     qrm_residual_orth!, qrm_residual_orth,
-    qrm_refine!, qrm_refine, qrm_set, qrm_get,
+    qrm_refine!, qrm_refine,
+    qrm_shift_spmat, qrm_update_shift_spmat!, 
+    qrm_golub_riley, qrm_golub_riley!,
+    qrm_set, qrm_get,
     qrm_user_permutation!
 
 @doc raw"""
@@ -522,6 +525,67 @@ function qrm_refine! end
 """
 function qrm_refine end
 
+@doc raw"""
+    qrm_shifted_spmat = qrm_shift_spmat(spmat, α = 0)
+
+Given a `spmat` structure representing some sparse matrix `A`, return the block matrix `(A  √α)` as a `qrm_shifted_spmat` type. See `qrm_shifted_spmat` for more information.
+This can be especially useful when `A` is rank deficient, as choosing `α > 0` acts as a regularization.
+
+### Input Arguments
+* `spmat`: the input matrix
+* `α ≥ 0`: the regularization parameter (note the square root in the block matrix representation).
+"""
+function qrm_shift_spmat end 
+
+@doc raw"""
+    qrm_update_shift_spmat!(shifted_spmat, α)
+
+Given a `shifted` block matrix of the form `(A  √α)`, update the parameter `α` in the matrix.
+### Input Arguments
+* `shifted_spmat`: the input matrix of the qrm_shifted_spmat type. See `qrm_shifted_spmat` for more information.
+* `α ≥ 0`: the regularization parameter (note the square root in the block matrix representation).
+"""
+function qrm_update_shift_spmat! end
+
+@doc raw"""
+    x = qrm_golub_riley(spmat, b)
+
+### Input Arguments : 
+* `spmat`: the input matrix of the ill-conditionned system `Ax = b`.
+* `b`: the RHS of the ill-conditionned system. 
+"""
+function qrm_golub_riley end 
+
+@doc raw"""
+    qrm_golub_riley!(shifted_spmat, spfct, x, b, Δx, y; α = ϵm, max_iter = 50, tol = ϵm, transp = 'n')
+
+This method implements the Golub-Riley iteration.
+Given a (possibly ill-conditionned or rank deficient) system `Ax = b` where `A` can have any shape `m×n`, compute `x = A†b = Aᵀ(AAᵀ)†b` where `A†` is the Moore-Penrose pseudoinverse of `A`.
+
+### Input Arguments :
+
+* `shifted_spmat`: a `qrm_shifted_spmat` type representing the matrix `(A  √α)` where `α` is a regularization parameter for the Golub-Riley iteration. See `qrm_shift_spmat` for more information.
+* `spfct`: a sparse factorization object of type `qrm_spfct`.
+* `x`: the approximate solution vector x := A†b = Aᵀy, the size of this vector is n+m, its values are overwritten by this function and it can be uninitialized when the function is called.
+* `b`: the RHS vector of the linear system, the size of this vector is m.
+* `Δx`: an auxiliary vector used to compute the solution, the size of this vector is n+m and can be uninitialized when the function is called.
+* `y`: the vector used to compute y := (AAᵀ)†b, the size of this vector is n and can be uninitialized when the function is called.
+* `Δy`: an auxiliary vector used to compute the solution, the size of this vector is n and can be uninitialized when the function is called.
+
+The definition of the `shifted_spmat` may look odd at first sight. We use such a block matrix as an argument of `qrm_golub_riley!` because the QR factorization of this matrix has the property that
+`RᵀR = AᵀA + αI` which is the system we need to repeatedly solve from for the Golub-Riley iteration. 
+
+The Golub-Riley method works as follows:
+
+* Given A and b, choose a fixed regularization parameter `α > 0` and let x := 0.
+* At each iteration, compute Δx := Aᵀ(AAᵀ + αI)⁻¹ (b - Ax), and update x := x + Δx.
+
+This method has the advantage of only costing one QR-factorization of the `qrm_shifted_spmat` matrix.
+For more details, see 
+A. Dax and L. Eldén, Approximating minimum norm solutions of rank-deficient least squares problems, Numerical Linear Algebra with Applications, 5, pp. 79-99, 1998, DOI 10.1002/(SICI)1099-1506(199803/04)5:2<79::AID-NLA126>3.0.CO;2-4
+"""
+function qrm_golub_riley! end
+
 @doc raw"""
     qrm_set(str, val)
     qrm_set(spfct, str, val)

src/utils.jl

@@ -191,3 +191,107 @@ function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfc
 
   qrm_refine!(spmat, spfct, x, z, Δx, y)
 end
+
+function qrm_shift_spmat(spmat :: qrm_spmat{T}, α :: T = T(0)) where T # Given a m×n matrix A, return the matrix Aₐ = [A √aI].
+  @assert real(α) ≥ 0
+  @assert imag(α) == 0
+
+  m = spmat.mat.m
+  n = spmat.mat.n
+  irn = copy(spmat.irn)
+  jcn = copy(spmat.jcn)
+  val = copy(spmat.val)
+  append!(irn, eltype(irn)(1):eltype(irn)(m))
+  append!(jcn, eltype(irn)(1 + n):eltype(irn)(m + n))
+  append!(val,  sqrt(α)*ones(T,m))
+  shifted_spmat = qrm_spmat_init(m, n + m, irn, jcn, val)
+
+  return qrm_shifted_spmat(shifted_spmat, α)
+end
+
+function qrm_update_shift_spmat!(shifted_spmat :: qrm_shifted_spmat{T}, α :: T) where T
+  shifted_spmat.α = α
+  shifted_spmat.spmat.val[shifted_spmat.spmat.mat.nz - shifted_spmat.spmat.mat.m + 1:end] .= sqrt(α)
+end
+
+function qrm_golub_riley(spmat :: qrm_spmat{T}, b :: AbstractVector{T}; α :: T = T(eps(real(T))), max_iter :: Int = 50, tol :: Real = eps(real(T)), transp :: Char = 'n') where T
+  shifted_spmat = qrm_shift_spmat(spmat, α)
+  spfct = qrm_spfct_init(shifted_spmat.spmat)
+  n = shifted_spmat.spmat.mat.n
+  m = shifted_spmat.spmat.mat.m
+
+  x = similar(b, n)
+  Δx = similar(b, n)
+  y = similar(b)
+  Δy = similar(b) 
+
+  qrm_golub_riley!(
+    shifted_spmat,
+    spfct,
+    x,
+    b,
+    Δx,
+    y,
+    Δy,
+    α = α,
+    max_iter = max_iter,
+    tol = tol,
+    transp = transp
+  )
+  return x[1:n-m]
+end
+
+function qrm_golub_riley!(
+  shifted_spmat :: qrm_shifted_spmat{T},
+  spfct :: qrm_spfct{T},  
+  x :: AbstractVector{T},
+  b :: AbstractVector{T}, 
+  Δx ::AbstractVector{T},
+  y :: AbstractVector{T},
+  Δy :: AbstractVector{T};  
+  α :: T = T(eps(real(T))), 
+  max_iter :: Int = 50, 
+  tol :: Real = eps(real(T)),
+  transp :: Char = 'n'
+  ) where T
+  
+  spmat = shifted_spmat.spmat
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+
+  @assert real(α) ≥ 0
+  @assert imag(α) == 0
+  @assert length(b) == spmat.mat.m
+  @assert length(x) == spmat.mat.n
+  @assert length(y) == spmat.mat.m
+  @assert length(Δy) == spmat.mat.m
+  @assert length(Δx) == spmat.mat.n
+
+  qrm_update_shift_spmat!(shifted_spmat, α)
+  qrm_spfct_init!(spfct, spmat)
+ 
+  qrm_set(spfct, "qrm_keeph", 0)
+  qrm_analyse!(spmat, spfct, transp = transp)
+  qrm_factorize!(spmat, spfct, transp = transp)
+
+  k = 0
+  solved = false
+  x .= T(0)
+  y .= T(0)
+  while k < max_iter && !solved
+
+    qrm_spmat_mv!(spmat, T(1), x, T(0), Δy, transp = ntransp)
+    @. Δy = b - Δy
+
+    qrm_solve!(spfct, Δy, Δx, transp = t)
+    qrm_solve!(spfct, Δx, Δy, transp = 'n')
+
+    qrm_spmat_mv!(spmat, T(1), Δy, T(0), Δx, transp = transp)
+
+    @. x = x + Δx
+    @. y = y + Δy
+    
+    solved = norm(Δx) ≤ tol*norm(x)
+    k = k + 1
+  end
+end

src/wrapper/qr_mumps_common.jl

@@ -35,6 +35,16 @@ mutable struct qrm_spmat{T} <: AbstractSparseMatrix{T, Cint}
   end
 end
 
+@doc raw"""
+This data type represents a "shifted" matrix. When one wants to solve a "regularized" problem of the form `(AᵀA + αI)x = b`,
+one can use a QR factorization of the block matrix `QR = (A  √α)ᵀ` and note that `RᵀR = AᵀA + αI`. This type of problem is especially useful when `A` is poorly conditioned or rank deficient.
+It only contains a `qrm_spmat` matrix representing the above block matrix and the regularization parameter `α`.
+"""
+mutable struct qrm_shifted_spmat{T} <: AbstractSparseMatrix{T, Cint}
+  spmat :: qrm_spmat{T}
+  α     :: T
+end
+
 function Base.cconvert(::Type{Ref{c_spmat{T}}}, spmat :: qrm_spmat{T}) where T
     return spmat.mat
 end
@@ -45,12 +55,18 @@ function Base.unsafe_convert(::Type{Ref{c_spmat{T}}}, spmat :: c_spmat{T}) where
 end
 
 Base.size(spmat :: qrm_spmat) = (spmat.mat.m, spmat.mat.n)
+Base.size(shifted_spmat :: qrm_shifted_spmat) = (shifted_spmat.spmat.mat.m, shifted_spmat.spmat.mat.n)
 SparseArrays.nnz(spmat :: qrm_spmat) = spmat.mat.nz
+SparseArrays.nnz(shifted_spmat :: qrm_shifted_spmat) = shifted_spmat.spmat.mat.nz
 
 function Base.show(io :: IO, ::MIME"text/plain", spmat :: qrm_spmat)
   println(io, "Sparse matrix -- qrm_spmat of size ", size(spmat), " with ", nnz(spmat), " nonzeros.")
 end
 
+function Base.show(io :: IO, ::MIME"text/plain", shifted_spmat :: qrm_shifted_spmat) 
+  println(io, "Sparse matrix -- qrm_spmat of size ", size(shifted_spmat.spmat), " with ", nnz(shifted_spmat.spmat), " nonzeros.")
+end
+
 mutable struct c_spfct{T}
   m        :: Cint
   n        :: Cint

test/test_qrm.jl

@@ -325,7 +325,7 @@ end
       X = qrm_min_norm(spmat, B, seminormal = true)
       R = B - A * X
       @test norm(R) ≤ tol
-      
+
       X = qrm_min_norm_semi_normal(spmat, B)
       R = B - A * X
       @test norm(R) ≤ tol
@@ -337,7 +337,7 @@ end
       qrm_min_norm!(spmat, b, x)
       r = b - A * x
       @test norm(r) ≤ tol
-      
+
       Δx = similar(x)
       y = similar(b)
       qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y)
@@ -669,6 +669,19 @@ end
     x_refined = qrm_refine(spmat, spfct, x, b)
     @test norm(b - A'*(A*x)) ≥ norm(b - A'*(A*x_refined))
 
+    tol = (real(T) == Float32) ? 1e-2 : 1e-9
+    A = sprand(T, m, n, 0.3)
+    k = n ÷ 8
+    A[randperm(n)[1:k],:] .= T(0) # Make A rank deficient.
+    b = A*rand(T, n) # Make sure b is in the range of A
+
+    spmat = qrm_spmat_init(A)
+    x = qrm_golub_riley(spmat, b, transp = transp)
+    B = Matrix(A)
+
+    @test norm(A*x - b) ≤ tol
+    @test norm(x - pinv(B)*b) ≤ tol
+    @test norm(x - A'*pinv(B*B')*b) ≤ tol
   end
 end