Structured L-BFGS Methods

• Structured L-BFGS is a collection of limited-memory quasi-Newton methods that exploit problem-specific Hessian structure to enhance convergence in ill-conditioned, nonconvex optimization tasks.
• It integrates cheaply computed regularizer Hessians into a modified seed matrix, using two-loop recursions and diagonal or block-structured scaling to accurately capture curvature.
• Empirical studies show that structured L-BFGS methods reduce iteration counts and runtime, proving effective in large-scale applications like medical image registration and inverse problems.

Structured L-BFGS refers to a collection of limited-memory BFGS (L-BFGS) algorithms that exploit problem-specific structure in large-scale unconstrained optimization, particularly when the objective is a sum of two functions with distinct Hessian properties. These methods incorporate partial, often cheap, Hessian information—typically available from regularization terms or known block structures—directly into the quasi-Newton update and search direction calculation, yielding improved convergence speed and robustness on ill-conditioned and non-convex problems. Recent advances provide rigorous global and linear convergence results in non-convex Hilbert space settings and demonstrate substantial performance gains in domains such as large-scale inverse problems and medical image registration (Mannel et al., 2023, Mannel et al., 2024).

1. Problem Framework and Motivation

Structured L-BFGS methods address optimization problems of the form

$\min_{x \in \mathbb{R}^n} \Phi(x) = f(x) + R(x)$

where $f$ is a data-fidelity term with inaccessible or computationally expensive Hessian $\nabla^2 f(x)$, while $R$ is a regularizer with cheaply computable, structured, or sparse Hessian $\nabla^2 R(x)$. The overall Hessian, $\nabla^2 \Phi(x) = \nabla^2 f(x) + \nabla^2 R(x)$, thus has a structure that can be exploited for efficient optimization.

Typical applications include large-scale inverse problems such as medical image registration and scientific imaging, where the regularizer encodes known spatial or physical priors, and the data term is dense or costly. Exploiting $\nabla^2 R(x)$ as a partial surrogate for the true Hessian is crucial for accelerating convergence and improving search direction quality on ill-conditioned landscapes (Mannel et al., 2023).

2. Structured L-BFGS Algorithmic Methodology

2.1 Seed Matrix Construction

Classical L-BFGS initializes the quasi-Newton update recursion with a scalar times the identity: $B_k^{(0)} = \tau_k I$, where $\tau_k > 0$. Structured L-BFGS replaces this with

$B_k^{(0)} = \tau_k I + S_k$

where $f$0 is the regularizer Hessian (or a cheap surrogate), and $f$1 calibration uses a Barzilai–Borwein–type fit based on the secant equation involving $f$2 and $f$3. This ensures that curvature contributions from $f$4 precondition all memory updates (Mannel et al., 2023).

Recent enhancements further generalize $f$5 to a full diagonal matrix $f$6 (with $f$7), allowing more accurate modeling of data-fidelity curvature, in addition to the regularizer, leading to faster convergence (Mannel et al., 2024).

2.2 Two-Loop Recursion and Direction Computation

The core quasi-Newton search direction is built using the two-loop recursion. For each iteration:

1. Apply the classical or structured two-loop recursion to stored pairs $f$8 with $f$9.
2. Substitute the scalar inversion with a solve $\nabla^2 f(x)$0 (or $\nabla^2 f(x)$1 for diagonal-seeded variants), leveraging the sparse or structured operator. The inversion is never formed explicitly; only one linear solve per iteration is required, and the cost per iteration matches that of the linear solve plus $\nabla^2 f(x)$2 vector operations (for memory $\nabla^2 f(x)$3).
3. Safeguards ensure all iterates preserve positive definiteness and bounded condition numbers essential for global convergence (Mannel et al., 2023, Mannel et al., 2024).

2.3 Diagonal and Block-Structured Initiation

The diagonalization approach introduced in ROSE (Mannel et al., 2024) solves a per-coordinate least-squares system to set $\nabla^2 f(x)$4, with entries projected into a bounded interval for numerical stability. This produces a seed matrix that, even without explicit knowledge of $\nabla^2 f(x)$5, better fits the true local curvature than scalar scaling.

Block-structured L-BFGS extensions, relevant for deep neural networks, use Kronecker-factored or block-diagonal Hessian models, integrating structured curvature within each block, with per-block or per-layer L-BFGS updates for massive-scale models (Goldfarb et al., 2020).

3. Convergence Analysis and Theoretical Results

Structured L-BFGS methods achieve global convergence under mild conditions: $\nabla^2 f(x)$6 objective, bounded gradients, bounded $\nabla^2 f(x)$7, and line search satisfying Armijo or Wolfe criteria. The positivity of the seed matrix (and thus the Hessian approximations) is preserved using cautious updates and safeguarding on $\nabla^2 f(x)$8 or $\nabla^2 f(x)$9 scaling.

Under a Kurdyka–Łojasiewicz-type (or Polyak–Łojasiewicz, PL) inequality, which holds for a broad class of non-convex problems, the method attains Q-linear convergence in objective reduction and R-linear convergence for iterates and gradients (Mannel et al., 2023, Mannel et al., 2024):

$R$0

Local linear rate results hold whenever the function is strongly convex in a neighborhood of the limit point. These are, to date, the first global and linear convergence proofs for L-BFGS-type methods without global convexity in Hilbert spaces.

A plausible implication is that structured initialization not only improves practical performance but also strengthens theoretical guarantees in the non-convex and ill-posed regimes typical of inverse problems.

4. Algorithmic Variants and Practical Implementation

The core method admits several extensions:

• Diagonal-Scaling Seed (ROSE): Improves upon scalar scaling by solving for the diagonal entries of $R$1 using two distinct least-squares criteria (the "p-formula" and "g-formula"), each projecting into enforceable bounds. The linear system $R$2 is solved using preconditioned Krylov methods (e.g., MINRES). Early stopping in these inner solves can further trade slight accuracy for large computational gain, with no loss in global convergence (Mannel et al., 2024).
• Nonlinear Preconditioning: For problems with nonlinear preconditioners (e.g., alternating least squares in tensor decomposition), a generalized L-BFGS can be wrapped around the nonlinear mapping, substituting the gradient and secant difference with the preconditioned residual and its differences. This framework provides robust acceleration for slow-converging fixed-point iterations (Sterck et al., 2018).
• Compact Representations: Matrix-free representations with structured BFGS updates use compact factorizations tailored to the composite structure of the objective and facilitate efficient limited-memory truncation and initialization (2208.00057). These are essential for scaling to $R$3 dimensions.
• Trust-Region Extensions: Efficient recursions for solving $R$4 arise naturally in trust-region and barrier methods, maintaining $R$5 complexity (with $R$6 memory terms) using Sherman-Morrison-Woodbury-style updates (Erway et al., 2011).

The following table summarizes key structural variants:

Variant	Seed Matrix	Structure Exploited
TULIP	$R$7	Regularizer Hessian
ROSE	$R$8 (diagonal)	Diagonal + regularizer
Block/Kronecker	Block diagonal, per-layer	Layerwise curvature (DNNs)

5. Computational Complexity and Storage

The primary computational bottleneck is the structured linear solve involving $R$9 or $\nabla^2 R(x)$0, whose cost depends on the sparsity or factorizability of $\nabla^2 R(x)$1. For large sparse systems, per-iteration complexity is near-linear in $\nabla^2 R(x)$2, with total storage $\nabla^2 R(x)$3 for memory window $\nabla^2 R(x)$4.

Compared to classical L-BFGS:

• Structured seed: +1 linear solve per iteration (structured or preconditioned)
• Storage: identical at $\nabla^2 R(x)$5
• Overall scaling dramatically improves when $\nabla^2 R(x)$6 or $\nabla^2 R(x)$7 is sparse, block-diagonal, or admits fast solvers (Mannel et al., 2023, Mannel et al., 2024).

In applications with millions of variables, block-diagonal and Kronecker-factored approaches avoid ever storing dense matrices, making the method viable for extreme-scale optimization (Goldfarb et al., 2020).

6. Numerical Results and Applications

Empirical studies over large-scale, non-convex medical image registration tasks (n up to $\nabla^2 R(x)$8), using similarity metrics (SSD, MI, NGF) and regularizers (curvature, elasticity, hyperelasticity), demonstrate that adaptive structured L-BFGS methods:

• Solve $\nabla^2 R(x)$9 of problems within 1.5$\nabla^2 \Phi(x) = \nabla^2 f(x) + \nabla^2 R(x)$0 the best time, compared to $\nabla^2 \Phi(x) = \nabla^2 f(x) + \nabla^2 R(x)$1 for classical L-BFGS.
• Achieve lower target registration error (TRE) across almost all benchmark cases.
• Outperform S-LBFGS‐M/P in both run-time (2–3$\nabla^2 \Phi(x) = \nabla^2 f(x) + \nabla^2 R(x)$2 faster) and solution quality.

For ill-conditioned quadratic models, structured variants cut iteration counts by $\nabla^2 \Phi(x) = \nabla^2 f(x) + \nabla^2 R(x)$3 depending on the regularizer strength and yield search directions significantly closer (in angle and norm) to the Newton step.

ROSE with diagonal scaling achieves $\nabla^2 \Phi(x) = \nabla^2 f(x) + \nabla^2 R(x)$4 run-time reductions relative to previous structured methods, with comparable or better accuracy on imaging problems. Similar acceleration is reported in other composite quadratic and high-dimensional settings (Mannel et al., 2023, Mannel et al., 2024).

7. Extensions, Related Work, and Future Outlook

Structured L-BFGS generalizes to settings where only partial Hessian information is available:

• Tensor Decomposition and Nonlinear Preconditioning: L-BFGS may be used as a nonlinear accelerator for iterations with an inherent fixed-point structure, such as ALS or HOOI, yielding up to $\nabla^2 \Phi(x) = \nabla^2 f(x) + \nabla^2 R(x)$5 time savings over standalone methods (Sterck et al., 2018).
• Superlinear Convergence: Limited-memory Greedy L-BFGS yields explicit non-asymptotic superlinear rates with appropriate displacement aggregation and memory-basis selection (Gao et al., 2023).
• Trust-Region and Interior-Point Methods: Fast solvers for $\nabla^2 \Phi(x) = \nabla^2 f(x) + \nabla^2 R(x)$6 are fundamental for inner linear systems in constrained optimization and large-scale interior methods, provided by efficient structured L-BFGS recursions (Erway et al., 2011).
• Compact Representations: Compact structured BFGS formulas with block or diagonal initialization facilitate memory scaling and efficient line-search computation for structured objectives in scientific imaging and control (2208.00057).

Open directions include unifying the various diagonal and block-structured seed strategies, adaptive memory management, and systematizing theory for non-convex and non-Euclidean settings. These advances expand the scope of L-BFGS solvers to new classes of inverse, scientific, and machine learning problems not tractable by first-order or classical quasi-Newton methods (Mannel et al., 2023, Mannel et al., 2024, 2208.00057).