Minres qlp npc (#1026)

* Add linesearch and nonpositive curvature detection to MINRES-QLP

Introduces a 'linesearch' option to MINRES-QLP for inexact Newton methods, enabling detection and handling of nonpositive curvature directions. Updates workspace and statistics to support this feature, adds related documentation, and includes new tests for linesearch behavior and error handling.(not yet working )

* Improve nonpositive curvature detection in MINRES-QLP

Refactored and enabled nonpositive curvature checks in MINRES-QLP with linesearch, correcting the logic for first iteration and updating residual computation.

For Minres,  Need to divide by β, since v has not been normalized yet.

Author: farhadrclass

src/krylov_workspaces.jl

@@ -624,6 +624,7 @@ mutable struct MinresQlpWorkspace{T,FC,S} <: KrylovWorkspace{T,FC,S}
   wₖ         :: S
   M⁻¹vₖ₋₁    :: S
   M⁻¹vₖ      :: S
+  npc_dir    :: S
   x          :: S
   p          :: S
   vₖ         :: S
@@ -642,11 +643,12 @@ function MinresQlpWorkspace(kc::KrylovConstructor)
   wₖ      = similar(kc.vn)
   M⁻¹vₖ₋₁ = similar(kc.vn)
   M⁻¹vₖ   = similar(kc.vn)
+  npc_dir = similar(kc.vn_empty)
   x       = similar(kc.vn)
   p       = similar(kc.vn)
   vₖ      = similar(kc.vn_empty)
   stats = SimpleStats(0, false, false, false, 0, T[], T[], T[], 0.0, "unknown")
-  workspace = MinresQlpWorkspace{T,FC,S}(m, n, Δx, wₖ₋₁, wₖ, M⁻¹vₖ₋₁, M⁻¹vₖ, x, p, vₖ, false, stats)
+  workspace = MinresQlpWorkspace{T,FC,S}(m, n, Δx, wₖ₋₁, wₖ, M⁻¹vₖ₋₁, M⁻¹vₖ, npc_dir, x, p, vₖ, false, stats)
   return workspace
 end
 
@@ -661,9 +663,10 @@ function MinresQlpWorkspace(m::Integer, n::Integer, S::Type)
   x       = S(undef, n)
   p       = S(undef, n)
   vₖ      = S(undef, 0)
+  npc_dir  = S(undef, 0)
   S = isconcretetype(S) ? S : typeof(x)
   stats = SimpleStats(0, false, false, false, 0, T[], T[], T[], 0.0, "unknown")
-  workspace = MinresQlpWorkspace{T,FC,S}(m, n, Δx, wₖ₋₁, wₖ, M⁻¹vₖ₋₁, M⁻¹vₖ, x, p, vₖ, false, stats)
+  workspace = MinresQlpWorkspace{T,FC,S}(m, n, Δx, wₖ₋₁, wₖ, M⁻¹vₖ₋₁, M⁻¹vₖ, npc_dir, x, p, vₖ, false, stats)
   return workspace
 end
 

src/minres.jl

@@ -374,7 +374,8 @@ kwargs_minres = (:M, :ldiv, :linesearch ,:λ, :atol, :rtol, :etol, :conlim, :itm
       if linesearch
         # compute the residual vector and store it in npc_dir
         kscal!(n, sn * sn, npc_dir)  # npc_dir  = sn * sn * npc_dir 
-        kaxpy!(n, -ϕbar * cs, v, npc_dir)  # npc_dir  = npc_dir  - ϕbar * cs * v
+        # Need to divide by β, since v has not been normalized yet.
+        kaxpy!(n, -ϕbar * cs / β, v, npc_dir)  # npc_dir  = npc_dir  - ϕbar * cs * v
       end
 
       # Update x.

src/minres_qlp.jl

@@ -10,6 +10,10 @@
 #
 # S.-C. T. Choi and M. A. Saunders, Algorithm 937: MINRES-QLP for symmetric and Hermitian linear equations and least-squares problems.
 # ACM Transactions on Mathematical Software, 40(2), pp. 1--12, 2014.
+# 
+# Negative curvature detection follows
+# Liu, Yang, and Roosta, MINRES: from negative curvature detection to monotonicity properties,
+# SIAM Journal on Optimization, 32(4), pp. 2636--2661, 2022.
 #
 # Alexis Montoison, <[email protected]>
 # Montreal, September 2019.
@@ -19,7 +23,7 @@ export minres_qlp, minres_qlp!
 """
     (x, stats) = minres_qlp(A, b::AbstractVector{FC};
                             M=I, ldiv::Bool=false, Artol::T=√eps(T),
-                            λ::T=zero(T), atol::T=√eps(T),
+                            linesearch::Bool=false, λ::T=zero(T), atol::T=√eps(T),
                             rtol::T=√eps(T), itmax::Int=0,
                             timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
                             callback=workspace->false, iostream::IO=kstdout)
@@ -60,6 +64,9 @@ For an in-place variant that reuses memory across solves, see [`minres_qlp!`](@r
 * `atol`: absolute stopping tolerance based on the residual norm;
 * `rtol`: relative stopping tolerance based on the residual norm;
 * `Artol`: relative stopping tolerance based on the Aᴴ-residual norm;
+* `linesearch`: if `true`, indicate that the solution is to be used in an inexact Newton method with linesearch. If `true` and nonpositive curvature is detected, the behavior depends on the iteration:
+ – at iteration k = 1, the solver takes the right-hand side (i.e., the preconditioned negative gradient) as the current solution. The same search direction is returned in `workspace.npc_dir`, and `stats.npcCount` is set to 1;
+ – at iteration k > 1, the solver returns the solution from iteration k – 1, the residual from iteration k is a nonpositive curvature direction stored in `stats.npc_dir` and `stats.npcCount` is set to 1;
 * `itmax`: the maximum number of iterations. If `itmax=0`, the default number of iterations is set to `2n`;
 * `timemax`: the time limit in seconds;
 * `verbose`: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every `verbose` iterations;
@@ -77,6 +84,8 @@ For an in-place variant that reuses memory across solves, see [`minres_qlp!`](@r
 * S.-C. T. Choi, *Iterative methods for singular linear equations and least-squares problems*, Ph.D. thesis, ICME, Stanford University, 2006.
 * S.-C. T. Choi, C. C. Paige and M. A. Saunders, [*MINRES-QLP: A Krylov subspace method for indefinite or singular symmetric systems*](https://doi.org/10.1137/100787921), SIAM Journal on Scientific Computing, Vol. 33(4), pp. 1810--1836, 2011.
 * S.-C. T. Choi and M. A. Saunders, [*Algorithm 937: MINRES-QLP for symmetric and Hermitian linear equations and least-squares problems*](https://doi.org/10.1145/2527267), ACM Transactions on Mathematical Software, 40(2), pp. 1--12, 2014.
+* Y. Liu and F. Roosta, [*MINRES: From Negative Curvature Detection to Monotonicity Properties*](https://doi.org/10.1137/21M143666X), SIAM Journal on Optimization, 32(4), pp. 2636--2661, 2022.
+
 """
 function minres_qlp end
 
@@ -100,6 +109,7 @@ def_optargs_minres_qlp = (:(x0::AbstractVector),)
 
 def_kwargs_minres_qlp = (:(; M = I                        ),
                          :(; ldiv::Bool = false           ),
+                         :(; linesearch::Bool = false     ),
                          :(; λ::T = zero(T)               ),
                          :(; atol::T = √eps(T)            ),
                          :(; rtol::T = √eps(T)            ),
@@ -115,7 +125,7 @@ def_kwargs_minres_qlp = extract_parameters.(def_kwargs_minres_qlp)
 
 args_minres_qlp = (:A, :b)
 optargs_minres_qlp = (:x0,)
-kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :verbose, :history, :callback, :iostream)
+kwargs_minres_qlp = (:M, :ldiv, :linesearch, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :verbose, :history, :callback, :iostream)
 
 @eval begin
   function minres_qlp!(workspace :: MinresQlpWorkspace{T,FC,S}, $(def_args_minres_qlp...); $(def_kwargs_minres_qlp...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: AbstractVector{FC}}
@@ -129,6 +139,7 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve
     m == n || error("System must be square")
     length(b) == m || error("Inconsistent problem size")
     (verbose > 0) && @printf(iostream, "MINRES-QLP: system of size %d\n", n)
+    (workspace.warm_start && linesearch) && error("warm_start and linesearch cannot be used together")
 
     # Tests M = Iₙ
     MisI = (M === I)
@@ -139,6 +150,10 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve
 
     # Set up workspace.
     allocate_if(!MisI, workspace, :vₖ, S, workspace.x)  # The length of vₖ is n
+    allocate_if(linesearch, workspace, :npc_dir , S, workspace.x)  # The length of npc_dir is n
+    if linesearch
+      npc_dir = workspace.npc_dir
+    end
     wₖ₋₁, wₖ, M⁻¹vₖ₋₁, M⁻¹vₖ = workspace.wₖ₋₁, workspace.wₖ, workspace.M⁻¹vₖ₋₁, workspace.M⁻¹vₖ
     Δx, x, p, stats = workspace.Δx, workspace.x, workspace.p, workspace.stats
     warm_start = workspace.warm_start
@@ -161,6 +176,7 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve
     # β₁v₁ = Mb
     MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv)
     βₖ = knorm_elliptic(n, vₖ, M⁻¹vₖ)
+    linesearch && kcopy!(n, npc_dir , vₖ)  # npc_dir  ← v; contain the preconditioned initial residual
     if βₖ ≠ 0
       kdiv!(n, M⁻¹vₖ, βₖ)
       MisI || kdiv!(n, vₖ, βₖ)
@@ -217,6 +233,7 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve
     status = "unknown"
     user_requested_exit = false
     overtimed = false
+    stats.indefinite = false
 
     while !(solved || tired || inconsistent || ill_cond_mach || breakdown || user_requested_exit || overtimed)
       # Update iteration index.
@@ -284,6 +301,32 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve
       end
       iter == 1 && (λbarₖ = αₖ)
 
+      # Check for nonpositive curvature
+      if linesearch
+        # if the first iteration, we need to check c₁ * λbar₁ = -α₁
+        if iter == 1
+          cγ = -αₖ
+        else
+          cγ = cₖ₋₁ * λbarₖ
+        end
+        if cγ ≥ 0
+          if iter == 1
+            kcopy!(n, x, b)
+          end
+          # Nonpositive curvature detected.
+          (verbose > 0) && @printf(iostream, "nonpositive curvature detected:  cₖ * λbarₖ = %e\n", cγ)
+          stats.solved = true
+          stats.npcCount = 1
+          stats.niter = iter
+          stats.inconsistent = false
+          stats.timer = start_time |> ktimer
+          stats.status = "nonpositive curvature"
+          workspace.warm_start = false
+          stats.indefinite = true
+          return workspace
+        end
+      end
+
       # Compute and apply current Givens reflection Qₖ.ₖ₊₁
       # [cₖ  sₖ] [λbarₖ] = [λₖ]
       # [sₖ -cₖ] [βₖ₊₁ ]   [0 ]
@@ -297,6 +340,12 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve
       ζₖ      = cₖ * ζbarₖ
       ζbarₖ₊₁ = sₖ * ζbarₖ
 
+      # check for nonpositive curvature
+      if linesearch
+        kscal!(n, sₖ * sₖ, npc_dir)  # npc_dir  = sₖ * sₖ * npc_dir
+        kaxpy!(n, -ζbarₖ₊₁ * cₖ , vₖ₊₁, npc_dir)  # npc_dir  = npc_dir  - ζbarₖ * cₖ * vₖ, this is rₖ
+      end
+      
       # Update the LQ factorization of Rₖ = LₖPₖ.
       # [ λ₁ γ₁ ϵ₁ 0  •  •  0  ]   [ μ₁   0    •    •     •      •      0  ]
       # [ 0  λ₂ γ₂ •  •     •  ]   [ ψ₁   μ₂   •                        •  ]

test/test_minres_qlp.jl

@@ -84,6 +84,45 @@
       @test workspace.stats.status == "user-requested exit"
       @test cb_n2(workspace)
 
+      # Test linesearch
+      A, b = symmetric_indefinite(FC=FC)
+      workspace = MinresQlpWorkspace(A, b)
+      minres_qlp!(workspace, A, b, linesearch=true)
+      x, stats, npc_dir = workspace.x, workspace.stats, workspace.npc_dir
+      @test stats.status == "nonpositive curvature"
+      @test stats.indefinite == true
+      # Verify that the returned direction indeed exhibits nonpositive curvature.
+      # For both real and complex cases, ensure to take the real part.
+      @test real(dot(npc_dir, A * npc_dir)) <= 0
+    
+      # Test Linesearch which would stop on the first call since A is negative definite
+      A, b = symmetric_indefinite(FC=FC; shift = 5)
+      workspace = MinresQlpWorkspace(A, b)
+      minres_qlp!(workspace, A, b, linesearch=true)
+      x, stats, npc_dir = workspace.x, workspace.stats, workspace.npc_dir
+      @test stats.status == "nonpositive curvature"
+      @test stats.niter == 1 
+      @test all(x .== b)
+      @test stats.solved == true
+      @test stats.indefinite == true
+      @test stats.npcCount == 1
+      @test real(dot(npc_dir, A * npc_dir)) <= 0      
+
+      # Test when b^TAb=0 and linesearch is true
+      A, b = system_zero_quad(FC=FC)
+      workspace = MinresQlpWorkspace(A, b)
+      minres_qlp!(workspace, A, b, linesearch=true)
+      x, stats, npc_dir = workspace.x, workspace.stats, workspace.npc_dir
+      @test stats.status == "nonpositive curvature"
+      @test all(x .== b)
+      @test stats.solved == true
+      @test stats.indefinite == true
+      @test real(dot(npc_dir, A * npc_dir)) ≈ 0.0
+
+      # Test if warm_start and linesearch are both true, it should throw an error
+      A, b = symmetric_indefinite(FC=FC)
+      @test_throws MethodError minres_qlp(A, b, warm_start = true, linesearch = true)
+
       @test_throws TypeError minres_qlp(A, b, callback = workspace -> "string", history = true)
     end
   end