Nonlinear Preconditioning and Iterate Limiting (PCNR)

Newton-type methods can fail or crawl on systems whose residuals are violently nonlinear — the classic example being the exponential I-V curves of semiconductor devices. Two classical remedies are nonlinear preconditioning (transform the residual so Newton sees a tamer function) and iterate limiting (clip each Newton update to a physically trusted move). NonlinearSolve.jl exposes both as solver options:

\[G\big(f(H(\tilde{u}, u_k), p), \tilde{u}, p\big) = 0\]

• precondition — a left preconditioner G(fu, u, p) applied to the residual. The solver replaces the residual with the root-equivalent composition u -> G(f(u, p), u, p)everywhere: function evaluations, automatic-differentiation Jacobians, line-search merit functions, and termination criteria. G must be root-preserving: G(r, u, p) = 0 if and only if r = 0.
• postcondition — a right preconditioner / corrector H(u_proposed, u_prev, p, cache) applied to every iterate a solver is about to accept, before the residual is evaluated or convergence is tested there (the initial guess is corrected once as H(u0, u0, p, nothing)). H must leave solutions fixed: H(u, u, p, cache) = u at any root.

Both are ordinary solve keywords, like abstol or termination_condition. They can be passed at solve/init time, or carried on the problem and forwarded like any other option — a solve-time value wins, matching how alias behaves:

solve(prob, NewtonRaphson(); postcondition = H)          # late-bound, problem untouched
NonlinearProblem(f, u0, p; postcondition = H)            # carried on the problem

This is the solver-composition viewpoint of Brune, Knepley, Smith & Tu, Composing Scalable Nonlinear Algebraic Solvers (SIAM Review 57(4), 2015): G and H may depend on the current iterate, and the solver freezes that dependence within each iteration. postcondition corresponds to PETSc SNES's SNESLineSearchSetPostCheck hook, and it is exactly the corrector phase of the Predictor/Corrector Newton-Raphson (PCNR) method of Aadithya, Keiter & Mei (Sandia) used to replace ad-hoc limiting in circuit simulators.

Left preconditioning: taming an exponential residual

Consider solving the diode equation for the junction voltage v at a target current:

\[I_s\big(e^{v/V_t} - 1\big) - I_{\text{target}} = 0\]

From an initial guess of v = 2 volts, the residual is astronomically large ($\sim 10^{20}$) and its slope is even larger, so every Newton step retreats by only $\mathcal{O}(V_t) = 25\,\text{mV}$ and the solve creeps:

using NonlinearSolve

p = (; Is = 1.0e-14, Vt = 0.025, It = 1.0e-2)
f_diode(v, p) = p.Is * expm1(v / p.Vt) - p.It

prob = NonlinearProblem(f_diode, 2.0, p)
sol_plain = solve(prob, NewtonRaphson())
sol_plain.retcode, sol_plain.u, sol_plain.stats.nsteps

(SciMLBase.ReturnCode.Success, 0.6907755278982387, 57)

The fix is a residual compression. asinh behaves like log for large arguments and like the identity near zero, is odd and strictly monotone — so asinh(F(v)) has exactly the same root, but is nearly affine in v wherever the exponential dominates:

G(fu, u, p) = asinh(fu)

sol_G = solve(prob, NewtonRaphson(); precondition = G)
sol_G.retcode, sol_G.u, sol_G.stats.nsteps

(SciMLBase.ReturnCode.Success, 0.6907755278983617, 9)

Newton now converges in a handful of steps instead of dozens, and prob itself is untouched — the option applies to this solve only. Since the composition is what the solver differentiates (via AD), the Newton step, the line-search merit function, and the convergence test are all consistently formulated in the preconditioned residual — note that sol_G.resid is the preconditioned residual G(f(u)).

If you provide jac, jvp, vjp, jac_prototype, or sparsity on the NonlinearFunction together with precondition, they must describe the derivatives and structure of the composed map, not of the raw f.

For NonlinearLeastSquaresProblems, precondition re-weights the least-squares objective to ‖G(f(u))‖²: for consistent (zero-residual) systems the solution is unchanged, but for genuinely overdetermined fits it changes the minimizer — which makes the option a natural place for residual weighting.

Iterate limiting for a circuit model: the PCNR method

Now the classic circuit-simulation problem, following Aadithya, Keiter & Mei, Predictor/Corrector Newton-Raphson (PCNR): A Simple, Flexible, Scalable, Modular, and Consistent Replacement for Limiting in Circuit Simulation, Scientific Computing in Electrical Engineering (2020): a voltage source Vs in series with a resistor R feeding a diode to ground. SPICE-family simulators solve the nodal equations with Newton-Raphson plus junction-voltage limiting (pnjlim): a proposed update to a diode voltage is clipped to a logarithmic move relative to the previous iterate, since a volt-sized overshoot puts exp(v/Vt) far outside the region where the linearization means anything.

The PCNR formulation makes the limited quantity an explicit unknown. Our unknowns are u = [v, vj] — the node voltage and the junction voltage — tied together by a consistency equation:

cp = (; Vs = 5.0, R = 1.0e3, Is = 1.0e-14, Vt = 0.025)

function circuit!(r, u, p)
    v, vj = u[1], u[2]
    r[1] = (v - p.Vs) / p.R + p.Is * expm1(vj / p.Vt)  # KCL at the node
    r[2] = vj - v                                       # junction consistency
    return nothing
end

prob_c = NonlinearProblem(circuit!, zeros(2), cp)
sol_c_plain = solve(prob_c, NewtonRaphson(); maxiters = 1000)
sol_c_plain.retcode, sol_c_plain.stats.nsteps

(SciMLBase.ReturnCode.Success, 179)

Plain Newton eventually gets there, but it takes a couple hundred iterations of millivolt-sized creep after the first step overshoots the junction voltage to several volts. Now add the classic SPICE3 pnjlim limiter as a postcondition. It sees the proposed iterate and the previous accepted iterate — precisely the two pieces of information limiting needs:

function pnjlim(vnew, vold, vt, vcrit)
    if vnew > vcrit && abs(vnew - vold) > 2vt
        if vold > 0
            arg = 1 + (vnew - vold) / vt
            vnew = arg > 0 ? vold + vt * log(arg) : vcrit
        else
            vnew = vt * log(vnew / vt)
        end
    end
    return vnew
end

vcrit = cp.Vt * log(cp.Vt / (sqrt(2) * cp.Is))

# corrector: limit the junction voltage update, leave the node voltage alone
H!(up, uprev, p, cache) = (up[2] = pnjlim(up[2], uprev[2], p.Vt, vcrit); nothing)

sol_c_lim = solve(prob_c, NewtonRaphson(); postcondition = H!, maxiters = 1000)
sol_c_lim.retcode, sol_c_lim.stats.nsteps

(SciMLBase.ReturnCode.Success, 11)

An order of magnitude fewer iterations, on the same problem object. This is the PCNR method:

1. Predictor — the solver's ordinary Newton step on the augmented system.
2. Corrector — the postcondition applies the limiting functions to the proposed iterate.
3. Consistency — the framework re-evaluates the residual and Jacobian at the corrected iterate, so the next linearization matches the state the devices were actually evaluated at. (Traditional device-level limiting breaks exactly this property, which is the inconsistency PCNR was designed to remove.)

Since the residual is evaluated at the limited iterate, the two options compose. Here is the same circuit with selective asinh compression of the exponential KCL row on top of the limiter:

Gsel(fu, u, p) = (fu[1] = asinh(fu[1]); nothing)
sol_both = solve(
    prob_c, NewtonRaphson(); precondition = Gsel, postcondition = H!, maxiters = 1000
)
sol_both.retcode, sol_both.stats.nsteps

(SciMLBase.ReturnCode.Success, 11)

Solver-state-aware correctors

A corrector always receives the solver cache as its fourth argument (analogous to PETSc's post-check receiving the SNES object); the correctors above simply ignore it. That argument is what staged limiting needs — Xyce-style simulators relax or disable limiting as the iteration proceeds:

H_staged! = function (up, uprev, p, cache)
    # `cache` is `nothing` for the initial-guess correction (no cache exists yet);
    # afterwards it is the solver cache, queried via the public accessors only
    if cache === nothing || NonlinearSolveBase.get_nsteps(cache) < 20
        up[2] = pnjlim(up[2], uprev[2], p.Vt, vcrit)
    end
    return nothing
end

sol_staged = solve(prob_c, NewtonRaphson(); postcondition = H_staged!, maxiters = 1000)
sol_staged.retcode, sol_staged.stats.nsteps

(SciMLBase.ReturnCode.Success, 11)

Treat the cache as read-only through its public accessors — NonlinearSolveBase.get_u, NonlinearSolveBase.get_fu (the residual at the previous accepted iterate at the time the corrector runs), NonlinearSolveBase.get_nsteps, NonlinearSolveBase.get_abstol, and NonlinearSolveBase.get_reltol — everything else on the cache is internal and subject to change without notice.

Constraints and projections via postcondition

Because H may enforce state exactly, the same option covers projection-style corrections: clamping iterates into a physical domain (positivity for concentrations, saturations in [0, 1]) or pinning Dirichlet-type values so the remaining equations act as the condensed problem. For example, protecting a log from a Newton overshoot into the negative domain:

flog(u, p) = log.(u) .- p
Hpos(up, uprev, p, cache) = clamp.(up, 1.0e-8, Inf)

sol_pos = solve(NonlinearProblem(flog, [10.0], -2.0), NewtonRaphson(); postcondition = Hpos)
sol_pos.retcode, sol_pos.u

(SciMLBase.ReturnCode.Success, [0.1353352832366127])

Correctors on bounded problems

For simple box bounds, prefer the native lb/ub support described in the bound constraints tutorial. With a postcondition, the bounded default uses BoundedTrustRegion, which applies the corrector in the original coordinates. If you explicitly choose an algorithm such as NewtonRaphson, bounds are handled by reparameterizing the iterate. You can then choose which coordinates the corrector sees:

prob_bounded = NonlinearProblem(
    (u, p) -> [u[1] - 1, u[2]^2 - u[1] - 3], [5.0, 5.0];
    lb = [0.0, 0.0], ub = [10.0, 10.0]
)
Hpin(up, uprev, p, cache) = [1.0, up[2]]   # pin the first component to 1

sol_orig = solve(prob_bounded, NewtonRaphson(); postcondition = Hpin, maxiters = 100)
sol_transformed = solve(
    prob_bounded, NewtonRaphson(); maxiters = 100,
    postcondition = PostconditionSpecifier(
        Hpin; space = PostconditionSpace.Transformed
    )
)
sol_orig.u[1], sol_transformed.u[1]

(0.9999999999999998, 7.310585786300049)

By default (space = PostconditionSpace.Original) H acts on your original bounded variable: each iterate is mapped back through the bounds transform, corrected, and mapped forward again. A limiting rule written for a physical quantity — pnjlim on a junction voltage in volts, a saturation clamped into [0, 1] — therefore means what it says, and the pin above lands on 1. Under space = PostconditionSpace.Transformed the same corrector pins the unconstrained coordinate to 1, which is the physical value lb + (ub - lb) * logistic(1) ≈ 7.31.

One consequence of PostconditionSpace.Transformed mode is that the initial-guess correction is skipped: the initial guess is still in the original variable at the point where it would run.

A correction that lands exactly on a bound sits at infinity in the transformed variable. Rather than committing an infinite iterate, NonlinearSolve nudges it into the interior by a relative eps^(3/4) — the same nudge the bounds transform applies to u0 — so clamping onto a bound is safe.

Semantics and caveats

• Where H acts.postcondition is applied to accepted iterates. Line searches and trust regions evaluate their merit/reduction models at the unlimited trial points, matching PETSc post-check semantics. Limiting therefore pairs most naturally with plain NewtonRaphson; a trust region will fight a strongly active limiter (it converges, but slowly), and quasi-Newton secant updates can be degraded by aggressively clipped steps.
• Convergence theory. Near a root the limiter must deactivate (H → identity), recovering Newton's local quadratic convergence; pnjlim and projections satisfy this. The Jacobian intentionally does not chain through H — it is a corrector between steps, not part of the residual.
• Solver support. The first-order (NewtonRaphson, TrustRegion, LevenbergMarquardt, ...), quasi-Newton (Broyden, Klement, ...), and spectral (DFSane) families and poly-algorithms composed of them support postcondition; unsupported solvers (e.g. SimpleNonlinearSolve or external wrappers) report the unapplied corrector through the unsupported_postcondition verbosity toggle, which raises by default rather than silently ignoring it and can be turned down through SciMLLogging when that is deliberate. precondition is a problem transformation and works with every solver that consumes the problem function, including SimpleNonlinearSolve.

References

• P. R. Brune, M. G. Knepley, B. F. Smith, X. Tu, Composing Scalable Nonlinear Algebraic Solvers, SIAM Review 57(4), 2015.
• K. V. Aadithya, E. R. Keiter, T. Mei, Predictor/Corrector Newton-Raphson (PCNR): A Simple, Flexible, Scalable, Modular, and Consistent Replacement for Limiting in Circuit Simulation, Scientific Computing in Electrical Engineering, 2020.
• L. W. Nagel, SPICE2: A Computer Program to Simulate Semiconductor Circuits, UC Berkeley, 1975 (the original pnjlim).