Scalable adaptive cubic regularization methods

Adaptive cubic regularization (ARC) methods for unconstrained optimization compute steps from linear systems involving a shifted Hessian in the spirit of the Levenberg-Marquardt and trust-region methods. The standard approach consists in performing an iterative search for the shift akin to solving the secular equation in trust-region methods. Such search requires computing the Cholesky factorization of a tentative shifted Hessian at each iteration, which limits the size of problems that can be reasonably considered. We propose a scalable implementation of ARC named ARC\(_q\)K in which we solve a set of shifted systems concurrently by way of an appropriate modification of the Lanczos formulation of the conjugate gradient (CG) method. At each iteration of ARC\(_q\)K to solve a problem with \(n\) variables, a range of \(m \ll n\) shift parameters is selected. The computational overhead in CG beyond the Lanczos process is thirteen scalar operations to update five vectors of length \(m\), and two \(n\)-vector updates for each value of the shift. The CG variant only requires one Hessian-vector product and one dot product per iteration, independently of \(m\). Solves corresponding to inadequate shift parameters are interrupted early. All shifted systems are solved inexactly. Such modest cost makes our implementation scalable and appropriate for large-scale problems. We provide a new analysis of the inexact ARC method including its worst case evaluation complexity, global and asymptotic convergence. We describe our implementation and report numerical experience that confirms that our implementation of ARC\(_q\)K outperforms a classic Steihaug-Toint trust-region method, and the ARC method of the GALAHAD library. The latter solves the subproblem in nested Krylov subspaces by a Lanczos-based method, which requires the storage of a dense matrix that can be comparable to or larger than the two dense arrays required by our approach if the problem is large or requires many Lanczos iterations. Finally, we generalize our convergence results to inexact Hessians and nonlinear least-squares problems.

1. An iterate was generated where the objective evaluated to \(-\infty \).

Authors and Affiliations

1. GERAD and Département d’Informatique, Université de Sherbrooke, Sherbrooke, QC, Canada

   Jean-Pierre Dussault

2. GERAD and Mathematics and Industrial Engineering Department, Polytechnique Montréal, Montréal, QC, Canada

   Tangi Migot & Dominique Orban

Corresponding author

Correspondence to Dominique Orban.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Research partially supported by an NSERC Discovery Grant, Research supported by IVADO and the Canada First Research Excellence Fund / Apogée, Research partially supported by an NSERC Discovery Grant.

Appendix A. Example of an ill-conditioned problem

We examine the iterations of the ARC\(_q\)K method on problem sbrybnd from the CUTEst test set [20], which has unfavorable eigenvalues. The problem is a Broyden banded system of nonlinear equations with 5000 variables, solved in the least squares sense. The smallest eigenvalue is of order \(-10^{12}\) at the initial point. With the parameters of Sect. 4, ARC\(_q\)K and ARC both stop because they attain the time limit. Both solvers stop close to stationary point in the sense that \(\Vert \nabla f(x_k)\Vert \approx {\tilde{\epsilon }} \Vert \nabla f(x_0)\Vert \) with \({\tilde{\epsilon }} = 2.26 \cdot 10^{-5}\) after \(148\) iterations for ARC\(_q\)K, and \({\tilde{\epsilon }} = 2.45 \cdot 10^{-4}\) after \(299\) iterations for ARC.

Figure 5 reports the shift \(\lambda _k\) used by ARC\(_q\)K and the absolute value of the smallest eigenvalue \(\sigma _k\) of \(\nabla ^2 f(x_k)\) as a function of the iterations (left), and the norm of the search direction \(d_k\) (right). As expected, at iterations where \(\sigma _k\) is very negative, \(\lambda _k\) is very large, and \(\Vert d_k\Vert \) is small, which induces slow progresses. Figure 6 shows the corresponding quantities for ARC, and the conclusions are similar.

Problem sbrybnd illustrates the necessity to make provision for potentially large shifts on problems with unfavorable eigenvalue distribution. Both methods make slow progress in regions dominated by very negative eigenvalues. They both manage to solve the problem to an acceptable level, though not to the prescribed stopping tolerance. Of course, large shifts may be unnecessary on many problems. One can easily imagine a refined implementation in which large shifts are only employed when smaller shifts fail to provide a satisfactory search direction. Evidently, the issue would be mitigated by an appropriate problem scaling, but that is out of the scope of the present analysis.

Left: shifts \(\lambda _k\) used by ARC\(_q\)K and the absolute value of the smallest eigenvalue \(\sigma _k\) of \(\nabla ^2 f(x_k)\) as a function of the iteration. Right: norm of the direction \(d_k\) as a function of the iteration

Full size image

Left: shifts used by ARC estimated as the ratio \(\alpha _k / \Vert d_k\Vert \) where \(\alpha _k\) is ARC’s regularization parameter and \(d_k\) s the search direction. Right: norm of the direction \(d_k\) as a function of the iteration

Full size image

Appendix B Technical results

Denote by \(\sigma _{1,k}\le ...\le \sigma _{n,k}\) the eigenvalues of \(\nabla ^2 f(x_k)\), and consider an eigendecomposition of \(\nabla ^2 f(x_k) = U_k^T\Sigma _k U_k\) where \(\Sigma _k = diag (\sigma _{1,k}, \ldots , \sigma _{n,k})\), and \(U_k\) is an orthonormal matrix of eigenvectors.

As described in the preamble of Sect. 3, at each iteration k, the sequence of shifts is defined such that

where \(\lambda _0:= \psi ^{-e_L}\) and \(\psi > 1\).

Assume that there exists \({\bar{\imath }}_k\) the smallest nonnegative integer such that

Note that such an \({\bar{\imath }}_k\) always exists if the eigenvalues are bounded and m is sufficiently large. The motivation of this section is to prove that there exists a constant \(\beta \) such that for some \(i \ge {\bar{\imath }}_k\), (2.2c) is satisfied,

assuming that for all \(\lambda ^k_{i}\), there exists an oracle returning \(d(\lambda ^k_{i})\) satisfying (2.2a) with \(r^k_i\). First, the following lemma shows that \(\phi \) is monotonic.

Lemma B.1

Assume that \(\nabla f(x_k) \ne 0\). Let \(\Gamma _k > 0\) be a constant such that

Then, the sequence \(i \mapsto \phi (\lambda ^k_{i})\) is strictly increasing for \(i \ge {\bar{\imath }}_k\).

Conditions (B.1) and (B.2) impose a constraint on the precision used to solve the linear system, which is asymptotically less stringent than Assumption 2.6 used in the proofs.

Proof

The beginning of the proof follows the same line as [8, Section 7.3.1]. Using the eigendecomposition of the \(\nabla ^2 f(x_k)\), it follows that

where \(\gamma _{j,k} = [U_k \nabla f(x_k)]_j\) is the j-th component of \(U_k \nabla f(x_k)\). Additionally, using the triangle inequality and that \(\Vert a - b\Vert \ge \Vert a\Vert - \Vert b\Vert \), it holds for all \(i \ge {\bar{\imath }}_k\) that

We now prove that \(\Vert d_k(\lambda _{i+1}^k)\Vert _2^2 < \Vert d_k(\lambda ^k_i)\Vert _2^2\), which would conclude the proof. For \(i \ge {\bar{\imath }}_k\), using previous inequality, we get

Considering the first part of the right-hand side in previous inequality yields

The proof is complete after combining this relation with (B.2). \(\square \)

The following lemma gives an upper bound independent of i of

Lemma B.2

For all \(i \ge {\bar{\imath }}_k\),

Proof

The case \(i = {\bar{\imath }}_k\) is straightforward. Let \(i > {\bar{\imath }}_k\). Using \(\sigma _{1,k} + \psi ^{{\bar{\imath }}_k}\lambda _0 > 0\), then

On the other hand,

Combining the two previous inequalities into the definition of \({\overline{\kappa }}_{k}(i)\) yields

and the result follows as \(\psi ^{i - {\bar{\imath }}_k} \ge \psi \). \(\square \)

Proposition B.1

Let the assumptions of Lemma B.1 hold. Moreover, let \(\beta _k \ge 1\) be chosen such that

where \(\Gamma _k\) and \({\overline{\kappa }}\) are respectively the constants defined in Lemma B.1 and Lemma B.2. Then, there exists \(i \ge {\bar{\imath }}_k\) such that (2.2c) holds.

Proof

We can further assume that \(\phi (\lambda ^k_{{\bar{\imath }}_k}) < \frac{1}{\beta _k}\), as otherwise \(\phi (\lambda ^k_{{\bar{\imath }}_k}) \ge \frac{1}{\beta _k}\) and the result would follow by (B.3). The condition (2.2c) can be reformulated as

so the targeted interval is of length \(2\log (\beta _k)\). By Lemma B.1, the sequence \(i \mapsto \phi (\lambda ^k_{i})\) is strictly increasing, \(\log (\phi (\lambda ^k_{{\bar{\imath }}_k})) < \log (\beta _k)\) from (B.3), and we now prove that

which would conclude the proof.

Let \(i \ge {\bar{\imath }}_k\), it follows from (2.2a) that

Taking the norm on both sides, and using the triangle inequality yields

Additionally, \(\lambda ^k_i = \psi ^i \lambda _0\) gives

Therefore, using that \(\log \) is increasing and that \(x \mapsto \log (1+x)\) is subadditive,

which proves (B.5), and the result. \(\square \)

Proposition B.2

Let \(\{x^k\}\) be a sequence generated by Algorithm 3. Assume that the eigenvalues of \(\nabla ^2 f(x_k)\) are bounded for all k. Assume further that there exists a constant \(\phi ^{\max }\) such that

Then, there exists a constant \(\beta \ge 1\) such that (B.5) holds for all k. Furthermore, for all k, there exists a \(i \ge {\bar{\imath }}_k\) such that (2.2c) holds.

Proof

Using Lemma B.2 and that the eigenvalues of \(\nabla ^2 f(x_k)\) are bounded for all k, it is clear that \(\beta \) can be chosen greater than \(\phi ^{\max }\) and \(\sqrt{\psi {\overline{\kappa }}_k (1 + \Gamma _k)}\) for all k. \(\square \)

The condition (B.6) further imposes that the shifts \(\lambda ^k_i\) covers sufficiently small values.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions