Release v0.3.2: algorithm-verification fixes

CHANGELOG.md

@@ -5,6 +5,41 @@ All notable changes to **DFMethods.jl** will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.3.2] — 2026-05-25
+
+### Added
+
+- **`SpectralThreeTerm`**: two new keyword arguments `alpha_min` and
+  `alpha_max` (defaults `1e-10` and `1e30`) that clamp the spectral
+  coefficient ``ϑ_k^I`` to the interval `[alpha_min, alpha_max]`.
+  The clamp underwrites the strict-descent property
+  ``F(w_k)' d_k ≤ -\alpha_{\min} \|F(w_k)\|^2`` and the trust-region
+  bound ``\|d_k\| ≤ (α_{\max} + 2/\bar α_1)\,\|F(w_k)\|``. In the
+  degenerate case `y_{k-1} = 0` (which forced `ϑ_k^I = 0` previously,
+  giving zero descent), ``ϑ_k^I`` now falls back to `alpha_min`.
+- **`SolodovSvaiterProjection`**: new keyword argument `γ::Float64 =
+  1.0` (relaxation factor; must lie in `(0, 2)`). The projection
+  target becomes `w − γ·λ_k·F(z_k)`, and the Dykstra tolerance scales
+  to `ε_k = (ζ²/2) γ² λ_k² ‖F(z_k)‖²` so Dykstra stops at a constant
+  fractional accuracy of the projection step across all `γ` values.
+  `γ = 1` preserves the prior release behavior byte-for-byte. For
+  `X = RealSpace` the halfspace projection cancels `γ ≤ 1`; the
+  relaxation has effect for `γ > 1` and for non-trivial constraint
+  sets.
+
+### Changed
+
+- **`SpectralThreeTerm.direction!`**: the assembled spectral
+  coefficient is now `clamp(s'y / y'y, alpha_min, alpha_max)` rather
+  than the unconstrained `s'y / y'y`. Default knobs are wide enough
+  (`[1e-10, 1e30]`) that any well-behaved trajectory produces a
+  bit-identical direction to the prior release; only the degenerate
+  fallback (`y_{k-1} = 0`) changes output, replacing the previous
+  zero-descent failure with a strict-descent step.
+- **`SolodovSvaiterProjection`** is now `Base.@kwdef`'d with a single
+  `γ::Float64 = 1.0` field. The zero-argument `SolodovSvaiterProjection()`
+  constructor continues to work and constructs the default rule.
+
 ## [0.3.1] — 2026-05-24
 
 Documentation enhancement release. Adds a new **Tutorial** page to the
@@ -325,6 +360,7 @@ for constrained nonlinear equations $F(x) = 0$ on a closed convex set $X$.
 - `SciMLBase` v2.x
 - `CommonSolve` v0.2.x
 
+[0.3.2]: https://github.com/mmogib/DFMethods.jl/releases/tag/v0.3.2
 [0.3.1]: https://github.com/mmogib/DFMethods.jl/releases/tag/v0.3.1
 [0.3.0]: https://github.com/mmogib/DFMethods.jl/releases/tag/v0.3.0
 [0.2.1]: https://github.com/mmogib/DFMethods.jl/releases/tag/v0.2.1

Project.toml

@@ -1,7 +1,7 @@
 name = "DFMethods"
 uuid = "5fd0b45f-fdf3-4567-a8b1-5e033765ff5d"
 authors = ["Mohammed Alshahrani <mshahrani@kfupm.edu.sa>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"

docs/src/extending.md

@@ -17,7 +17,7 @@ Each subsection points to a runnable file in `examples/` that
 demonstrates a custom subtype end-to-end. Copy one of those files as a
 template for your own extension.
 
-## 0. Contract surface: `ctx`, `cache`, and per-solve state
+## Preliminaries: contract surface (`ctx`, `cache`, per-solve state)
 
 DFMethods's seven extension points use three different conventions for
 passing per-iteration state to user code. This section is the canonical
@@ -252,6 +252,25 @@ implements one method:
 DFMethods.direction!(d::AbstractVector, rule::MyDirection, ctx) -> d
 ```
 
+### Convergence contract
+
+Custom directions must produce $d_k$ satisfying two bounds, which DFMethods
+does **not** check at runtime — they are the rule author's responsibility:
+
+- **Sufficient descent**: $-F(w_k)^\top d_k \ge c\, \|F(w_k)\|^2$ for some
+  constant $c > 0$ independent of $k$.
+- **Bounded growth**: $\|d_k\| \le \bar c\, \|F(w_k)\|$ for some constant
+  $\bar c > 0$ independent of $k$.
+
+The built-in [`SpectralThreeTerm`](@ref) satisfies both: the
+$[\alpha_{\min}, \alpha_{\max}]$ clamp on its spectral coefficient
+$\vartheta_k^I$ gives $-F(w_k)^\top d_k = \vartheta_k^I \|F(w_k)\|^2 \ge
+\alpha_{\min} \|F(w_k)\|^2$, and the $v_k$ denominator together with the
+same clamp yields $\|d_k\| \le (\alpha_{\max} + 2/\bar\alpha_1)\|F(w_k)\|$.
+Violating either bound forfeits the convergence guarantee (see §8).
+
+### `direction!` `ctx`
+
 The context `ctx` is a `NamedTuple` carrying the per-iteration inputs:
 
 | Field                  | Description                                                  |
@@ -290,10 +309,8 @@ function DFMethods.direction!(d, rule::AbubakarNHSCG, ctx)
 end
 ```
 
-Convergence requires the rule to produce a direction with sufficient
-descent ($-F(w)^\top d \geq c\|F(w)\|^2$) and bounded norm
-($\|d\| \leq \bar c\|F(w)\|$). See §8 for the consequences of relaxing
-either property.
+(The convergence contract for `d_k` is stated at the top of this section.
+See §8 for the consequences of violating it.)
 
 ## 3. Iterate update
 
@@ -305,8 +322,8 @@ DFMethods.update_iterate!(x_new::AbstractVector, rule::MyUpdate, ctx) -> x_new
 ```
 
 The `ctx` NamedTuple here has 11 fields (`w`, `d`, `α`, `z`, `Fw`, `Fz`,
-`set`, `k`, `ζ`, `inner_maxiter`, `state`); see §0 for the full table
-and the per-iteration ordering.
+`set`, `k`, `ζ`, `inner_maxiter`, `state`); see *Preliminaries* above for
+the full table and the per-iteration ordering.
 
 Three built-in strategies ship:
 
@@ -324,6 +341,32 @@ $\beta(k) = 1/(k+2)$ for the classical
 [Halpern (1967)](https://doi.org/10.1090/S0002-9904-1967-11864-0)
 iteration.
 
+### Implicit contract for `SolodovSvaiterProjection`
+
+The hyperplane-projection update relies on the multiplier
+
+```math
+\lambda_k \;=\; \frac{F(z_k)^\top (w_k - z_k)}{\|F(z_k)\|^2}
+```
+
+being strictly positive — otherwise the target $w_k - \gamma\lambda_k F(z_k)$
+lies on the wrong side of the separating hyperplane $H_k$ and the
+geometry flips silently. Positivity follows whenever the line search
+enforces an Armijo-style separation
+
+```math
+-F(z_k)^\top d_k \;\ge\; \sigma\, \alpha_k\, \gamma_k\, \|d_k\|^2 \;>\; 0
+```
+
+which together with $z_k = w_k + \alpha_k d_k$ gives
+$F(z_k)^\top (w_k - z_k) = -\alpha_k F(z_k)^\top d_k > 0$. All three
+built-in line searches ([`ConstantBacktrack`](@ref),
+[`ResidualNormBacktrack`](@ref), [`AdaptiveClampedBacktrack`](@ref))
+enforce this. **Custom line searches that don't satisfy the separation
+will produce $\lambda_k \le 0$ silently** — there is no runtime check;
+debug by logging `λ` from inside a custom `update_iterate!` or by
+verifying the separation condition holds at acceptance.
+
 ### Custom iterate-update strategy
 
 A custom rule supplies its own `update_iterate!`; if it needs per-solve

src/iterate_updates.jl

@@ -22,16 +22,44 @@ init_state(::AbstractIterateUpdate, prob, x0, alg) = nothing
 # ============================================================================
 
 """
-    SolodovSvaiterProjection()
+    SolodovSvaiterProjection(; γ = 1.0)
 
 Default iterate-update strategy. Implements the hyperplane projection
 scheme of Solodov & Svaiter (1999): given the trial point `z` with
 residual `F(z)`, construct the separating hyperplane
 `H_k = {x : F(z)' (x − z) ≤ 0}` and project the target
-`w − λ_k F(z)` onto `X ∩ H_k` to tolerance
-`ε_k = (ζ²/2) ‖λ_k F(z)‖²`.
+`w − γ·λ_k F(z)` onto `X ∩ H_k` to tolerance
+`ε_k = (ζ²/2) ‖γ·λ_k F(z)‖² = (ζ²/2) γ² λ_k² ‖F(z)‖²`.
+
+The 1999 scheme assumes exact projection onto `X ∩ H_k`. DFMethods
+realizes it through the inexact-projection refinement standard in the
+broader Solodov–Svaiter inexact-projection lineage: a Dykstra inner
+loop terminated at the tolerance `ε_k` above (driven by `approx_project_X_halfspace!`).
+For `X = RealSpace` the inner solver collapses to a single exact
+halfspace projection regardless of `ε_k`.
+
+The relaxation factor `γ ∈ (0, 2)` controls the step length past the
+separating hyperplane. The convergence bound for this family carries a
+factor of `γ(2 − γ)` (maximized at `γ = 1`); some derivative-free
+projection methods set `γ` larger (e.g. `1.6`–`1.8`) to trade theoretical
+contraction for empirical speed on test problems.
+
+# Parameters
+- `γ::Float64 = 1.0`: relaxation factor; must lie in `(0, 2)`.
+
+# Caveat on under-relaxation (γ < 1)
+For `X = RealSpace` the target `w − γ·λ·F(z)` with `γ < 1` lies on the
+strict-violation side of `H_k`; the exact halfspace projection then
+brings it back to the boundary, recovering the `γ = 1` iterate. In
+short, **`γ ≤ 1` cancels in `RealSpace`**. The relaxation has effect for
+`γ > 1` (target sits inside `H_k`, projection no-op) and for non-trivial
+`X` (joint projection onto `X ∩ H_k` is shaped by both constraints). The
+`ε_k` tolerance is scaled by `γ²` so Dykstra stops to a constant
+fractional accuracy of the projection step regardless of `γ`.
 """
-struct SolodovSvaiterProjection <: AbstractIterateUpdate end
+Base.@kwdef struct SolodovSvaiterProjection <: AbstractIterateUpdate
+    γ::Float64 = 1.0
+end
 
 """
     SolodovSvaiterState
@@ -60,14 +88,15 @@ init_state(::SolodovSvaiterProjection, prob, x, alg) =
     SolodovSvaiterState{eltype(x)}(length(x))
 
 function update_iterate!(x_new::AbstractVector,
-                          ::SolodovSvaiterProjection, ctx)
+                          rule::SolodovSvaiterProjection, ctx)
     w, d, z      = ctx.w, ctx.d, ctx.z
     Fz           = ctx.Fz
     set          = ctx.set
     ζ            = ctx.ζ
     inner_maxiter = ctx.inner_maxiter
     state        = ctx.state
     T            = eltype(w)
+    γ            = T(rule.γ)
 
     # Compute λ_k = F(z)' (w − z) / ‖F(z)‖² and ‖F(z)‖²
     inner_wz   = zero(T)
@@ -80,13 +109,15 @@ function update_iterate!(x_new::AbstractVector,
     end
     λ = inner_wz / Fz_norm_sq
 
-    # target = w - λ·F(z)
+    # target = w - γ·λ·F(z)
     @inbounds @simd for i in eachindex(w)
-        state.proj_target[i] = w[i] - λ * Fz[i]
+        state.proj_target[i] = w[i] - γ * λ * Fz[i]
     end
 
-    # ε_k = (ζ²/2) ‖λ F(z)‖². ζ is Float64 (algorithm parameter); coerce to T.
-    ε_k = T(0.5) * T(ζ)^2 * λ * λ * Fz_norm_sq
+    # ε_k = (ζ²/2) ‖γ·λ F(z)‖² = (ζ²/2) γ² λ² ‖F(z)‖²
+    # — scales with the actual projection-step magnitude so Dykstra
+    # stops to a constant fractional accuracy across γ.
+    ε_k = T(0.5) * γ * γ * T(ζ)^2 * λ * λ * Fz_norm_sq
 
     # Project onto X ∩ H_k
     approx_project_X_halfspace!(x_new, state.proj_target,

src/search_directions.jl

@@ -54,7 +54,7 @@ init_state(::AbstractSearchDirection, prob, x, alg) = nothing
 # ============================================================================
 
 """
-    SpectralThreeTerm(; r=0.1, alpha_bar=1.0)
+    SpectralThreeTerm(; r=0.1, alpha_bar=1.0, alpha_min=1e-10, alpha_max=1e30)
 
 Spectral three-term derivative-free direction. The formula is:
 
@@ -68,23 +68,34 @@ with
 ```
 y_{k-1} = F(w_k) - F(w_{k-1})
 s_{k-1} = (w_k - w_{k-1}) + r · y_{k-1}
-ϑ_k^I   = (s_{k-1}' y_{k-1}) / (y_{k-1}' y_{k-1})
+ϑ_k^I   = clamp((s_{k-1}' y_{k-1}) / (y_{k-1}' y_{k-1}), alpha_min, alpha_max)
 v_k     = max(alpha_bar · ‖d_{k-1}‖ · ‖y_{k-1}‖, ‖F(w_{k-1})‖²)
 β_k     = (F(w_k)' y_{k-1}) / v_k
 ϑ_k^II  = (F(w_k)' d_{k-1}) / v_k
 ```
 
-Satisfies the sufficient-descent and boundedness properties used in
-the convergence proofs for this class of algorithms.
+In the degenerate case `y_{k-1} = 0`, ``ϑ_k^I`` falls back to
+`alpha_min`, preserving the strict-descent property ``F(w_k)' d_k ≤
+-\\alpha_{\\min} \\, \\|F(w_k)\\|^2``.
+
+The uniform bound `alpha_min ≤ ϑ_k^I ≤ alpha_max` is what
+underwrites the sufficient-descent and trust-region properties of
+this class of algorithms.
 
 # Parameters
 - `r`: spectral parameter in the definition of `s_{k-1}` (default 0.1).
 - `alpha_bar`: the parameter ``\\bar{α}_1`` in the definition of `v_k`
   (default 1.0).
+- `alpha_min`: lower clamp on the spectral coefficient ``ϑ_k^I``
+  (default `1e-10`); also the degenerate-case fallback.
+- `alpha_max`: upper clamp on the spectral coefficient ``ϑ_k^I``
+  (default `1e30`).
 """
 Base.@kwdef struct SpectralThreeTerm <: AbstractSearchDirection
     r::Float64         = 0.1
     alpha_bar::Float64 = 1.0
+    alpha_min::Float64 = 1e-10
+    alpha_max::Float64 = 1e30
 end
 
 """
@@ -138,7 +149,10 @@ function direction!(d::AbstractVector, rule::SpectralThreeTerm, ctx)
     v_k = max(T(rule.alpha_bar) * dd_norm * yy_norm, Fwm_sq)
     v_k = max(v_k, eps(typeof(v_k)))   # guard against v_k = 0
 
-    ϑ_I  = yy_sq > 0 ? sy / yy_sq : zero(T)
+    α_min_T = T(rule.alpha_min)
+    α_max_T = T(rule.alpha_max)
+    ϑ_I_raw = yy_sq > 0 ? sy / yy_sq : α_min_T
+    ϑ_I     = clamp(ϑ_I_raw, α_min_T, α_max_T)
     β_k  = Fw_y / v_k
     ϑ_II = Fw_d / v_k