From 1b5a746c0f0bd5d99c861e990c734d12dbba5256 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Thu, 12 Jun 2025 13:17:17 +0200
Subject: [PATCH 1/5] minor fixes in shiftedCompositeNormL2

---
 src/compositeNormL2.jl        |  6 ++++--
 src/shiftedCompositeNormL2.jl | 16 ++++++++++++----
 2 files changed, 16 insertions(+), 6 deletions(-)

diff --git a/src/compositeNormL2.jl b/src/compositeNormL2.jl
index 8056fd4..2e79130 100644
--- a/src/compositeNormL2.jl
+++ b/src/compositeNormL2.jl
@@ -35,19 +35,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
 
diff --git a/src/shiftedCompositeNormL2.jl b/src/shiftedCompositeNormL2.jl
index 368cb1f..cb48ab1 100644
--- a/src/shiftedCompositeNormL2.jl
+++ b/src/shiftedCompositeNormL2.jl
@@ -33,6 +33,7 @@ mutable struct ShiftedCompositeNormL2{
   c!::F0
   J!::F1
   A::M
+  A_prev::Union{Nothing, M}  # (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 +42,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,6 +62,8 @@ 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)
@@ -68,6 +73,7 @@ mutable struct ShiftedCompositeNormL2{
       c!,
       J!,
       A,
+      A_prev,
       shifted_spmat,
       spfct,
       b,
@@ -76,6 +82,7 @@ mutable struct ShiftedCompositeNormL2{
       dq,
       p,
       dp,
+      false
     )
   end
 end
@@ -94,7 +101,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"
@@ -110,6 +117,7 @@ function prox!(
   atol = eps(T)^0.3,
   max_time = 180.0,
 ) where {T <: Real, F0 <: Function, F1 <: Function, M <: AbstractMatrix{T}, V <: AbstractVector{T}}
+  @assert ν > zero(T)
   start_time = time()
   θ = T(0.8)
   α = zero(T)
@@ -134,8 +142,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!(

From 2a452ee3c2fefac5f1733f66e293fc4a218c9eb1 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sun, 30 Nov 2025 17:27:25 -0500
Subject: [PATCH 2/5] add unit tests for store previous Jacobian option

---
 test/runtests.jl | 29 ++++++++++++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 0c8470f..77a7590 100644
--- a/test/runtests.jl
+++ b/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,31 @@ 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 different types
     h = Op(Float32(λ))
@@ -101,6 +127,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)

From 8cd80132d52da43b2906664cb34ffd7b52fab2f3 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Sun, 30 Nov 2025 17:33:46 -0500
Subject: [PATCH 3/5] add store_previous_jacobian doc

---
 src/compositeNormL2.jl        | 4 +++-
 src/shiftedCompositeNormL2.jl | 5 ++++-
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/compositeNormL2.jl b/src/compositeNormL2.jl
index 2e79130..17dc675 100644
--- a/src/compositeNormL2.jl
+++ b/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,
diff --git a/src/shiftedCompositeNormL2.jl b/src/shiftedCompositeNormL2.jl
index cb48ab1..94c8d64 100644
--- a/src/shiftedCompositeNormL2.jl
+++ b/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,6 +21,9 @@ 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,

From 6a03b5df608c7ea77de51ac1e33fc5582b35da3c Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Mon, 1 Dec 2025 12:25:37 -0500
Subject: [PATCH 4/5] update previous Jacobian on shift!

---
 src/ShiftedProximalOperators.jl | 1 +
 test/runtests.jl                | 6 ++++++
 2 files changed, 7 insertions(+)

diff --git a/src/ShiftedProximalOperators.jl b/src/ShiftedProximalOperators.jl
index de750c1..b21055a 100644
--- a/src/ShiftedProximalOperators.jl
+++ b/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 ψ
diff --git a/test/runtests.jl b/test/runtests.jl
index 54ee48b..ed17918 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -89,6 +89,12 @@ for (op, composite_op, shifted_op) ∈
     @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(λ))

From 362053c3a8c7348c453173c4cae0af6ec6e92944 Mon Sep 17 00:00:00 2001
From: Maxence Gollier <mgollier45@gmail.com>
Date: Tue, 9 Dec 2025 09:27:32 -0500
Subject: [PATCH 5/5] remove Union type for A_prev and make it parametric

---
 src/shiftedCompositeNormL2.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/shiftedCompositeNormL2.jl b/src/shiftedCompositeNormL2.jl
index 94c8d64..5661f91 100644
--- a/src/shiftedCompositeNormL2.jl
+++ b/src/shiftedCompositeNormL2.jl
@@ -30,13 +30,14 @@ mutable struct ShiftedCompositeNormL2{
   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::Union{Nothing, M}  # (Optional) can be used to store the previous Jacobian, useful for quasi-Newton approximations
+  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
@@ -71,7 +72,7 @@ mutable struct ShiftedCompositeNormL2{
     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!,
@@ -113,13 +114,13 @@ 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)