Extended-Krylov-subspace methods for trust-region and norm-regularization subproblems

1. Problem Statement

The paper addresses the efficient solution of Trust-Region (TR) and Norm-Regularization subproblems, which are fundamental components in large-scale nonlinear optimization algorithms.

• Trust-Region Subproblem (TRS): Given a symmetric Hessian matrix $A \in \mathbb{R}^{n \times n}$, a vector $b \in \mathbb{R}^{n}$, and a radius $\Delta > 0$, the goal is to find:
  $$x^* = \arg \min_{x \in \mathbb{R}^n} \left( \frac{1}{2}x^T A x - b^T x \right) \quad \text{subject to} \quad \|x\| \le \Delta$$
• Norm-Regularization Subproblem: Similar to TRS but with an objective function including a regularization term $\frac{\rho}{r}\|x\|^r$ (where $r \ge 2$).

Key Challenge:
Existing methods generally fall into two categories:

1. Factorization-based methods (e.g., Moré-Sorensen): Require multiple matrix factorizations (e.g., Cholesky or LDLT) as the shift parameter $\sigma$ changes. This is computationally expensive for large $n$.
2. Krylov-subspace methods (e.g., GLTR): Use standard Krylov subspaces ($\text{span}\{b, Ab, A^2b, \dots\}$) and avoid factorizations but may require many iterations, especially when the solution lies near the boundary of the trust region or when the shift $\sigma$ is close to the smallest eigenvalue of $A$.

The paper posits that the set of all possible solutions $x(\sigma)$ for varying $\sigma$ lies in a very low-dimensional subspace, but standard Krylov methods do not efficiently capture this structure for all $\sigma$.

2. Methodology: The TREK and NREK Algorithms

The authors propose a new approach based on Extended Krylov Subspaces (EKS).

A. Theoretical Insight: Low-Dimensional Manifold

The authors demonstrate that the solution manifold $\mathcal{M}(A, b) = \{x(\sigma) \mid \sigma \ge \sigma_{\min}\}$ is approximately spanned by a very small number of vectors (often $< 40$ even for $n=10^4$).

• They analyze the matrix $S$ formed by solutions at discrete shifts, showing it behaves like a generalized Cauchy matrix.
• Using results from Beckermann and Townsend, they prove that the $\epsilon$-rank of this matrix is logarithmically dependent on the condition number and precision, implying a very low effective rank.

B. The Extended Krylov Subspace

To capture the solution manifold efficiently, the method constructs a subspace that includes both positive and negative powers of $A$:
$$\mathcal{K}_{2m}(A, b) = \text{span}\{b, A^{-1}b, Ab, A^{-2}b, \dots, A^{-m}b\}$$
This is motivated by the fact that:

• For large $\sigma$ (interior solutions), $x(\sigma) \approx (A+\sigma I)^{-1}b$ can be approximated by polynomials in $A^{-1}$.
• For small $\sigma$ (boundary solutions), it can be approximated by polynomials in $A$.
• The EKS captures both regimes simultaneously.

C. The EKS Algorithm (Algorithm 4.1)

The authors adapt the Jagels-Reichel algorithm to generate an orthonormal basis $V_{2m+1}$ for the extended subspace.

• Process: It iteratively computes vectors $v_k$ and $v_{-k}$ by applying $A$ and $A^{-1}$ to the current basis vectors and orthogonalizing them.
• Cost: Requires only one matrix factorization (to solve systems with $A$) at the beginning. Subsequent steps involve only matrix-vector products ($Av$) and triangular solves ($A^{-1}v$).
• Structure: The projected Hessian $P = V^T A V$ becomes a pentadiagonal symmetric matrix, allowing for very cheap updates and eigenvalue computations.

D. The TREK Algorithm (Trust-Region)

1. Initialization: Check for the "interior case" (where the unconstrained minimizer satisfies the radius constraint). If not, compute $A^{-1}b$ (one factorization).
2. Iteration:
   • Expand the EKS basis $V_k$.
   • At each step, solve the small projected subproblem (dimension $\approx 2k$) using the spectral decomposition of the pentadiagonal $P_k$.
   • Compute the residual norm $\|r_k\|$ using the recurrence relations derived from the EKS structure.
   • Stopping Criterion: Terminate when the residual is below a tolerance $\epsilon$.
3. Handling Indefinite Matrices: If $A$ is not positive definite, the algorithm shifts $A$ by a Gershgorin-based shift $\sigma_S$ to make it positive definite, solves the shifted problem, and maps the solution back.

E. The NREK Algorithm (Regularization)

The same framework is applied to norm-regularization problems. The only difference is that the subproblem minimization includes the regularization term $\frac{\rho}{r}\|y\|^r$ instead of a hard constraint $\|y\| \le \Delta$.

3. Key Contributions

1. Low-Rank Approximation Proof: Theoretical justification that the solution manifold for TRS/Regularization problems is spanned by a very low-dimensional subspace, independent of the problem dimension $n$.
2. Single Factorization: The method requires only one matrix factorization (for $A$ or a shifted $A$) regardless of the number of iterations, contrasting with multi-factorization methods.
3. Efficient Residual Computation: Derivation of formulas to compute the residual norm $\|r_k\|$ using only the small projected matrices and basis vectors, avoiding expensive full-vector residual calculations at every step.
4. Warm-Start Capability: The algorithm naturally supports "resolves" (solving for a smaller radius $\Delta$ with the same $A$ and $b$) by continuing the iteration from the previous stopping point, making it highly efficient for optimization sequences.
5. Generalization: Extension to indefinite Hessians (via diagonal shifting) and elliptical norms (via preconditioning).

4. Numerical Results

The authors tested TREK (and NREK) against:

• TRS: A multi-factorization method (GALAHAD library).
• GLTR: A standard factorization-free Krylov method.

Experimental Setup:

• 93 large-scale unconstrained test problems from CUTEst ($n \ge 1000$).
• Tested with decreasing trust-region radii ($\Delta = 10 \to 1 \to 0.1$).

Findings:

• Performance: TREK is competitive with or superior to TRS and GLTR in many cases.
• Factorization vs. Iteration: TREK excels when factorization is expensive but the solution requires a shift $\sigma$ near the smallest eigenvalue (where standard Krylov methods like GLTR converge slowly).
• Resolves: TREK is significantly faster for sequences of solves with decreasing radii because it reuses the subspace, whereas TRS often requires re-factorization.
• Iteration Count: The method typically converges in very few iterations (often $< 20$), confirming the low-rank hypothesis.
• Failure Cases: TREK failed in 6 cases (mostly due to iteration limits or singularity issues in specific edge cases), comparable to the failure rate of TRS.

5. Significance and Conclusion

This paper introduces a robust alternative to the standard "factorization vs. matrix-free" dichotomy in optimization subproblem solvers.

• Hybrid Efficiency: It combines the stability of factorization-based methods (via the initial solve) with the low iteration count of rational approximation methods.
• Scalability: By relying on a single factorization and low-dimensional projections, it is highly scalable for large $n$.
• Practical Impact: The implementation is available in the GALAHAD library, providing a new tool for large-scale nonlinear optimization solvers that rely on trust-region or cubic-regularization steps.

The work demonstrates that exploiting the specific structure of the solution manifold (via Extended Krylov Subspaces) yields a method that is often faster than multi-factorization approaches and more robust than standard Krylov methods, particularly in the "hard" cases where the shift parameter is small.