Descent Subroutines

NewtonDescent(; linsolve = nothing)

Compute the descent direction as $J δu = -fu$. For non-square Jacobian problems, this is commonly referred to as the Gauss-Newton Descent.

See also Dogleg, SteepestDescent, DampedNewtonDescent.

SteepestDescent(; linsolve = nothing)

Compute the descent direction as $δu = -Jᵀfu$. The linear solver and preconditioner are only used if J is provided in the inverted form.

See also Dogleg, NewtonDescent, DampedNewtonDescent.

DampedNewtonDescent(;
    linsolve = nothing, initial_damping, damping_fn,
    min_norm_mode = :auto
)

A Newton descent algorithm with damping. The damping factor is computed using the damping_fn function. The descent direction is computed as $(JᵀJ + λDᵀD) δu = -fu$. For non-square Jacobians, we default to solving for Jδx = -fu and √λ⋅D δx = 0 simultaneously. If the linear solver can't handle non-square matrices, we use the normal form equations $(JᵀJ + λDᵀD) δu = Jᵀ fu$. Note that this factorization is often the faster choice, but it is not as numerically stable as the least squares solver.

For underdetermined systems (more unknowns than equations, length(u) > length(fu)) with a normal form damping function, we default to a minimum-norm formulation: solve $(JJᵀ + λD̃ᵀD̃) z = -fu$ and set $δu = Jᵀz$, where $D̃ᵀD̃$ is the damping the damping_fn produces for the JJᵀ system — for Levenberg-Marquardt this is the running elementwise maximum of diag(JJᵀ), so the damping scales with $‖J‖²$ and the step is invariant to rescaling the problem. As $λ → 0$ this recovers the minimum-norm step solving the linearized equations, and it keeps the linear system small (m × m where m = length(fu) instead of n × n). Like the normal form equations, this formulation squares the conditioning of the Jacobian, so it trades some numerical robustness for speed; pass min_norm_mode = :disabled to use the QR-based least squares formulation instead. JJᵀ + λD̃ᵀD̃ is symmetric positive definite, so when no linsolve is specified this mode defaults to a Cholesky factorization (when LinearSolve is loaded).

The damping factor returned must be a non-negative number.

Keyword Arguments

• initial_damping: the initial damping factor to use

• damping_fn: the function to use to compute the damping factor. This must satisfy the NonlinearSolveBase.AbstractDampingFunction interface.

• min_norm_mode: controls the minimum-norm formulation for underdetermined systems:

  • :auto (default): use the minimum-norm formulation for underdetermined systems when the damping function returns normal form damping
  • :minimum_norm: force the minimum-norm formulation regardless of system dimensions
  • :disabled: never use the minimum-norm formulation

Dogleg(; linsolve = nothing)

Switch between Newton's method and the steepest descent method depending on the size of the trust region. The trust region is specified via keyword argument trust_region to solve!.

See also SteepestDescent, NewtonDescent, DampedNewtonDescent.

MoreTrustRegionDescent(;
    linsolve = nothing, scaling = TrustRegionScaling.None, min_damping_D = 1e-8
)

Compute the descent direction by solving the trust-region subproblem min ‖J δu + fu‖ subject to ‖D δu‖ ≤ Δ nearly exactly. A safeguarded Newton iteration finds the damping parameter λ such that the least-squares problem

min ‖ [J; √λ D] δu - [-fu; 0] ‖

— equivalent to the damped normal equations (JᵀJ + λDᵀD) δu = -Jᵀfu — has solution with ‖D δu‖ = Δ; if the Gauss-Newton step already lies inside the region it is the subproblem solution and λ = 0 is returned directly. This is the algorithm of Moré (MINPACK lmpar), also described in Nocedal & Wright §10.3.

Solving the augmented least-squares system keeps the condition number of J unsquared and handles rank-deficient Jacobians natively: [J; √λ D] has full column rank for every λ > 0. When the chosen linsolve only supports square systems the equivalent damped normal equations are used instead, and for matrix-free (operator) Jacobians the augmented system is applied as a stacked operator, which Krylov least-squares solvers handle without forming J.

Unlike Dogleg, which follows a two-piece polygonal approximation of the solution curve, this descent follows the true solution δu(λ), which is substantially more robust on ill-conditioned least-squares problems. Pair it with RadiusUpdateSchemes.More to get Moré's step-following radius update as well.

Keyword Arguments

• linsolve: the linear solver used for the subproblem solves. On the default dense path the Jacobian is factorized once per iteration and the damping iteration runs against the stored factors (MINPACK lmpar/qrsolv), so no linsolve is invoked; on wide Jacobians (2m ≤ n) the iteration instead works on the m × m Gram system (J D⁻² Jᵀ + λI) w = -fu. An explicit choice is honored through the rectangular augmented system [J; √λD] p = [-fu; 0], so it should handle least-squares problems; solvers requiring square systems such as LUFactorization are routed through the normal equations automatically.
• scaling: the diagonal scaling matrix D, a TrustRegionScaling. TrustRegionScaling.None uses D = I; TrustRegionScaling.Jacobian uses Moré's scaling Dᵢᵢ = max(Dᵢᵢ, ‖J[:, i]‖), which never decreases across iterations and makes the trust region scale-covariant. TrustRegionScaling.Auto engages that same scaling only when the problem needs it: if the nonzero column norms of the first Jacobian span a ratio maxⱼ‖J[:, j]‖ / minⱼ‖J[:, j]‖ above 20 the solve proceeds exactly as Jacobian, and otherwise exactly as None. The decision is taken once per solve (re-decided after reinit!), and a matrix-free Jacobian simply never engages. TrustRegionScaling.Jacobian requires a concrete Jacobian.
• min_damping_D: lower bound for the entries of DᵀD under Jacobian or active Auto scaling.

TrustRegionSubproblem

Choice of trust-region subproblem solver, selected with TrustRegion(subproblem = ...):

• TrustRegionSubproblem.More: solve min ‖J δu + fu‖ subject to ‖D δu‖ ≤ Δ nearly exactly via Moré's safeguarded iteration on the damping parameter (MINPACK lmpar), through MoreTrustRegionDescent.
• TrustRegionSubproblem.Dogleg: the classical two-piece polygonal dogleg approximation to the subproblem solution curve, through Dogleg.

A custom AbstractDescentDirection can also be passed directly to use a user-defined subproblem solver.

TrustRegionScaling

Choice of the diagonal scaling D in the trust-region subproblem min ‖J δu + fu‖ subject to ‖D δu‖ ≤ Δ, selected with MoreTrustRegionDescent(scaling = ...):

• TrustRegionScaling.None: D = I.
• TrustRegionScaling.Jacobian: Moré's scaling Dᵢᵢ = max(Dᵢᵢ, ‖J[:, i]‖), which never decreases across iterations and makes the trust region scale-covariant. Requires a concrete Jacobian.
• TrustRegionScaling.Auto: engages Jacobian scaling only when the problem needs it — the decision is taken once per solve from the first Jacobian and a matrix-free Jacobian never engages.

GeodesicAcceleration(; descent, finite_diff_step_geodesic, α)

Uses the descent algorithm to compute the velocity and acceleration terms for the geodesic acceleration method. The velocity and acceleration terms are then combined to compute the descent direction.

This method in its current form was developed for LevenbergMarquardt. Performance for other methods are not theoretically or experimentally verified.

Keyword Arguments

• descent: the descent algorithm to use for computing the velocity and acceleration.
• finite_diff_step_geodesic: the step size used for finite differencing used to calculate the geodesic acceleration. Defaults to 0.1 which means that the step size is approximately 10% of the first-order step. See Section 3 of [1].
• α: a factor that determines if a step is accepted or rejected. To incorporate geodesic acceleration as an addition to the Levenberg-Marquardt algorithm, it is necessary that acceptable steps meet the condition $\frac{2||a||}{||v||} \le \alpha_{\text{geodesic}}$, where $a$ is the geodesic acceleration, $v$ is the Levenberg-Marquardt algorithm's step (velocity along a geodesic path) and α_geodesic is some number of order 1. For most problems α_geodesic = 0.75 is a good value but for problems where convergence is difficult α_geodesic = 0.1 is an effective choice. Defaults to 0.75. See Section 3 of Transtrum and Sethna [1].