API Reference

ConstrainedNonlinearProblem(inner::NonlinearProblem, set::AbstractConstraintSet)
ConstrainedNonlinearProblem(f, u0, p = NullParameters();
                            set::AbstractConstraintSet,
                            lb = nothing, ub = nothing, kwargs...)

Wraps a SciMLBase.NonlinearProblem with a feasibility set the algorithm must respect. The constraint becomes part of the problem (where it belongs), not the algorithm — so one DFProjection instance solves many problems with different feasibility sets.

Use this for non-box constraints (HalfSpace, CappedBox, Intersection, UserSet). For pure box constraints, prefer NonlinearProblem(f, u0; lb, ub) directly — SciMLBase's native support.

AbstractDFProjection

Supertype for all derivative-free projection algorithms in DFMethods. Concrete subtypes (e.g., DFProjection in algorithm.jl) wire together a search direction, a line search, an inertial rule, and a constraint set.

Subtypes SciMLBase.AbstractNonlinearAlgorithm, so every concrete algorithm participates in solve(prob::NonlinearProblem, alg; kwargs...) out of the box.

Convergence assumptions (e.g., monotonicity of $F$, closed convex $X$) are documented on the Algorithm page of the manual. The library does not enforce them at runtime — algorithms iterate regardless and either converge or terminate via a stopping criterion.

DFProjection(; direction, linesearch, inertial, iterate_update,
               abstol, reltol, maxiters, maxtime, stopping, ζ,
               inner_maxiter, maxbt, callbacks)

Concrete derivative-free projection algorithm with pluggable components. The constraint set lives on the problem, not on the algorithm — see ConstrainedNonlinearProblem and the lb/ub keywords of NonlinearProblem.

Fields

• direction::AbstractSearchDirection — search-direction rule. Default: SpectralThreeTerm().
• linesearch::LineSearch.AbstractLineSearchAlgorithm — line search. Default: ResidualNormBacktrack().
• inertial::AbstractInertialRule — inertial rule. Default: Inertial(0.25).
• iterate_update::AbstractIterateUpdate — post-line-search iterate update strategy. Default: SolodovSvaiterProjection(). Alternatives: DirectUpdate(), HalpernUpdate(β).
• abstol::Float64 — absolute residual tolerance used to build the default stopping. Default: 1e-6.
• reltol::Float64 — relative residual tolerance (target ‖F(z_k)‖ ≤ reltol·‖F(x_0)‖) used to build the default stopping. 0 disables it. Default: 0.0.
• maxiters::Int — outer-iteration cap used to build the default stopping. Default: 2000.
• maxtime::Float64 — wall-clock budget in seconds used to build the default stopping. Inf disables it. Default: Inf.
• stopping::AbstractStoppingCriterion — full stopping rule. If not supplied, built from the knobs above as AnyOf(AbsResidualTol(abstol)[, RelResidualTol(reltol)], MaxIters(maxiters)[, MaxTime(maxtime)]) — the bracketed criteria appear only when reltol > 0 / maxtime is finite.
• ζ::Float64 — approximate-projection tolerance factor used by SolodovSvaiterProjection. Default: 0.5.
• inner_maxiter::Int — max inner-projection iterations (Dykstra). Default: 500.
• maxbt::Int — max line-search backtracks per iteration. Default: 50.
• callbacks::Vector{<:AbstractCallback} — observer callbacks fired during the solve. Default: empty.

Stopping criteria

abstol, reltol, maxiters, and maxtime are convenience knobs that build the default stopping rule (mirroring the SciML common-solver options of the same name). For composite or domain-specific criteria, pass stopping = AnyOf(StepNormTol(...), DirectionNormTol(...), …); the knob fields are still stored (for introspection and SciMLBase kwarg overrides) but step! ignores them in favour of the supplied stopping. When a stopping rule is supplied explicitly, solve-time tolerance keywords (abstol/reltol/maxiters/maxtime) cannot be re-applied and are reported via a verbose warning.

DFProjectionCache{T, Alg, F, S, DirState, LSCache, IUpState, StopState}

Mutable per-solve state. Created via init_cache(F, x0, alg). All per-iteration buffers are pre-allocated; step! does not allocate.

Parametric on element type T <: AbstractFloat, derived from eltype(x0) (or Float64 if eltype(x0) is non-floating, e.g. Int).

Common state (vectors)

• x, x_prev, w, w_prev, z, d, d_prev — iteration vectors, Vector{T}.
• Fw, Fw_prev, Fz — function values, Vector{T}.
• x_new, scratch1 — generic scratch, Vector{T}.

Pluggable component state (typed via parameters)

• direction_state::DirState — init_state(alg.direction, ...). Default nothing.
• line_search_cache::LSCache — populated by the line-search component. Default nothing (Stage 2 transition state; filled by Stage 3 once LineSearch.jl is adopted).
• iterate_update_state::IUpState — init_state(...) for the iterate-update strategy. In Stage 2 this is always a SolodovSvaiterState (holding the projection scratch and Dykstra buffers).
• stopping_state::StopState — init_state(alg.stopping, ...). Default nothing.

Common scalars

• k::Int — iteration counter (0-based).
• n_evals::Int — total F evaluations.
• converged::Bool, done::Bool, retcode::Symbol, resid::T.
• F0_norm::T — ‖F(x_0)‖, used by RelResidualTol.
• t_start::Float64 — wall-clock seconds at solve start (for MaxTime).
• α_prev::T — previous line-search step size, surfaced to direction! via ctx.α_prev.

init_cache(F, x0, alg) -> DFProjectionCache

Allocate per-solve buffers and prepare the initial state. Projects x0 onto set to ensure the feasibility invariant x_0 ∈ X. Evaluates F once at the projected x_0 to populate F0_norm (used by RelResidualTol) and starts the wall-clock for MaxTime.

DFSciMLCache

Wrapper cache returned by CommonSolve.init(prob, alg; …). Carries the original NonlinearProblem, the user-facing algorithm, the inner Phase 2 DFProjectionCache, and the resolved verbose flag (used by solve! to decide whether to warn on a non-Success exit). CommonSolve.solve! drives the inner cache to termination and packages the result as a NonlinearSolution.

CommonSolve.init(prob::NonlinearProblem, alg::DFProjection;
                 abstol, reltol, maxiters, maxtime, verbose, kwargs...)
    -> DFSciMLCache

Build a cache for solve(prob, alg; …). Honors the SciML common-solver options that map onto DFMethods' callback-based stopping system:

• abstol → AbsResidualTol(abstol)
• reltol → RelResidualTol(reltol) (target ‖F(z_k)‖ ≤ reltol·‖F(x_0)‖)
• maxiters → MaxIters(maxiters)
• maxtime → MaxTime(maxtime) seconds (nothing ⇒ no limit)
• verbose → toggles the non-convergence warning emitted by solve!

Each defaults to the matching field on alg; any that is supplied overrides it by rebuilding the default stopping rule. If alg was built with an explicit stopping=, these tolerance keywords cannot be applied and (when verbose) a warning is emitted.

All other SciML common-solver keywords — termination_condition, internalnorm, alias_u0, show_trace, store_trace, trace_level — are accepted and absorbed without error. DFMethods uses its own callback-based termination and the Euclidean residual norm, so solve(prob, ::DFProjection; any_standard_kwarg…) never throws.

CommonSolve.step!(cache::DFSciMLCache) -> cache

Advance the inner cache by one outer iteration of the algorithm. No-op if the inner cache has already terminated. Returns the wrapper cache so calls chain.

Usage:

using NonlinearSolve, DFMethods
cache = init(prob, alg)
while !cache.inner.done
    step!(cache)
end
sol = solve!(cache)

For batched runs, solve(prob, alg) (or equivalently solve!(init(prob, alg))) is the higher-level entry point that handles the loop internally.

CommonSolve.solve!(cache::DFSciMLCache) -> NonlinearSolution

Drive the inner cache to termination, then build and return a SciMLBase.NonlinearSolution.

AbstractSearchDirection

Supertype for search-direction rules. Concrete subtypes implement the single 3-argument method

direction!(d, rule, ctx)

where the output direction d is written in place and ctx is a NamedTuple carrying the per-iteration cache state. The vector fields alias into the DFProjectionCache buffers, parametric on element type T, where T = eltype(x0) <: AbstractFloat (with Float64 fallback if eltype(x0) is not floating). The fields are:

A rule that doesn't need a given field simply ignores it. New cache fields can be added in future versions without breaking existing direction methods — they will just be unused by older code.

For ctx.k == 0, only ctx.Fw is meaningful; the other fields may be uninitialized. Rules typically set d = -Fw at k = 0.

The convergence theory for this class of algorithms requires sufficient descent (-Fw' d ≥ c‖Fw‖²) and boundedness (‖d‖ ≤ c̄‖Fw‖). It is the rule author's responsibility to ensure these hold.

SpectralThreeTerm(; r=0.1, alpha_bar=1.0, alpha_min=1e-10, alpha_max=1e30)

Spectral three-term derivative-free direction. The formula is:

d_0 = -F(w_0)
d_k = -ϑ_k^I · F(w_k) + β_k · d_{k-1} - ϑ_k^II · y_{k-1}        (k ≥ 1)

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   = 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

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).

direction!(d, rule, ctx)

In-place compute $d = d_k$ per rule. See AbstractSearchDirection for the ctx NamedTuple shape.

ConstantBacktrack(; σ=0.01, ρ=0.6, maxbt=50)

Backtracking line search with γ_k ≡ 1 — Armijo's classical form.

ResidualNormBacktrack(; σ=0.01, ρ=0.6, maxbt=50)

Backtracking with γ_k = ‖F(z_k)‖. Default line search for DFProjection().

AdaptiveClampedBacktrack(; σ=0.01, ρ=0.6, lo=1e-4, Δ_init=10.0, maxbt=50)

Backtracking with γ_k = clamp(‖F(z_k)‖, lo, lo + Δ). The Δ upper-bound range is held on the cache (mutable; could adapt across solves in a future extension; not currently updated).

AbstractIterateUpdate

Supertype for the iterate-update strategy: given the inputs (w, d, α, z, Fw, Fz, set, k, ζ, inner_maxiter, state), produce x_{k+1}. Concrete subtypes implement

update_iterate!(x_new, rule::AbstractIterateUpdate, ctx)

writing into x_new and returning it. ctx is a NamedTuple; see SolodovSvaiterProjection, DirectUpdate, and HalpernUpdate for built-in strategies.

The strategy may declare its per-solve state via init_state(rule, prob, x0, alg). The state is held on DFProjectionCache.iterate_update_state and surfaced to update_iterate! via ctx.state.

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)‖² = (ζ²/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 γ.

DirectUpdate()

Stateless update: x_{k+1} = project(z, set). Simplest possible scheme — take the trial point and project it onto the feasible set. Useful as a baseline for comparison and as a building block in higher-level methods.

HalpernUpdate(β)

Halpern-style averaging with feasibility maintenance:

x_{k+1} = P_X(β · x_0 + (1 − β) · z_k)

where P_X is projection onto the problem's feasibility set X. The trial point z_k = w_k + α_k d_k is generally infeasible, so the final projection is required to maintain the framework invariant x_k ∈ X.

The parameter β is either a scalar (constant across iterations) or a callable β(k) -> Float64 (per-iteration schedule, e.g. k -> 1 / (k + 2) for the classical schedule with strong convergence to the nearest fixed point).

Reference

Halpern, B. (1967). Fixed Points of Nonexpanding Maps. Bulletin of the AMS 73(6): 957–961. doi:10.1090/S0002-9904-1967-11864-0.

update_iterate!(x_new, rule::AbstractIterateUpdate, ctx) -> x_new

In-place: write the next iterate into x_new per rule's logic. ctx is a NamedTuple with fields w, d, α, z, Fw, Fz, set, k, ζ, inner_maxiter, state. See AbstractIterateUpdate.

init_state(component, prob, x0, alg) -> state

Allocate per-solve state for a pluggable algorithm component (search direction, line search, iterate-update strategy, stopping criterion). The returned state is held on DFProjectionCache in the slot matching the component's role. Components that are stateless (e.g. SpectralThreeTerm) return nothing; components that need per-iteration scratch (e.g. the default Solodov–Svaiter iterate-update strategy with its Dykstra buffers) return a small struct holding their allocations.

Default: nothing. Concrete components override.

This is an internal helper (not exported). The DFProjectionCache constructor calls it once per pluggable component during init_cache.

AbstractInertialRule

Supertype for inertial-extrapolation rules. A rule produces the inertial coefficient $θ_k$ from the iteration index, current iterate, and previous iterate. The inertial point is then

w_k = x_k + θ_k (x_k - x_{k-1}).

Inertial(θ = 0.25)

Standard summable-step inertial-extrapolation rule:

θ_k = min(θ, 1 / (k² · ‖x_k - x_{k-1}‖))    if x_k ≠ x_{k-1}
      θ                                       otherwise

θ ∈ (0, 1); 0.25 is a common default value in the literature.

The 1/k² cap ensures $\sum_k θ_k \|x_k - x_{k-1}\| \le \sum_k 1/k^2 < ∞$, which is the summability condition that drives the convergence proof for inertial schemes of this class. Lineage: Alvarez & Attouch (2001) introduced the heavy-ball / inertial idea for monotone operators; Maingé (2008) gave the modern form for inertial KM-type algorithms; the specific 1/k² cap used here is the form adopted by recent derivative-free projection methods.

julia> Inertial().θ
0.25

julia> Inertial(0.4).θ
0.4

julia> inertial_coef(Inertial(0.25), 0, [1.0], [0.0])    # k=0: no inertia
0.0

NoInertial()

Disables inertia: $w_k = x_k$. Useful for ablation comparisons against the inertial variant.

julia> inertial_coef(NoInertial(), 5, [1.0, 2.0], [0.0, 1.0])
0.0

inertial_coef(rule, k, xk, xkm1) -> θ_k::Float64

Compute the inertial coefficient. k is the iteration index (0-based); xk is the current iterate; xkm1 is the previous iterate.

apply_inertial!(w, rule, k, xk, xkm1) -> (w, θ_k)

Compute $w = x_k + θ_k (x_k - x_{k-1})$ in place, returning $θ_k$ too.

AbstractCallback

Supertype for solve-time callbacks. Implement on_event!(cb, cache, event::Symbol) to observe state at specific lifecycle events. Default behavior: observe-only (return nothing).

Stopping criteria are a special case — see AbstractStoppingCriterion, which subtypes AbstractCallback and returns (stopped::Bool, retcode::Symbol) from on_event! to signal termination.

Events fired in v0.2:

The minimal 4-event set is intentional. More events (:pre_iter, :post_direction, :post_iterate_update) may be added in a future release if real need emerges.

on_event!(cb, cache, event::Symbol)

Called by the algorithm at each instrumented lifecycle event. For AbstractCallback subtypes that observe only (e.g. HistoryCallback, LoggingCallback), return value is ignored. For AbstractStoppingCriterion subtypes, return (stopped::Bool, retcode::Symbol).

Default (observer): nothing.

HistoryCallback(; fields = HISTORY_FIELDS)
HistoryCallback(extractor::Function)

Per-iteration history accumulator. Two constructor paths:

• Pre-declared field names (the common case): pick any subset of HISTORY_FIELDS. Each iter, on_event!(cb, cache, :post_iter) appends a NamedTuple with those fields to cb.history.
• Custom extractor: HistoryCallback(extractor = cache -> NamedTuple) bypasses fields and uses the user's function to produce each row. Useful for derived quantities (e.g. max(abs, F), condition numbers).

Constructor errors if both fields and extractor are supplied.

Example

hist = HistoryCallback(fields = (:k, :F_norm))
sol = solve(prob, DFProjection(callbacks = [hist]))
hist.history  # Vector{NamedTuple} with .k and .F_norm per row

LoggingCallback(; io = stdout, columns = HISTORY_FIELDS, every = 1,
                  header_every = 50, footer = true)

Per-iteration table logger. Prints:

• A header row at :initialize
• One row per iteration at :post_iter, respecting every (1 = every iter, 10 = every 10th iter)
• Header re-prints every header_every rows (set to 0 to disable)
• A footer summary at :terminate (skip via footer = false)

HISTORY_FIELDS

Available scalar fields for HistoryCallback and the default column set for LoggingCallback.

julia> HISTORY_FIELDS
(:k, :F_norm, :d_norm, :α, :n_evals, :resid, :elapsed)

AbstractStoppingCriterion <: AbstractCallback

Supertype for stopping criteria. Concrete subtypes implement on_event!(crit, cache, event::Symbol) -> (Bool, Symbol) for the events they care about.

The default fallback returns (false, :Default) — i.e., never stops. Concrete criteria override on_event! for their relevant event(s).

AbsResidualTol(abstol)

Stop when $\|F(z_k)\| \le \text{abstol}$ (post line search). Retcode :Success.

RelResidualTol(rtol; abstol=0)

Stop when $\|F(z_k)\| \le \text{abstol} + \text{rtol} \cdot \|F(x_0)\|$ (post line search). Retcode :Success.

StepNormTol(xtol)

Stop when $\|x_k - x_{k-1}\| \le \text{xtol}$ at end of iteration. Catches stalls. Retcode :Stalled.

DirectionNormTol(dtol)

Stop when $\|d_k\| \le \text{dtol}$. Retcode :Stalled. Pair with AbsResidualTol (via AnyOf) for a complete picture — a vanishing direction is not the same as a vanishing residual.

MaxIters(maxiters)

Stop when cache.k >= maxiters. Retcode :MaxIters.

MaxTime(maxtime)

Stop when elapsed wall-clock time exceeds maxtime seconds. Retcode :MaxTime.

MaxFEvals(maxevals)

Stop when total function evaluations reach maxevals. Retcode :MaxFEvals.

UserStop(f)

User-supplied predicate. f(cache) -> (stopped::Bool, retcode::Symbol). Called at :post_iter by default.

AnyOf(criteria...)

Composite criterion: stops when any constituent stops. Delegates each event to all sub-criteria in declaration order; the first to fire wins and its retcode propagates.

AbstractConstraintSet

Supertype for all constraint sets. Concrete subtypes implement project!(y, x, set) that writes the orthogonal projection of x onto the set into y and returns y.

RealSpace()

The whole of $\mathbb{R}^n$ — projection is the identity. Use to disable constraints (e.g., for unconstrained comparison with NonlinearSolve.jl's SimpleDFSane).

julia> project([1.0, -2.0, 3.0], RealSpace())
3-element Vector{Float64}:
  1.0
 -2.0
  3.0

BoxSet(lower, upper)

Box constraint $\{x : \text{lower}_i \le x_i \le \text{upper}_i\}$. lower and upper are vectors of equal length with lower[i] ≤ upper[i] element-wise. Use BoxSet(fill(a, n), fill(b, n)) for uniform scalar bounds.

julia> project([2.0, -3.0, 0.5], BoxSet([-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]))
3-element Vector{Float64}:
  1.0
 -1.0
  0.5

HalfSpace(a, c)

Half-space $\{x : a^\top x \le c\}$. a is the normal vector (must be nonzero), c the offset. Projection:

y = x                          if a' x ≤ c
y = x − ((a' x − c)/‖a‖²) a   otherwise

The hyperplane $H_k$ in the Solodov–Svaiter step is a HalfSpace.

Intersection(set1, set2; maxiter=200, tol=1e-12)

Intersection of two closed convex sets. Projection uses Dykstra's algorithm (alternating projection with two correction sequences). For $\text{Box} \cap \text{HalfSpace}$ (the polyhedral set $\Omega$ realized by CappedBox), Dykstra converges geometrically.

For workloads that need general polyhedral projection, a direct QP solver (JuMP + HiGHS) will be added in Phase 2 as an alternative backend.

CappedBox(a, b, c)

The polyhedral set $\Omega(a, b, c) = \{x \in \mathbb{R}^n : a \le x_i \le b,\ \sum_i x_i \le c\}$ — a uniform-bound box intersected with a single budget-style halfspace. a, b, c are scalars; the dimension n is inferred from the input vector at projection time. This set is a common experimental domain in the convex-constrained nonlinear-equations literature.

Projection

Closed-form via bisection on a 1D Lagrange multiplier:

1. Compute the box projection $y_B = \mathrm{clamp}(x, a, b)$.
2. If $\mathbf{1}^\top y_B \le c$ → return $y_B$ (constraint inactive).
3. Otherwise find $λ^* > 0$ such that $\sum_i \mathrm{clamp}(x_i - λ^*, a, b) = c$, then return $\mathrm{clamp}(x - λ^* \mathbf{1}, a, b)$.

The sum is piecewise-linear monotone-decreasing in λ, so bisection converges in O(log(1/tol)) steps.

Feasibility

The set is empty iff n·a > c. project! throws ArgumentError in that case (the feasibility check requires n, so it cannot happen at construction time).

UserSet(proj!)

Constraint set defined by a user-supplied in-place projection proj!(y, x) that writes the projection of x onto the set into y. The user is responsible for ensuring proj! projects onto a closed convex set — otherwise the algorithm's convergence guarantees do not hold.

project!(y::AbstractVector, x::AbstractVector, set::AbstractConstraintSet) -> y

In-place orthogonal projection of x onto set. Writes the projection into y and returns y. Each concrete AbstractConstraintSet subtype provides its own method; see the individual set docstrings for the projection formula or algorithm used.

For an allocating wrapper, see project.

project(x, set) -> Vector

Allocating out-of-place projection. For tests and small problems. Inside the algorithm hot loop, always use project!.

approx_project_X_halfspace!(out, target, X, a, c, ε,
                            p_buf, q_buf, scratch, out_prev;
                            maxiter=500) -> out

Approximate projection of target onto $X \cap H$ with $H = \{x : a^\top x \le c\}$, via Dykstra. Inner buffers (p_buf, q_buf, scratch, out_prev) are pre-allocated by the caller (typically the algorithm cache). Stops when $\|y - y_{\text{prev}}\|^2 \le \epsilon$ or maxiter reached.