cjones6/cubic_reg

Implementation of Nesterov and Polyak's (2006) cubic regularization algorithm and Cartis et al's (2011) adaptive cubic regularization algorithm

cubic_reg

This code implements two algorithms:

1. Nesterov and Polyak's (2006) cubic regularization algorithm; and
2. Cartis et al's (2011) adaptive cubic regularization algorithm.

Cubic regularization solves unconstrained minimization problems by minimizing a cubic upper bound to the function at each iteration.

See the example.py file for an example of how to use them and the comments in cubic_reg.py for all of the possible input options. Briefly, you can load the file and numpy using

import numpy as np
import src.cubic_reg

specify your function, the gradient, Hessian, and initial point (the gradient and Hessian can be None)

f = lambda x: x[0] ** 2 * x[1] ** 2 + x[0] ** 2 + x[1] ** 2
grad = lambda x: np.asarray([2 * x[0] * x[1] ** 2 + 2 * x[0], 2 * x[0] ** 2 * x[1] + 2 * x[1]])
hess = lambda x: np.asarray([[2 * x[1] ** 2 + 2, 4 * x[0] * x[1]], [4 * x[0] * x[1], 2 * x[0] ** 2 + 2]])
x0 = np.array([1, 2]

and then use cubic regularization by running

cr = src.cubic_reg.CubicRegularization(x0, f=f, gradient=grad, hessian=hess, conv_tol=1e-4)
x_opt, intermediate_points, n_iter, flag = cr.cubic_reg()

To run adaptive cubic regularization instead, you can set

cr = src.cubic_reg.AdaptiveCubicReg(x0, f=f, gradient=grad, hessian=hess, hessian_update_method='broyden', conv_tol=1e-4)
x_opt, intermediate_points, n_iter, flag = cr.adaptive_cubic_reg()

There are many other options you can specify and parameters you can control.

References:

• Nesterov, Y., & Polyak, B. T. (2006). Cubic regularization of Newton method and its global performance. Mathematical Programming, 108(1), 177-205.
• Cartis, C., Gould, N. I., & Toint, P. L. (2011). Adaptive cubic regularisation methods for unconstrained optimization. Part I: motivation, convergence and numerical results. Mathematical Programming, 127(2), 245-295.
• Conn, A. R., Gould, N. I., & Toint, P. L. (2000). Trust region methods (Vol. 1). Siam.
• Gould, N. I., Lucidi, S., Roma, M., & Toint, P. L. (1999). Solving the trust-region subproblem using the Lanczos method. SIAM Journal on Optimization, 9(2), 504-525.

cubic_reg's Issues

How to solve the subproblem

What approach do you take to solving subproblems in the literature? Is this an approximate way of solving subproblems or an exact way of solving them? It the article is "Trust region methods (Vol. 1). Siam." ?

Possible incorrect handling of f_xold

In cubic_reg.py, f_xold (line 302) is calculated only once, then used as f_x for every calculation of rho going forward. Is this correct?

Possible incorrect handling of hard case

Hey again! Thanks for your quick response on my previous issue. I may have stumbled across another one. Consider the following:

N = 50

B = np.random.rand(N, N)
B += B.T
g = np.random.rand(N)
lambda_nplus = max(-scipy.linalg.eigh(B, eigvals_only=True, eigvals=(0, 0)), 0)
eigs, V = np.linalg.eigh(B)
# make sure indefinite
assert min(eigs) < 0 and max(eigs) > 0
v1 = V[:, 0]

# make gamma_1 = Ug_1 == 0
g = v1 - g * (v1.T @ v1)/(g.T @ v1)
assert abs(v1 @ g) < 1e-10
assert (V.T @ g)[0] < 1e-10
p = _AuxiliaryProblem(None, g, B, 10, lambda_nplus, 1e-4, 1e4)

print(p.solve())

This should trigger the hard-case (B indefinite, Ug_1 == 0) [1], but it does not reliably. Sometimes, it will trigger the hard case code, but then the output will be imaginary.

I'd love to submit a pull-request with a fix, but I'm actually not quite sure on how to fix this one!

[1]: https://people.maths.ox.ac.uk/cartis/papers/ARCpI.pdf -- Page 27, just below formula 6.5.