Add option to store the previous Jacobian in ShiftedComposite

src/ShiftedProximalOperators.jl

@@ -79,6 +79,7 @@ function shift!(ψ::ShiftedProximableFunction, shift::AbstractVector{R}) where {
 end
 
 function shift!(ψ::ShiftedCompositeProximableFunction, shift::AbstractVector{R}) where {R <: Real}
+  !isnothing(ψ.A_prev) && (ψ.A_prev.vals .= ψ.A.vals) # Update previous Jacobian if necessary
   ψ.c!(ψ.b, shift)
   ψ.J!(ψ.A, shift)
   return ψ

src/compositeNormL2.jl

@@ -2,7 +2,7 @@
 export CompositeNormL2
 
 @doc raw"""
-    CompositeNormL2(λ, c!, J!, A, b)
+    CompositeNormL2(λ, c!, J!, A, b; store_previous_jacobian::Bool = false)
 
 Returns function `c` composed with the `ℓ₂` norm:
 ```math
@@ -22,6 +22,8 @@ such that `J` is the Jacobian of `c`. It is expected that `m ≤ n`.
     c!(b <: AbstractVector{Real}, xk <: AbstractVector{Real})
     J!(A <: AbstractSparseMatrixCOO{Real, Integer}, xk <: AbstractVector{Real})
 ```
+Moreover, if you want shifted instances of the operator to store the previous Jacobian on each shift, you can specify `store_previous_jacobian = true`.
+This is particularly useful for quasi-Newton updates in the context of constrained optimization.
 """
 mutable struct CompositeNormL2{
   T <: Real,
@@ -35,19 +37,21 @@ mutable struct CompositeNormL2{
   J!::F1
   A::M
   b::V
+  store_previous_jacobian::Bool
 
   function CompositeNormL2(
     λ::T,
     c!::Function,
     J!::Function,
     A::AbstractMatrix{T},
-    b::AbstractVector{T},
+    b::AbstractVector{T};
+    store_previous_jacobian::Bool = false
   ) where {T <: Real}
     λ > 0 || error("CompositeNormL2: λ should be positive")
     length(b) == size(A, 1) || error(
       "Composite Norm L2: Wrong input dimensions, the length of c(x) should be the same as the number of rows of J(x)",
     )
-    new{T, typeof(c!), typeof(J!), typeof(A), typeof(b)}(NormL2(λ), c!, J!, A, b)
+    new{T, typeof(c!), typeof(J!), typeof(A), typeof(b)}(NormL2(λ), c!, J!, A, b, store_previous_jacobian)
   end
 end
 

src/shiftedCompositeNormL2.jl

@@ -1,6 +1,6 @@
 export ShiftedCompositeNormL2
 @doc raw"""
-    ShiftedCompositeNormL2(h, c!, J!, A, b)
+    ShiftedCompositeNormL2(h, c!, J!, A, b; store_previous_jacobian::Bool = false)
 
 Returns the shift of a function `c` composed with the `ℓ₂` norm (see CompositeNormL2.jl).
 Here, `c` is linearized i.e, `c(x + s) ≈ c(x) + J(x)s`. 
@@ -21,18 +21,23 @@ such that `J` is the Jacobian of `c`. It is expected that `m ≤ n`.
 c!(b <: AbstractVector{Real}, xk <: AbstractVector{Real})
 J!(A <: AbstractSparseMatrixCOO{Real, Integer}, xk <: AbstractVector{Real})
 ```
+Moreover, if you want shifted instances of the operator to store the previous Jacobian on each shift, you can specify `store_previous_jacobian = true`.
+In this case, each time a shift is performed, the previous Jacobian is stored in the `A_prev` field.
+This is particularly useful for quasi-Newton updates in the context of constrained optimization.
 """
 mutable struct ShiftedCompositeNormL2{
   T <: Real,
   F0 <: Function,
   F1 <: Function,
   M <: AbstractMatrix{T},
+  N <: Union{Nothing, M},
   V <: AbstractVector{T},
 } <: ShiftedCompositeProximableFunction
   h::NormL2{T}
   c!::F0
   J!::F1
   A::M
+  A_prev::N # (Optional) can be used to store the previous Jacobian, useful for quasi-Newton approximations
   shifted_spmat::qrm_shifted_spmat{T}
   spfct::qrm_spfct{T}
   b::V
@@ -41,12 +46,14 @@ mutable struct ShiftedCompositeNormL2{
   dq::V # Preallocated vector to refine the q solution.
   p::V  # Preallocated vector used to compute s(α)ᵀ∇s(α) for the secular equation.
   dp::V # Preallocated vector used to refine the p vector.
+  full_row_rank::Bool # Boolean that tells whether A has full row rank or not. Is updated on each call to `prox!`
   function ShiftedCompositeNormL2(
     λ::T,
     c!::Function,
     J!::Function,
     A::AbstractMatrix{T},
-    b::AbstractVector{T},
+    b::AbstractVector{T};
+    store_previous_jacobian::Bool = false
   ) where {T <: Real}
     p = similar(b, A.n + A.m)
     dp = similar(b, A.n + A.m)
@@ -59,15 +66,18 @@ mutable struct ShiftedCompositeNormL2{
       )
     end
 
+    A_prev = store_previous_jacobian ? copy(A) : nothing
+
     spmat = qrm_spmat_init(A; sym = false)
     shifted_spmat = qrm_shift_spmat(spmat)
     spfct = qrm_spfct_init(spmat)
 
-    new{T, typeof(c!), typeof(J!), typeof(A), typeof(b)}(
+    new{T, typeof(c!), typeof(J!), typeof(A), typeof(A_prev), typeof(b)}(
       NormL2(λ),
       c!,
       J!,
       A,
+      A_prev,
       shifted_spmat,
       spfct,
       b,
@@ -76,6 +86,7 @@ mutable struct ShiftedCompositeNormL2{
       dq,
       p,
       dp,
+      false
     )
   end
 end
@@ -94,7 +105,7 @@ shifted(
   ψ.c!(b, xk)
   A = similar(ψ.A)
   ψ.J!(A, xk)
-  ShiftedCompositeNormL2(ψ.h.lambda, ψ.c!, ψ.J!, A, b)
+  ShiftedCompositeNormL2(ψ.h.lambda, ψ.c!, ψ.J!, A, b, store_previous_jacobian = ψ.store_previous_jacobian)
 end
 
 fun_name(ψ::ShiftedCompositeNormL2) = "shifted `ℓ₂` norm"
@@ -103,13 +114,14 @@ fun_params(ψ::ShiftedCompositeNormL2) = "c(xk) = $(ψ.b)\n" * " "^14 * "J(xk) =
 
 function prox!(
   y::AbstractVector{T},
-  ψ::ShiftedCompositeNormL2{T, F0, F1, M, V},
+  ψ::ShiftedCompositeNormL2{T, F0, F1, M, N, V},
   q::AbstractVector{T},
   ν::T;
   max_iter = 10,
   atol = eps(T)^0.3,
   max_time = 180.0,
-) where {T <: Real, F0 <: Function, F1 <: Function, M <: AbstractMatrix{T}, V <: AbstractVector{T}}
+) where {T <: Real, F0 <: Function, F1 <: Function, M <: AbstractMatrix{T}, N <: Union{Nothing, M}, V <: AbstractVector{T}}
+  @assert ν > zero(T)
   start_time = time()
   θ = T(0.8)
   α = zero(T)
@@ -134,8 +146,8 @@ function prox!(
   _obj_dot_grad!(spmat, spfct, ψ.p, ψ.q, ψ.g, ψ.dq)
 
   # Check full-rankness
-  full_row_rank = (qrm_get(spfct, "qrm_rd_num") == 0)
-  if !full_row_rank
+  ψ.full_row_rank = (qrm_get(spfct, "qrm_rd_num") == 0)
+  if !ψ.full_row_rank
     # QRMumps cannot factorize rank-deficient matrices; use the Golub-Riley iteration instead
     α = αmin
     qrm_golub_riley!(

test/runtests.jl

@@ -48,7 +48,8 @@ for (op, composite_op, shifted_op) ∈
     y = similar(x)
     ν = 0.1056
     prox!(y, ϕ, x, ν)
-
+    @test ϕ.full_row_rank == true
+
     if "$op" == "NormL2"
       y_true = [0.24545429, 0.75250248, -0.66619752, 1.19372286]
       norm = Op(1.0)
@@ -63,6 +64,37 @@ for (op, composite_op, shifted_op) ∈
     @test ϕ.A == SparseMatrixCOO(Float64[2 0 0 -1; 0 1 1 0])
     @test ϕ(ones(Float64, 4)) ==
           h([1.0, 2.0] + SparseMatrixCOO(Float64[2 0 0 -1; 0 1 1 0])*ones(Float64, 4))
+
+    # test store previous Jacobian option
+
+    function c_store_prev!(z, x)
+      z[1] = 2*x[1]^2 - x[4]
+      z[2] = x[2] + x[3]
+    end
+    function J_store_prev!(z, x)
+      z.vals .= Float64[4.0*x[1], 1.0, 1.0, -1.0]
+    end
+    ψ_store_prev = CompositeOp(λ, c_store_prev!, J_store_prev!, A, b, store_previous_jacobian = true)
+    @test ψ_store_prev.store_previous_jacobian == true
+
+    xk = [1.0, 0.0, 0.0, 0.0]
+    ϕ_store_prev = shifted(ψ_store_prev, xk)
+    @test !isnothing(ϕ_store_prev.A_prev)
+
+    @test all(ϕ_store_prev.A_prev .== ϕ_store_prev.A)  # A_prev should be a copied version of A at this point
+    @test all(ϕ_store_prev.A .== SparseMatrixCOO([4.0 0.0 0.0 -1.0; 0.0 1.0 1.0 0.0]))
+
+    xk = [2.0, 0.0, 0.0, 0.0]
+    shift!(ϕ_store_prev, xk)
+    @test all(ϕ_store_prev.A_prev .== SparseMatrixCOO([4.0 0.0 0.0 -1.0; 0.0 1.0 1.0 0.0]))
+    @test all(ϕ_store_prev.A .== SparseMatrixCOO([8.0 0.0 0.0 -1.0; 0.0 1.0 1.0 0.0]))
+
+    # test reshifting
+    xk = [0.5, 0.0, 0.0, 0.0]
+    shift!(ϕ_store_prev, xk)
+    @test all(ϕ_store_prev.A_prev .== SparseMatrixCOO([8.0 0.0 0.0 -1.0; 0.0 1.0 1.0 0.0]))
+    @test all(ϕ_store_prev.A .== SparseMatrixCOO([2.0 0.0 0.0 -1.0; 0.0 1.0 1.0 0.0]))
+
 
     # test different types
     h = Op(Float32(λ))
@@ -101,6 +133,7 @@ for (op, composite_op, shifted_op) ∈
     y = similar(x)
     ν = 0.1056
     prox!(y, ϕ, x, ν)
+    @test ϕ.full_row_rank == false
     if "$op" == "NormL2"
       y_true = [0.33642, 1.1287, -0.29, 1.14824]
       norm = Op(1.0)