solve_shifted_system! method for LBFGS solving step in recursive way

README.md

@@ -78,6 +78,7 @@ Function           | Description
 `size`             | Return the size of a linear operator
 `symmetric`        | Determine whether the operator is symmetric
 `normest`          | Estimate the 2-norm
+`solve_shifted_system!`          | Solves linear system $(B + \sigma I) x = b$, where $B$ is a forward L-BFGS operator and $\sigma \geq 0$.
 
 
 ## Other Operations on Operators

src/lbfgs.jl

@@ -16,6 +16,9 @@ mutable struct LBFGSData{T, I <: Integer}
   b::Vector{Vector{T}}
   insert::I
   Ax::Vector{T}
+  shifted_p::Matrix{T} # Temporary matrix used in the computation solve_shifted_system!
+  shifted_v::Vector{T}
+  shifted_u::Vector{T}
 end
 
 function LBFGSData(
@@ -43,6 +46,9 @@ function LBFGSData(
     inverse ? Vector{T}(undef, 0) : [zeros(T, n) for _ = 1:mem],
     1,
     Vector{T}(undef, n),
+    zeros(T, n, 2*mem),
+    zeros(T, 2*mem),
+    zeros(T, n)
   )
 end
 

src/utilities.jl

@@ -1,4 +1,4 @@
-export check_ctranspose, check_hermitian, check_positive_definite, normest
+export check_ctranspose, check_hermitian, check_positive_definite, normest, solve_shifted_system!
 
 """
   normest(S) estimates the matrix 2-norm of S.
@@ -145,3 +145,73 @@ end
 
 check_positive_definite(M::AbstractMatrix; kwargs...) =
   check_positive_definite(LinearOperator(M); kwargs...)
+
+
+  """
+  solve_shifted_system!(B, 
+                        x, 
+                        σ, 
+                        γ_inv, 
+                        inv_Cx, 
+                        p,
+                        v, 
+                        u)
+
+Solves linear system (B + σI) x = b, where B is a forward L-BFGS operator and σ ≥ 0.
+
+### Parameters
+
+- `B::LBFGSOperator{T,I,F1,F2,F3}`: forward L-BFGS operator.
+- `x::AbstractVector{T}`: The vector representing `-b`.
+- `σ::T`: Nonnegative shift.
+- `inv_Cx::AbstractVector{T}`: A preallocated vector used to store the solution.
+
+### Returns
+
+- `inv_Cx::AbstractVector{T}`: The solution vector `x`, the solution to `(B + σI) x = b`.
+
+### Method
+
+The method uses a two-loop recursion-like approach with modifications to handle the shift `σ`.
+
+### References
+Erway, J. B., Jain, V., & Marcia, R. F. (2014). Shifted L-BFGS Systems. Optimization Methods and Software, 29(5), 992-1004.
+"""
+
+function solve_shifted_system!(
+  B::LBFGSOperator{T, I, F1, F2, F3},
+  x::AbstractVector{T},
+  σ::T,
+  inv_Cx::AbstractVector{T},
+  ) where {T, I, F1, F2, F3}
+  data = B.data
+  insert = data.insert
+
+  γ_inv = 1 / data.scaling_factor
+  inv_c0 = 1 / (γ_inv + σ)
+  @. inv_Cx = inv_c0 * x
+
+  max_i = 2 * data.mem
+  sign_i = 1
+  for i = 1:max_i
+    j = (i + 1) ÷ 2
+    k = mod(insert + j - 1, data.mem) + 1
+    data.shifted_u .= ((i % 2) == 0 ? data.b[k] : data.a[k])
+
+    @. data.shifted_p[:, i] = inv_c0 * data.shifted_u
+
+    sign_t = 1
+    for t = 1:(i - 1)
+      c0 = dot(view(data.shifted_p, :, t), data.shifted_u)
+      c1= sign_t .*data.shifted_v[t]
+      c2 = c1 * c0
+      view(data.shifted_p, :, i) .+= c2 .* view(data.shifted_p, :, t)
+      sign_t = -sign_t
+    end
+
+    data.shifted_v[i] = 1 / (1 + ((i % 2) == 0 ? 1 : -1) * dot(data.shifted_u, view(data.shifted_p, :, i)))
+    inv_Cx .+= sign_i  *data.shifted_v[i] * (view(data.shifted_p, :, i)' * x) .* view(data.shifted_p, :, i)
+    sign_i = -sign_i
+  end
+  return inv_Cx
+end

test/runtests.jl

@@ -15,3 +15,4 @@ include("test_deprecated.jl")
 include("test_normest.jl")
 include("test_diag.jl")
 include("test_chainrules.jl")
+include("test_solve_shifted_system.jl")

test/test_solve_shifted_system.jl

@@ -0,0 +1,43 @@
+using Test
+using LinearOperators
+using LinearAlgebra
+
+function setup_test_val(; M = 5, n = 100, scaling = false)
+  σ = 0.1
+  B = LBFGSOperator(n, mem = 5, scaling = scaling)
+
+  for _ = 1:10
+    s = rand(n)
+    y = rand(n)
+    push!(B, s, y)
+  end
+
+  x = randn(n)
+  z = B * x + σ .* x # so we know the true answer is x
+
+  return B, z, σ, zeros(n), x
+end
+
+function test_solve_shifted_system()
+  @testset "solve_shifted_system! Default setup test" begin
+    # Setup Test Case 1: Default setup from setup_test_val
+    B, z, σ, inv_Cz, x = setup_test_val(n = 100, M = 5)
+
+    result = solve_shifted_system!(B, z, σ, inv_Cz)
+
+    # Test 1: Check if result is a vector of the same size as z
+    @test length(result) == length(z)
+
+    # Test 2: Verify that inv_Cz (result) is modified in place
+    @test result === inv_Cz
+
+    # Test 3: Check if the function produces finite values
+    @test all(isfinite, result)
+
+    # Test 4: Check if inv_Cz is close to the known solution x
+    @test isapprox(inv_Cz, x, atol = 1e-6, rtol = 1e-6)
+  end
+
+end
+
+test_solve_shifted_system()