diff --git a/CHANGELOG.md b/CHANGELOG.md index 885517f..eb7a018 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Project.toml b/Project.toml index 90e1b7e..767038e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DFMethods" uuid = "5fd0b45f-fdf3-4567-a8b1-5e033765ff5d" authors = ["Mohammed Alshahrani "] -version = "0.3.1" +version = "0.3.2" [deps] CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" diff --git a/docs/src/extending.md b/docs/src/extending.md index 6093e7d..c526752 100644 --- a/docs/src/extending.md +++ b/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 diff --git a/src/iterate_updates.jl b/src/iterate_updates.jl index 6a0e010..b449b71 100644 --- a/src/iterate_updates.jl +++ b/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,7 +88,7 @@ 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 @@ -68,6 +96,7 @@ function update_iterate!(x_new::AbstractVector, 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, diff --git a/src/search_directions.jl b/src/search_directions.jl index b5d7d65..0337b69 100644 --- a/src/search_directions.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 73d4a9f..3cda2a7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -197,6 +197,131 @@ end direction!(d, rule, _ctx(Fw, Fw_prev, w, w_prev, d_prev, 1)) @test all(isfinite, d) end + + @testset "constructor: alpha_min / alpha_max defaults + overrides" begin + r_default = SpectralThreeTerm() + @test r_default.alpha_min == 1e-10 + @test r_default.alpha_max == 1e30 + + r_custom = SpectralThreeTerm(; alpha_min = 0.5, alpha_max = 2.0) + @test r_custom.alpha_min == 0.5 + @test r_custom.alpha_max == 2.0 + end + + @testset "ϑ_I clamp: lower bound active" begin + # Force the raw spectral coefficient s'y/y'y to be very small + # (≈ 0) by choosing y_{k-1} ≂̸ 0 with s ⟂ y. Set alpha_min = 0.5 + # so the clamp lifts ϑ_I from ~0 to 0.5. With Fw_prev = 0, β_k + # and ϑ_II vanish, so d_k = -0.5 · F(w_k) exactly. + rng = MersenneTwister(7) + w = randn(rng, n) + y = randn(rng, n) # y_{k-1} + Fw_prev = randn(rng, n) + Fw = Fw_prev .+ y # so F(w_k) - F(w_{k-1}) = y + # Choose w_prev so that w - w_prev = -r·y + e, with e ⟂ y; then + # s = (w - w_prev) + r·y = e ⟂ y, hence s'y = 0. + e = randn(rng, n) + e .-= (dot(e, y) / dot(y, y)) .* y # project e onto y^⊥ + r_val = 0.1 + w_prev = w .+ r_val .* y .- e + d_prev = zeros(n) # so β_k and ϑ_II annihilate + # Disable Fw_prev's effect by zeroing it (v_k stays positive via + # alpha_bar branch with zero d_prev → fallback to eps guard). + # We instead set Fw_prev so that ‖Fw_prev‖² > 0; β_k and ϑ_II + # depend on F_w · y and F_w · d_prev. With d_prev = 0, ϑ_II = 0. + # β_k is generally nonzero unless F_w ⟂ y. Force F_w ⟂ y: + Fw .-= (dot(Fw, y) / dot(y, y)) .* y # now F_w ⟂ y → β_k = 0 + Fw_prev .= Fw .- y # preserve y = F_w - F_w_prev + + rule_lo = SpectralThreeTerm(; r = r_val, alpha_min = 0.5, alpha_max = 1e30) + d = similar(w) + direction!(d, rule_lo, _ctx(Fw, Fw_prev, w, w_prev, d_prev, 1)) + @test d ≈ -0.5 .* Fw atol = 1e-10 + end + + @testset "ϑ_I clamp: upper bound active" begin + # Mirror of the previous test, but force s'y/y'y to be very + # large by making y small and aligning s with y. Then alpha_max + # caps ϑ_I. + rng = MersenneTwister(13) + w = randn(rng, n) + y = 1e-6 .* randn(rng, n) # tiny y → ‖y‖² very small + Fw_prev = randn(rng, n) + # F_w must satisfy F_w - F_w_prev = y AND F_w ⟂ y. + Fw = Fw_prev .+ y + Fw .-= (dot(Fw, y) / dot(y, y)) .* y # F_w ⟂ y → β_k = 0 + Fw_prev .= Fw .- y + r_val = 0.1 + # Pick s = (w - w_prev) + r·y aligned with y so s'y/y'y is huge. + scale = 1e6 + s = scale .* y # s ∥ y, s'y/y'y = scale + w_prev = w .+ r_val .* y .- s + d_prev = zeros(n) + + rule_hi = SpectralThreeTerm(; r = r_val, alpha_min = 1e-10, alpha_max = 2.0) + d = similar(w) + direction!(d, rule_hi, _ctx(Fw, Fw_prev, w, w_prev, d_prev, 1)) + @test d ≈ -2.0 .* Fw atol = 1e-8 + end + + @testset "degenerate y = 0: ϑ_I falls back to alpha_min" begin + # When y_{k-1} = 0 (so yy_sq = 0), ϑ_I should fall back to + # alpha_min, NOT to zero — the pre-fix `zero(T)` fallback would + # have produced d = β·d_prev (no descent contribution from F). + # With Fw_prev = Fw the y vector is zero. v_k stays positive + # via Fwm_sq, so β_k = (F_w · 0)/v_k = 0 and ϑ_II = (F_w·d_prev)/v_k. + # Then d = -alpha_min · F_w + 0 - ϑ_II · 0 = -alpha_min · F_w + # provided d_prev ⟂ F_w; we force that. + n_local = 4 + Fw = [1.0, 0.0, 1.0, 0.0] + Fw_prev = copy(Fw) # → y = 0 + w = randn(MersenneTwister(99), n_local) + w_prev = randn(MersenneTwister(100), n_local) + d_prev = [0.0, 1.0, 0.0, 1.0] # ⟂ Fw + + rule_deg = SpectralThreeTerm(; alpha_min = 0.25, alpha_max = 1e30) + d = similar(Fw) + direction!(d, rule_deg, _ctx(Fw, Fw_prev, w, w_prev, d_prev, 1)) + @test d ≈ -0.25 .* Fw atol = 1e-10 + end + + @testset "default alpha_min/alpha_max do not perturb generic case" begin + # Regression: with non-pathological inputs, sy/yy_sq lands in + # [1e-10, 1e30] and the clamp is a no-op. Output matches the + # pre-fix formula byte-equivalently. + rng = MersenneTwister(2026) + w = randn(rng, n) + w_prev = randn(rng, n) + d_prev = randn(rng, n) + Fw = randn(rng, n) + Fw_prev = randn(rng, n) + + # Recompute ϑ_I, β_k, ϑ_II, d the way the pre-fix code did + # (no clamp) and confirm bitwise-equal output for default knobs. + r_val = 0.1 + ab = 1.0 + y = Fw .- Fw_prev + s = (w .- w_prev) .+ r_val .* y + yy_sq = dot(y, y) + sy = dot(s, y) + Fw_y = dot(Fw, y) + Fw_d = dot(Fw, d_prev) + dd = norm(d_prev) + yyn = sqrt(yy_sq) + Fwm_sq = dot(Fw_prev, Fw_prev) + v_k = max(ab * dd * yyn, Fwm_sq) + v_k = max(v_k, eps(typeof(v_k))) + ϑ_I_ref = sy / yy_sq # raw, no clamp + β_ref = Fw_y / v_k + ϑ_II_ref = Fw_d / v_k + # Sanity: under random inputs the raw ϑ_I lands well inside [1e-10, 1e30] + @test 1e-10 < abs(ϑ_I_ref) < 1e30 + d_ref = -ϑ_I_ref .* Fw .+ β_ref .* d_prev .- ϑ_II_ref .* y + + d = similar(w) + direction!(d, SpectralThreeTerm(), _ctx(Fw, Fw_prev, w, w_prev, d_prev, 1)) + @test d ≈ d_ref atol = 1e-12 + end end # ======================================================================== @@ -321,6 +446,84 @@ end @test norm(sol.u) <= 1e-5 end + @testset "SolodovSvaiterProjection: γ default + override" begin + @test SolodovSvaiterProjection().γ == 1.0 + @test SolodovSvaiterProjection(; γ = 1.8).γ == 1.8 + @test SolodovSvaiterProjection(; γ = 0.5).γ == 0.5 + end + + @testset "SolodovSvaiterProjection: γ = 1 byte-equivalent to v0.3.1 default" begin + # The default γ = 1.0 must reproduce v0.3.1's iterate stream. + # Use a deterministic problem and compare two solves: one with + # SolodovSvaiterProjection() (defaults), one with γ = 1.0 explicit. + f(u, p) = copy(u) + prob = SciMLBase.NonlinearProblem(f, [0.7, -0.3, 0.5]) + alg_default = DFProjection(; iterate_update = SolodovSvaiterProjection(), + linesearch = ConstantBacktrack(), + inertial = NoInertial()) + alg_explicit = DFProjection(; iterate_update = SolodovSvaiterProjection(; γ = 1.0), + linesearch = ConstantBacktrack(), + inertial = NoInertial()) + sol_d = solve(prob, alg_default) + sol_e = solve(prob, alg_explicit) + @test sol_d.u == sol_e.u + @test sol_d.retcode == sol_e.retcode + end + + @testset "SolodovSvaiterProjection: γ > 1 over-relaxation converges" begin + # γ = 1.8 (STTDFPM 2024 experimental default) should still solve + # cleanly on a strongly-monotone test problem. + f(u, p) = copy(u) + prob = SciMLBase.NonlinearProblem(f, [1.0, 1.0]) + alg = DFProjection(; + iterate_update = SolodovSvaiterProjection(; γ = 1.8), + linesearch = ConstantBacktrack(), + inertial = NoInertial(), + abstol = 1e-6, + maxiters = 500, + ) + sol = solve(prob, alg) + @test sol.retcode == SciMLBase.ReturnCode.Success + @test norm(sol.u) <= 1e-5 + end + + @testset "SolodovSvaiterProjection: γ ≤ 1 cancels in RealSpace (documented)" begin + # For X = RealSpace, the halfspace projection brings the target + # w − γ·λ·F(z) back to the H_k boundary for γ ≤ 1, recovering + # the γ = 1 iterate. Confirm by checking that solves with γ = 0.5 + # and γ = 1.0 produce the same iterate sequence. + f(u, p) = copy(u) + prob = SciMLBase.NonlinearProblem(f, [0.8, -0.4]) + alg_half = DFProjection(; iterate_update = SolodovSvaiterProjection(; γ = 0.5), + linesearch = ConstantBacktrack(), + inertial = NoInertial()) + alg_one = DFProjection(; iterate_update = SolodovSvaiterProjection(; γ = 1.0), + linesearch = ConstantBacktrack(), + inertial = NoInertial()) + sol_h = solve(prob, alg_half) + sol_o = solve(prob, alg_one) + @test sol_h.u ≈ sol_o.u atol = 1e-12 + @test sol_h.stats.nsteps == sol_o.stats.nsteps + end + + @testset "SolodovSvaiterProjection: γ > 1 with BoxSet still converges" begin + # Box-constrained smoke: γ > 1 with a non-trivial X. The H_k + # cancellation no longer applies (joint projection onto X ∩ H_k + # is shaped by both constraints), so γ has a genuine effect. + target = [0.3, -0.2] + f(u, p) = u .- target + prob = SciMLBase.NonlinearProblem(f, [1.0, -1.0]; lb = -ones(2), ub = ones(2)) + alg = DFProjection(; + iterate_update = SolodovSvaiterProjection(; γ = 1.6), + inertial = NoInertial(), + abstol = 1e-6, + maxiters = 500, + ) + sol = solve(prob, alg) + @test sol.retcode == SciMLBase.ReturnCode.Success + @test sol.u ≈ target atol = 1e-4 + end + @testset "DirectUpdate: x_{k+1} = project(z, set)" begin f(u, p) = copy(u) prob = SciMLBase.NonlinearProblem(f, [1.0, 1.0])