Optimization

Optimization algorithms for Aurora-GLM.

class aurora.core.optimization.OptimizationResult(x, fun, grad=None, hess=None, success=True, message='Optimization terminated successfully', nit=0, nfev=0, njev=0, nhev=0, backtrack_iterations=0, condition_number=None)[source]

Bases: object

Container for optimization results.

Variables:

  • x (Array) – Optimal parameter vector found by the optimiser.

  • fun (Scalar) – Objective function value at x.

  • grad (Array or None) – Gradient of the objective at x (if available).

  • hess (Array or None) – Hessian of the objective at x (if available).

  • success (bool, default=True) – Whether the optimiser converged successfully.

  • message (str) – Human-readable description of the termination reason.

  • nit (int, default=0) – Number of iterations performed.

  • nfev (int, default=0) – Number of objective function evaluations.

  • njev (int, default=0) – Number of gradient evaluations.

  • nhev (int, default=0) – Number of Hessian evaluations.

  • backtrack_iterations (int, default=0) – Number of backtracking (step-halving) reductions performed.

  • condition_number (float or None) – Maximum observed condition number of the weighted normal equations (IRLS) or the Hessian (Newton). None when not monitored.

Parameters:
x: ndarray | Any
fun: int | float | complex
grad: ndarray | Any | None = None
hess: ndarray | Any | None = None
success: bool = True
message: str = 'Optimization terminated successfully'
nit: int = 0
nfev: int = 0
njev: int = 0
nhev: int = 0
backtrack_iterations: int = 0
condition_number: float | None = None
class aurora.core.optimization.Optimizer[source]

Bases: ABC

Abstract base class for Aurora optimizers.

abstractmethod minimize(loss_fn, init_params, *, callback=None, max_iter=100, tol=1e-6, **kwargs)[source]

Execute the optimization routine.

Parameters:
Return type:

OptimizationResult

aurora.core.optimization.newton_raphson(loss_fn, init_params, *, backend=None, args=(), kwargs=None, max_iter=100, tol=1e-6, callback=None)[source]

Run the Newton-Raphson method with automatic differentiation support.

Minimises loss_fn using the classical Newton-Raphson algorithm. The Hessian is computed via automatic differentiation (PyTorch or JAX) when available, falling back to central finite differences otherwise.

Parameters:

  • loss_fn (callable) – Objective function to minimise. Must accept a parameter array as its first argument and return a scalar.

  • init_params (array-like, shape (p,)) – Initial parameter values.

  • backend (object, optional) – Array backend (NumPy, PyTorch, JAX). Defaults to JAX when available.

  • args (tuple, default=()) – Additional positional arguments forwarded to loss_fn.

  • kwargs (dict, optional) – Additional keyword arguments forwarded to loss_fn.

  • max_iter (int, default=100) – Maximum number of Newton iterations.

  • tol (float, default=1e-6) – Convergence tolerance on gradient norm and step size.

  • callback (callable, optional) – Function called as callback(iteration, params, loss_value) after each iteration.

Returns:

Result containing the solution, final loss, gradient, convergence status, and diagnostic counts.

Return type:

OptimizationResult

Notes

Convergence is quadratic near the solution when the Hessian is positive definite. The method may diverge from poor starting points; consider using modified_newton for robustness.

Raises:

numpy.linalg.LinAlgError – (Returned, not raised) When the Hessian is singular at the current iterate, the result is returned with success=False.

Parameters:
Return type:

OptimizationResult

See also

modified_newton

Newton-Raphson with Levenberg-Marquardt regularisation.

irls

IRLS for GLMs (preferred for canonical-link models).

lbfgs

L-BFGS quasi-Newton method for large-scale problems.

aurora.core.optimization.modified_newton(loss_fn, init_params, *, backend=None, args=(), kwargs=None, max_iter=100, tol=1e-6, lambda_init=1e-3, lambda_factor=10.0, lambda_max=1e8, callback=None)[source]

Modified Newton-Raphson with Levenberg-Marquardt regularization.

Mathematical Framework

When the Hessian H is indefinite or ill-conditioned, the standard Newton step p = -H⁻¹g may not be a descent direction. Modified Newton adds regularization:

(H + λI) p = -g

where λ > 0 ensures positive definiteness.

Levenberg-Marquardt Strategy

The regularization parameter λ is adjusted adaptively:

  1. Start with small λ: Allow near-Newton steps when H is well-behaved

  2. Increase λ if step rejected: When function doesn’t decrease sufficiently, multiply λ by lambda_factor to trust gradient more than curvature

  3. Decrease λ if step accepted: After successful steps, reduce λ to allow faster convergence via better Hessian approximation

Interpretation: - λ → 0: Pure Newton (quadratic convergence near solution) - λ → ∞: Steepest descent (slow but robust)

The algorithm smoothly interpolates between these extremes.

Guaranteeing Positive Definiteness

For symmetric H with eigenvalues λ₁ ≤ λ₂ ≤ … ≤ λₙ: - H + λI has eigenvalues λ₁ + λ, λ₂ + λ, …, λₙ + λ - If λ > -λ₁ (where λ₁ is the smallest eigenvalue), H + λI is positive definite

In practice, we use the Cholesky factorization attempt: - If Cholesky succeeds: H + λI is positive definite - If Cholesky fails: Increase λ and retry

Convergence Properties

Theorem (Dennis & Schnabel, 1996): Modified Newton with adaptive λ converges globally for any starting point if: 1. f is bounded below 2. ∇f is Lipschitz continuous 3. f is twice continuously differentiable

Local rate: Near the solution where H is positive definite: - If λ → 0 sufficiently fast: Quadratic convergence - If λ bounded away from 0: Linear convergence

param loss_fn:

Objective function to minimize

type loss_fn:

callable

param init_params:

Initial parameter values

type init_params:

array

param backend:

Array backend (NumPy, PyTorch, JAX)

type backend:

object, optional

param args:

Additional positional arguments for loss_fn

type args:

tuple

param kwargs:

Additional keyword arguments for loss_fn

type kwargs:

dict, optional

param max_iter:

Maximum number of iterations

type max_iter:

int, default=100

param tol:

Convergence tolerance for gradient norm

type tol:

float, default=1e-6

param lambda_init:

Initial regularization parameter

type lambda_init:

float, default=1e-3

param lambda_factor:

Factor to increase/decrease λ

type lambda_factor:

float, default=10.0

param lambda_max:

Maximum allowed λ (switches to pure gradient descent)

type lambda_max:

float, default=1e8

param callback:

Function called after each iteration

type callback:

callable, optional

returns:

Contains solution, final function value, convergence status

rtype:

OptimizationResult

References

[1] Levenberg, K. (1944). “A method for the solution of certain non-linear

problems in least squares.” Quarterly of Applied Mathematics, 2(2), 164-168.

[2] Marquardt, D. W. (1963). “An algorithm for least-squares estimation of

nonlinear parameters.” Journal of SIAM, 11(2), 431-441.

[3] Dennis, J. E., & Schnabel, R. B. (1996). Numerical Methods for

Unconstrained Optimization and Nonlinear Equations. SIAM.

[4] Nocedal, J., & Wright, S. J. (2006). Numerical Optimization (2nd ed.).

Springer. Chapter 4: Trust-Region Methods.

Examples

>>> def rosenbrock(x):
...     return (1 - x[0])**2 + 100*(x[1] - x[0]**2)**2
>>> result = modified_newton(rosenbrock, [-1.0, 1.0])
>>> print(result.x)  # Should be close to [1, 1]
Parameters:
Return type:

OptimizationResult

aurora.core.optimization.irls(loss_fn, init_params, *, backend=None, args=(), kwargs=None, max_iter=100, tol=1e-6, callback=None, design_matrix=None, response=None, link=None, variance_fn=None, offset=None)[source]

Run the IRLS procedure for generalized linear models.

This function automatically detects sparse design matrices and uses optimized sparse operations when beneficial.

Parameters:

  • loss_fn (callable) – Loss function taking (params, *args, **kwargs)

  • init_params (array-like (p,)) – Initial parameter estimates

  • backend (Backend, optional) – Computational backend (numpy, pytorch, jax)

  • args (tuple) – Positional arguments for loss_fn

  • kwargs (dict) – Keyword arguments for loss_fn

  • max_iter (int, default=100) – Maximum IRLS iterations

  • tol (float, default=1e-6) – Convergence tolerance on parameter change

  • callback (callable, optional) – Callback function(iteration, params, loss)

  • design_matrix (array-like (n, p)) – Design matrix X. May be dense or scipy.sparse.

  • response (array-like (n,)) – Response vector y

  • link (LinkFunction) – Link function with inverse() and derivative()

  • variance_fn (callable) – Variance function V(μ)

  • offset (array-like (n,), optional) – Known offset term

Returns:

Optimization result containing: - x: final parameters - fun: final loss value - success: convergence status - nit: number of iterations - message: status message

Return type:

OptimizationResult

Notes

For sparse design matrices (scipy.sparse), the algorithm automatically uses sparse matrix operations:

  1. Matrix-vector products X @ β use sparse BLAS

  2. Normal equations X^T W X are formed efficiently

  3. The p × p system is solved with Cholesky factorization

This is efficient when: - X is sparse (> 90% zeros) - Number of predictors p is moderate (< 10,000) - X^T W X fits in memory as dense

For very large sparse X^T W X, iterative solvers (CG) would be needed.

aurora.core.optimization.lbfgs(loss_fn, init_params, *, backend=None, args=(), kwargs=None, max_iter=100, tol=1e-6, m=10, line_search='strong-wolfe', callback=None, grad_fn=None)[source]

Run L-BFGS (Limited-memory BFGS) quasi-Newton optimization.

L-BFGS is a quasi-Newton method that approximates Newton’s method without explicitly forming or storing the full Hessian matrix. It uses the two-loop recursion algorithm to compute the search direction from the m most recent curvature pairs (s, y) and supports multi-backend computation (NumPy, PyTorch, and JAX).

Automatic differentiation is used when available (PyTorch / JAX); otherwise, central finite differences are used for gradient computation.

Parameters:

  • loss_fn (callable) – Objective function to minimise. Must accept the parameter vector as its first argument and return a scalar.

  • init_params (array-like of shape (p,)) – Initial parameter values.

  • backend (object, optional) – Computational backend (NumPy, PyTorch, or JAX). If None, defaults to the JAX backend.

  • args (tuple, default=()) – Additional positional arguments forwarded to loss_fn.

  • kwargs (dict, optional) – Additional keyword arguments forwarded to loss_fn.

  • max_iter (int, default=100) – Maximum number of L-BFGS iterations.

  • tol (float, default=1e-6) – Convergence tolerance on the gradient norm. The algorithm stops when ||grad|| < tol.

  • m (int, default=10) – Number of past curvature correction pairs (s, y) to store. Larger values give a better Hessian approximation at the cost of more memory (O(m * p)).

  • line_search ({'strong-wolfe', 'backtracking'}, default='strong-wolfe') – Line search method. 'strong-wolfe' satisfies the Armijo and curvature conditions (recommended for quasi-Newton methods). 'backtracking' uses a simpler Armijo-only backtracking search.

  • callback (callable, optional) – Function called after each iteration with signature (iteration, params, loss_value).

  • grad_fn (callable, optional) – Gradient function with the same signature as loss_fn that returns the gradient vector. When provided, this is used directly instead of calling backend.grad(loss_fn). This is useful for backends that do not support automatic differentiation (e.g. NumPy).

Returns:

Optimisation result containing:

  • x : final parameter vector

  • fun : final objective value

  • grad : final gradient vector

  • success : whether convergence was achieved

  • nit : number of iterations performed

  • nfev : number of function evaluations

  • njev : number of gradient evaluations

  • message : human-readable status message

Return type:

OptimizationResult

Notes

The algorithm maintains positive definiteness of the inverse Hessian approximation by discarding curvature pairs that violate the curvature condition y^T s > 0.

Examples

>>> from aurora.core.optimization import lbfgs
>>> result = lbfgs(lambda x: ((x - 3) ** 2).sum(), np.zeros(5))
>>> result.x  # Should be close to [3, 3, 3, 3, 3]

See also

newton_raphson

Classical Newton-Raphson method.

modified_newton

Newton-Raphson with Levenberg-Marquardt regularisation.

irls

IRLS for GLMs (preferred for canonical-link models).

aurora.core.optimization.optimize(loss_fn, init_params, method='lbfgs', backend='jax', **kwargs)[source]

Unified optimization interface.

Parameters:
aurora.core.optimization.solve_sparse_penalized_ls(B, z, weights, penalty, lambda_, method='auto', tol=1e-8, maxiter=None)[source]

Solve sparse penalized least squares: (B^T W B + λS)β = B^T W z.

Parameters:

  • B (array or sparse matrix, shape (n, k)) – Basis matrix (dense or sparse CSR format).

  • z (ndarray, shape (n,)) – Working response vector (pseudo-data).

  • weights (ndarray, shape (n,)) – Diagonal weight matrix entries (W = diag(weights)).

  • penalty (array or sparse matrix, shape (k, k)) – Penalty matrix S (typically sparse band matrix).

  • lambda (float) – Smoothing parameter (λ ≥ 0). Set to 0 for unpenalized LS.

  • method ({'auto', 'direct', 'cg', 'minres'}, default='auto') – Solver method: - ‘auto’: Choose based on problem size and sparsity - ‘direct’: Sparse Cholesky factorization (best for k ≤ 100) - ‘cg’: Conjugate Gradient (best for large k, well-conditioned) - ‘minres’: Minimum Residual (for nearly-singular problems)

  • tol (float, default=1e-8) – Convergence tolerance for iterative methods.

  • maxiter (int, optional) – Maximum iterations for iterative methods (default: k).

  • lambda_ (float)

Returns:

  • beta (ndarray, shape (k,)) – Solution vector β.

  • info (dict) – Solver information: - ‘method’: actual method used - ‘success’: whether solver converged - ‘iterations’: number of iterations (for iterative methods) - ‘residual’: final residual norm (for iterative methods)

Raises:

  • ImportError – If scipy is not available.

  • ValueError – If inputs have incompatible shapes or invalid method.

Return type:

tuple[ndarray, dict[str, Any]]

Notes

Problem: Minimize ||W^{1/2}(z - Bβ)||² + λ β^T S β

Normal equations: (B^T W B + λS) β = B^T W z

Method selection (‘auto’): - If k ≤ 50: use ‘direct’ (fast for small problems) - If k > 50 and λ > 1e-6: use ‘cg’ (well-conditioned) - If k > 50 and λ ≤ 1e-6: use ‘minres’ (may be ill-conditioned)

Sparsity: - If B is scipy.sparse, uses sparse operations throughout - If B is dense, converts to sparse if beneficial (k > 20 and B has <20% non-zeros) - Penalty S is assumed sparse (band matrix for difference penalties)

Preconditioning: Iterative methods use diagonal (Jacobi) preconditioning: M = diag(C) This improves convergence by approximately factor of 2-5×.

Examples

>>> # Solve with sparse B-spline basis
>>> from aurora.smoothing.splines import BSplineBasis
>>> x = np.linspace(0, 10, 100)
>>> knots = BSplineBasis.create_knots(x, n_basis=15, degree=3)
>>> basis = BSplineBasis(knots, degree=3)
>>> B_sparse = basis.basis_matrix(x, sparse=True)  # Sparse CSR matrix
>>>
>>> # Simulate data
>>> np.random.seed(42)
>>> z = np.sin(x) + 0.1 * np.random.randn(100)
>>> weights = np.ones(100)
>>>
>>> # Penalty matrix (second-order differences)
>>> S = basis.penalty_matrix(order=2)
>>>
>>> # Solve with automatic method selection
>>> beta, info = solve_sparse_penalized_ls(
...     B_sparse, z, weights, S, lambda_=0.1, method='auto'
... )
>>> print(f"Method: {info['method']}, Success: {info['success']}")
Method: direct, Success: True
>>>
>>> # Predictions
>>> fitted = B_sparse @ beta

IRLS

aurora.core.optimization.irls.irls(loss_fn, init_params, *, backend=None, args=(), kwargs=None, max_iter=100, tol=1e-6, callback=None, design_matrix=None, response=None, link=None, variance_fn=None, offset=None)[source]

Run the IRLS procedure for generalized linear models.

This function automatically detects sparse design matrices and uses optimized sparse operations when beneficial.

Parameters:

  • loss_fn (callable) – Loss function taking (params, *args, **kwargs)

  • init_params (array-like (p,)) – Initial parameter estimates

  • backend (Backend, optional) – Computational backend (numpy, pytorch, jax)

  • args (tuple) – Positional arguments for loss_fn

  • kwargs (dict) – Keyword arguments for loss_fn

  • max_iter (int, default=100) – Maximum IRLS iterations

  • tol (float, default=1e-6) – Convergence tolerance on parameter change

  • callback (callable, optional) – Callback function(iteration, params, loss)

  • design_matrix (array-like (n, p)) – Design matrix X. May be dense or scipy.sparse.

  • response (array-like (n,)) – Response vector y

  • link (LinkFunction) – Link function with inverse() and derivative()

  • variance_fn (callable) – Variance function V(μ)

  • offset (array-like (n,), optional) – Known offset term

Returns:

Optimization result containing: - x: final parameters - fun: final loss value - success: convergence status - nit: number of iterations - message: status message

Return type:

OptimizationResult

Notes

For sparse design matrices (scipy.sparse), the algorithm automatically uses sparse matrix operations:

  1. Matrix-vector products X @ β use sparse BLAS

  2. Normal equations X^T W X are formed efficiently

  3. The p × p system is solved with Cholesky factorization

This is efficient when: - X is sparse (> 90% zeros) - Number of predictors p is moderate (< 10,000) - X^T W X fits in memory as dense

For very large sparse X^T W X, iterative solvers (CG) would be needed.

Newton-Raphson

aurora.core.optimization.newton.newton_raphson(loss_fn, init_params, *, backend=None, args=(), kwargs=None, max_iter=100, tol=1e-6, callback=None)[source]

Run the Newton-Raphson method with automatic differentiation support.

Minimises loss_fn using the classical Newton-Raphson algorithm. The Hessian is computed via automatic differentiation (PyTorch or JAX) when available, falling back to central finite differences otherwise.

Parameters:

  • loss_fn (callable) – Objective function to minimise. Must accept a parameter array as its first argument and return a scalar.

  • init_params (array-like, shape (p,)) – Initial parameter values.

  • backend (object, optional) – Array backend (NumPy, PyTorch, JAX). Defaults to JAX when available.

  • args (tuple, default=()) – Additional positional arguments forwarded to loss_fn.

  • kwargs (dict, optional) – Additional keyword arguments forwarded to loss_fn.

  • max_iter (int, default=100) – Maximum number of Newton iterations.

  • tol (float, default=1e-6) – Convergence tolerance on gradient norm and step size.

  • callback (callable, optional) – Function called as callback(iteration, params, loss_value) after each iteration.

Returns:

Result containing the solution, final loss, gradient, convergence status, and diagnostic counts.

Return type:

OptimizationResult

Notes

Convergence is quadratic near the solution when the Hessian is positive definite. The method may diverge from poor starting points; consider using modified_newton for robustness.

Raises:

numpy.linalg.LinAlgError – (Returned, not raised) When the Hessian is singular at the current iterate, the result is returned with success=False.

Parameters:
Return type:

OptimizationResult

See also

modified_newton

Newton-Raphson with Levenberg-Marquardt regularisation.

irls

IRLS for GLMs (preferred for canonical-link models).

lbfgs

L-BFGS quasi-Newton method for large-scale problems.

L-BFGS

aurora.core.optimization.lbfgs.lbfgs(loss_fn, init_params, *, backend=None, args=(), kwargs=None, max_iter=100, tol=1e-6, m=10, line_search='strong-wolfe', callback=None, grad_fn=None)[source]

Run L-BFGS (Limited-memory BFGS) quasi-Newton optimization.

L-BFGS is a quasi-Newton method that approximates Newton’s method without explicitly forming or storing the full Hessian matrix. It uses the two-loop recursion algorithm to compute the search direction from the m most recent curvature pairs (s, y) and supports multi-backend computation (NumPy, PyTorch, and JAX).

Automatic differentiation is used when available (PyTorch / JAX); otherwise, central finite differences are used for gradient computation.

Parameters:

  • loss_fn (callable) – Objective function to minimise. Must accept the parameter vector as its first argument and return a scalar.

  • init_params (array-like of shape (p,)) – Initial parameter values.

  • backend (object, optional) – Computational backend (NumPy, PyTorch, or JAX). If None, defaults to the JAX backend.

  • args (tuple, default=()) – Additional positional arguments forwarded to loss_fn.

  • kwargs (dict, optional) – Additional keyword arguments forwarded to loss_fn.

  • max_iter (int, default=100) – Maximum number of L-BFGS iterations.

  • tol (float, default=1e-6) – Convergence tolerance on the gradient norm. The algorithm stops when ||grad|| < tol.

  • m (int, default=10) – Number of past curvature correction pairs (s, y) to store. Larger values give a better Hessian approximation at the cost of more memory (O(m * p)).

  • line_search ({'strong-wolfe', 'backtracking'}, default='strong-wolfe') – Line search method. 'strong-wolfe' satisfies the Armijo and curvature conditions (recommended for quasi-Newton methods). 'backtracking' uses a simpler Armijo-only backtracking search.

  • callback (callable, optional) – Function called after each iteration with signature (iteration, params, loss_value).

  • grad_fn (callable, optional) – Gradient function with the same signature as loss_fn that returns the gradient vector. When provided, this is used directly instead of calling backend.grad(loss_fn). This is useful for backends that do not support automatic differentiation (e.g. NumPy).

Returns:

Optimisation result containing:

  • x : final parameter vector

  • fun : final objective value

  • grad : final gradient vector

  • success : whether convergence was achieved

  • nit : number of iterations performed

  • nfev : number of function evaluations

  • njev : number of gradient evaluations

  • message : human-readable status message

Return type:

OptimizationResult

Notes

The algorithm maintains positive definiteness of the inverse Hessian approximation by discarding curvature pairs that violate the curvature condition y^T s > 0.

Examples

>>> from aurora.core.optimization import lbfgs
>>> result = lbfgs(lambda x: ((x - 3) ** 2).sum(), np.zeros(5))
>>> result.x  # Should be close to [3, 3, 3, 3, 3]

See also

newton_raphson

Classical Newton-Raphson method.

modified_newton

Newton-Raphson with Levenberg-Marquardt regularisation.

irls

IRLS for GLMs (preferred for canonical-link models).

Sparse Solvers

Sparse linear system solvers for penalized regression.

This module provides efficient solvers for large-scale penalized least squares problems commonly arising in GAM and GAMM fitting when using sparse basis matrices.

Mathematical Framework

Penalized Least Squares Problem

The penalized least squares problem arises in GAM fitting:

minimize: ||W^{1/2}(z - Bβ)||² + λ β^T S β

over β

where: - z ∈ ℝⁿ is the working response (pseudo-data) - B ∈ ℝⁿˣᵏ is the basis matrix (sparse for B-splines) - W ∈ ℝⁿˣⁿ is the diagonal weight matrix - S ∈ ℝᵏˣᵏ is the penalty matrix (sparse band matrix for difference penalties) - λ ≥ 0 is the smoothing parameter - β ∈ ℝᵏ are the spline coefficients

Normal Equations

Taking derivatives and setting to zero yields the penalized normal equations:

(B^T W B + λS) β = B^T W z (1)

Define: - C = B^T W B + λS (coefficient matrix) - d = B^T W z (right-hand side)

Then: C β = d

Structure: - C is symmetric positive definite (if S is positive semidefinite and λ > 0) - C is sparse when both B and S are sparse - Dimension: k × k where k = number of basis functions (typically 10-50)

Sparse Matrix Exploitation

For B-splines of degree p: - B has O(np) non-zeros (at most p+1 per row) - S has O(k) non-zeros (band matrix, bandwidth ≈ m for m-th order differences) - B^T W B is sparse with bandwidth ≈ 2p - C = B^T W B + λS is sparse, symmetric, positive definite

Storage: - Dense: k² elements (e.g., 30² = 900) - Sparse CSR: O(kp) elements (e.g., 30×4 = 120) - 6-8× memory reduction for typical GAMs

Solver Selection

Three solver strategies are provided:

  1. Direct Sparse Cholesky (default for small-medium problems) - Algorithm: Cholesky factorization C = LL^T using sparse ordering - Complexity: O(k³) flops but with much smaller constant than dense - Memory: O(k²) but exploits sparsity pattern - Best for: k ≤ 100, multiple solves with same C - Implementation: scipy.sparse.linalg.spsolve with ordering=’natural’ or ‘amd’

  2. Conjugate Gradient (CG) - Algorithm: Iterative Krylov subspace method for symmetric positive definite C - Complexity: O(k²) per iteration, typically converges in O(k) iterations - Memory: O(k) working space - Best for: Large k > 100, C is well-conditioned - Implementation: scipy.sparse.linalg.cg

  3. MINRES (Minimum Residual) - Algorithm: Iterative method for symmetric (not necessarily positive definite) - Complexity: Similar to CG - Best for: Nearly-singular C (e.g., λ ≈ 0, rank-deficient S) - Implementation: scipy.sparse.linalg.minres

Preconditioning: For iterative solvers, we use diagonal (Jacobi) preconditioning:

M = diag(C)

This improves convergence rate by 2-5× for typical GAM problems.

Computational Complexity

### Dense Solver (naive) - Form C: O(n k² + k²) ≈ O(nk²) flops (dense matrix-matrix product B^T W B) - Cholesky: O(k³) flops - Solve: O(k²) flops - Total: O(nk² + k³)

For n = 1000, k = 30: ≈ 900M + 27K ≈ 900M flops

### Sparse Direct Solver - Form C: O(np² + kp) ≈ O(np²) flops (exploiting sparse B) - Sparse Cholesky: O(k × bandwidth²) ≈ O(kp²) flops - Solve: O(k × bandwidth) ≈ O(kp) flops - Total: O(np² + kp²)

For n = 1000, k = 30, p = 3: ≈ 9K + 270 ≈ 9K flops Speedup: ~100× over dense

### Sparse Iterative Solver (CG) - Form C: O(np² + kp) flops - Each CG iteration: O(kp) flops (sparse matrix-vector product) - Iterations: O(√κ) where κ is condition number - Total: O(np² + √κ × kp)

For well-conditioned problems: ~5-10 iterations suffice Speedup: ~200× over dense for large k

Numerical Stability

Well-conditioned cases (λ > 0, full-rank problem): - Condition number: κ(C) ≈ max(eigenvalues of C) / min(eigenvalues of C) - With penalty: κ(C) ≈ O(λ + h⁻²) where h = knot spacing - Direct and iterative methods both stable

Ill-conditioned cases (λ ≈ 0 or rank-deficient S): - C may be nearly singular - Direct Cholesky may fail or be numerically unstable - Use MINRES with tight tolerance or add ridge penalty (λ_ridge + λ)

Preconditioning: Diagonal preconditioning M = diag(C) reduces condition number by:

κ(M^{-1}C) ≈ √κ(C)

References

Sparse linear systems: - Davis, T. A. (2006). Direct Methods for Sparse Linear Systems. SIAM.

Penalized regression splines: - Wood, S. N. (2017). Generalized Additive Models: An Introduction with R (2nd ed.).

CRC Press. Chapter 4: Smoothing and Penalized Least Squares.

Conjugate Gradient: - Hestenes, M. R., & Stiefel, E. (1952). “Methods of conjugate gradients for solving

linear systems.” Journal of Research of the National Bureau of Standards, 49(6).

Cholesky for sparse matrices: - George, A., & Liu, J. W. (1981). *Computer Solution of Large Sparse Positive

Definite Systems*. Prentice-Hall.

See also

aurora.core.optimization.irls

IRLS algorithm for GLM fitting

aurora.smoothing.penalties

Penalty matrix construction

aurora.models.gam.fitting

GAM fitting using penalized splines

Notes

For detailed algorithmic complexity analysis and benchmarks, see the repository documentation and the sparse algorithm complexity analysis document.

When λ = 0 (unpenalized least squares), the problem reduces to ordinary least squares and may be rank-deficient if k > rank(B). Add a small ridge penalty (e.g., λ = 1e-6) for numerical stability.

aurora.core.optimization.sparse_solvers.solve_sparse_penalized_ls(B, z, weights, penalty, lambda_, method='auto', tol=1e-8, maxiter=None)[source]

Solve sparse penalized least squares: (B^T W B + λS)β = B^T W z.

Parameters:

  • B (array or sparse matrix, shape (n, k)) – Basis matrix (dense or sparse CSR format).

  • z (ndarray, shape (n,)) – Working response vector (pseudo-data).

  • weights (ndarray, shape (n,)) – Diagonal weight matrix entries (W = diag(weights)).

  • penalty (array or sparse matrix, shape (k, k)) – Penalty matrix S (typically sparse band matrix).

  • lambda (float) – Smoothing parameter (λ ≥ 0). Set to 0 for unpenalized LS.

  • method ({'auto', 'direct', 'cg', 'minres'}, default='auto') – Solver method: - ‘auto’: Choose based on problem size and sparsity - ‘direct’: Sparse Cholesky factorization (best for k ≤ 100) - ‘cg’: Conjugate Gradient (best for large k, well-conditioned) - ‘minres’: Minimum Residual (for nearly-singular problems)

  • tol (float, default=1e-8) – Convergence tolerance for iterative methods.

  • maxiter (int, optional) – Maximum iterations for iterative methods (default: k).

  • lambda_ (float)

Returns:

  • beta (ndarray, shape (k,)) – Solution vector β.

  • info (dict) – Solver information: - ‘method’: actual method used - ‘success’: whether solver converged - ‘iterations’: number of iterations (for iterative methods) - ‘residual’: final residual norm (for iterative methods)

Raises:

  • ImportError – If scipy is not available.

  • ValueError – If inputs have incompatible shapes or invalid method.

Return type:

tuple[ndarray, dict[str, Any]]

Notes

Problem: Minimize ||W^{1/2}(z - Bβ)||² + λ β^T S β

Normal equations: (B^T W B + λS) β = B^T W z

Method selection (‘auto’): - If k ≤ 50: use ‘direct’ (fast for small problems) - If k > 50 and λ > 1e-6: use ‘cg’ (well-conditioned) - If k > 50 and λ ≤ 1e-6: use ‘minres’ (may be ill-conditioned)

Sparsity: - If B is scipy.sparse, uses sparse operations throughout - If B is dense, converts to sparse if beneficial (k > 20 and B has <20% non-zeros) - Penalty S is assumed sparse (band matrix for difference penalties)

Preconditioning: Iterative methods use diagonal (Jacobi) preconditioning: M = diag(C) This improves convergence by approximately factor of 2-5×.

Examples

>>> # Solve with sparse B-spline basis
>>> from aurora.smoothing.splines import BSplineBasis
>>> x = np.linspace(0, 10, 100)
>>> knots = BSplineBasis.create_knots(x, n_basis=15, degree=3)
>>> basis = BSplineBasis(knots, degree=3)
>>> B_sparse = basis.basis_matrix(x, sparse=True)  # Sparse CSR matrix
>>>
>>> # Simulate data
>>> np.random.seed(42)
>>> z = np.sin(x) + 0.1 * np.random.randn(100)
>>> weights = np.ones(100)
>>>
>>> # Penalty matrix (second-order differences)
>>> S = basis.penalty_matrix(order=2)
>>>
>>> # Solve with automatic method selection
>>> beta, info = solve_sparse_penalized_ls(
...     B_sparse, z, weights, S, lambda_=0.1, method='auto'
... )
>>> print(f"Method: {info['method']}, Success: {info['success']}")
Method: direct, Success: True
>>>
>>> # Predictions
>>> fitted = B_sparse @ beta

Result

class aurora.core.optimization.result.OptimizationResult(x, fun, grad=None, hess=None, success=True, message='Optimization terminated successfully', nit=0, nfev=0, njev=0, nhev=0, backtrack_iterations=0, condition_number=None)[source]

Bases: object

Container for optimization results.

Variables:

  • x (Array) – Optimal parameter vector found by the optimiser.

  • fun (Scalar) – Objective function value at x.

  • grad (Array or None) – Gradient of the objective at x (if available).

  • hess (Array or None) – Hessian of the objective at x (if available).

  • success (bool, default=True) – Whether the optimiser converged successfully.

  • message (str) – Human-readable description of the termination reason.

  • nit (int, default=0) – Number of iterations performed.

  • nfev (int, default=0) – Number of objective function evaluations.

  • njev (int, default=0) – Number of gradient evaluations.

  • nhev (int, default=0) – Number of Hessian evaluations.

  • backtrack_iterations (int, default=0) – Number of backtracking (step-halving) reductions performed.

  • condition_number (float or None) – Maximum observed condition number of the weighted normal equations (IRLS) or the Hessian (Newton). None when not monitored.

Parameters:
x: ndarray | Any
fun: int | float | complex
grad: ndarray | Any | None = None
hess: ndarray | Any | None = None
success: bool = True
message: str = 'Optimization terminated successfully'
nit: int = 0
nfev: int = 0
njev: int = 0
nhev: int = 0
backtrack_iterations: int = 0
condition_number: float | None = None