import some functions from stochopy to test the code
[1]:
# python lbfgsb works with numpy arrays
import numpy as np
from lbfgsb import minimize_lbfgsb
from scipy.optimize import minimize
import logging
from lbfgsb.benchmarks import (
ackley,
ackley_grad,
beale,
beale_grad,
griewank,
griewank_grad,
quartic,
quartic_grad,
rastrigin,
rastrigin_grad,
rosenbrock,
rosenbrock_grad,
sphere,
sphere_grad,
styblinski_tang,
styblinski_tang_grad,
)
from lbfgsb.utils import get_grad_projection_inf_norm
import matplotlib.pyplot as plt
[2]:
from typing import Callable
from lbfgsb.types import NDArrayFloat
def plot_2d(lb: NDArrayFloat, ub: NDArrayFloat, f: Callable) -> None:
x = np.linspace(lb, ub)
plt.contourf()
Define a logger for the display
[3]:
logging.info("Test")
main_logger = logging.getLogger(" -MAIN- ")
main_logger.setLevel(level=logging.INFO)
main_logger.info("test the main logger")
solver_logger = logging.getLogger("L-BFGS-B")
solver_logger.setLevel(level=logging.INFO)
solver_logger.info("test the solver logger")
INFO: -MAIN- :test the main logger
INFO:L-BFGS-B:test the solver logger
[4]:
iprint = -1
First example: min x^{\mathrm{T}}x such that x>=1. Optimal solution hits the boundary.
[5]:
def quad(x):
return 10 * x[0] ** 2 + x[1] ** 2
def grad_quad(x):
return np.array([20 * x[0], 2 * x[1]])
ftol = 1e-5
gtol = 1e-5
l = np.array([1.0, 1.0])
u = np.array([np.inf, np.inf])
bounds = np.array((l, u)).T
x0 = np.array([5.0, 5.0])
[6]:
bounds[0]
[6]:
array([ 1., inf])
New implementation#
[7]:
main_logger.info("====================== Quadratic example ======================")
opt_quad = minimize_lbfgsb(
x0=x0,
fun=quad,
jac=grad_quad,
bounds=bounds,
ftol=ftol,
gtol=gtol,
iprint=iprint,
logger=solver_logger,
)
main_logger.info("")
main_logger.info("")
xOpt = np.array([1, 1])
theoOpt = {"x": xOpt, "f": quad(xOpt), "df": grad_quad(xOpt)}
main_logger.info("Theoretical optimal value: ")
main_logger.info(theoOpt)
main_logger.info("Optimal value found: ")
main_logger.info(opt_quad)
INFO: -MAIN- :====================== Quadratic example ======================
INFO: -MAIN- :
INFO: -MAIN- :
INFO: -MAIN- :Theoretical optimal value:
INFO: -MAIN- :{'x': array([1, 1]), 'f': np.int64(11), 'df': array([20, 2])}
INFO: -MAIN- :Optimal value found:
INFO: -MAIN- : message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
success: True
status: 1
fun: 11.0
x: [ 1.000e+00 1.000e+00]
nit: 2
jac: [ 2.000e+01 2.000e+00]
nfev: 3
njev: 3
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
Scipy#
[8]:
minimize(
quad,
x0,
jac=grad_quad,
bounds=bounds,
method="l-bfgs-b",
options={"gtol": gtol, "ftol": ftol},
)
[8]:
message: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
success: True
status: 0
fun: 11.0
x: [ 1.000e+00 1.000e+00]
nit: 2
jac: [ 2.000e+01 2.000e+00]
nfev: 3
njev: 3
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
Second example: Ackley function#
[9]:
l = np.array([-2, -2, -1])
u = np.array([2, 2, 3])
bounds = np.array((l, u)).T
# x0 = np.array([0.12, 0.12])
x0 = np.array([-0.8, -1, 0.1])
[10]:
opt_akley = minimize_lbfgsb(
x0=x0,
fun=ackley,
jac=ackley_grad,
bounds=bounds,
maxcor=5,
ftol=ftol,
gtol=gtol,
iprint=iprint,
ftol_linesearch=1e-3,
gtol_linesearch=0.9,
xtol_linesearch=0.1,
logger=solver_logger,
)
main_logger.info("")
main_logger.info("")
theoOpt = {"x": np.array([1, 1, 2]), "f": 0, "df": ackley_grad(np.array([1, 1, 3]))}
main_logger.info("Theoretical optimal value: ")
main_logger.info(opt_akley)
main_logger.info("Optimal value found: ")
INFO: -MAIN- :
INFO: -MAIN- :
INFO: -MAIN- :Theoretical optimal value:
INFO: -MAIN- : message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FTOL
success: True
status: 0
fun: 3.4493521599753656
x: [-8.000e-01 -1.000e+00 1.000e-01]
nit: 0
jac: [-3.232e+00 -1.550e+00 -1.076e+00]
nfev: 21
njev: 21
hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
INFO: -MAIN- :Optimal value found:
[11]:
minimize(
x0=x0,
fun=ackley,
jac=ackley_grad,
bounds=bounds,
method="l-bfgs-b",
options=dict(
maxcor=5,
ftol=ftol,
gtol=gtol,
),
)
[11]:
message: ABNORMAL:
success: False
status: 2
fun: 3.4493521599753656
x: [-8.000e-01 -1.000e+00 1.000e-01]
nit: 0
jac: [-3.232e+00 -1.550e+00 -1.076e+00]
nfev: 21
njev: 21
hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
Second example : min Rosenbrock function#
[12]:
l = np.array([-2, -2])
u = np.array([2, 2])
bounds = np.array((l, u)).T
# x0 = np.array([0.12, 0.12])
x0 = np.array([-0.8, -1])
New implementation#
[13]:
opt_rosenbrock = minimize_lbfgsb(
x0=x0,
fun=rosenbrock,
jac=rosenbrock_grad,
bounds=bounds,
maxcor=5,
ftol=ftol,
gtol=gtol,
iprint=iprint,
ftol_linesearch=1e-3,
gtol_linesearch=0.9,
xtol_linesearch=0.1,
logger=solver_logger,
)
main_logger.info("")
main_logger.info("")
theoOpt = {"x": np.array([1, 1]), "f": 0, "df": rosenbrock_grad(np.array([1, 1]))}
main_logger.info("Theoretical optimal value: ")
main_logger.info(theoOpt)
main_logger.info("Optimal value found: ")
INFO: -MAIN- :
INFO: -MAIN- :
INFO: -MAIN- :Theoretical optimal value:
INFO: -MAIN- :{'x': array([1, 1]), 'f': 0, 'df': array([0., 0.])}
INFO: -MAIN- :Optimal value found:
[14]:
main_logger.info(opt_rosenbrock)
INFO: -MAIN- : message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FTOL
success: True
status: 0
fun: 3.8738696325664974e-07
x: [ 9.994e-01 9.988e-01]
nit: 18
jac: [-1.603e-03 1.794e-04]
nfev: 24
njev: 24
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
[15]:
rosenbrock_grad(opt_rosenbrock.x)
[15]:
array([-0.00160319, 0.00017937])
Scipy#
[16]:
opt_rosenbrock_sp = minimize(
rosenbrock,
x0,
jac=rosenbrock_grad,
bounds=bounds,
method="l-bfgs-b",
options={"gtol": gtol, "ftol": ftol, "iprint": iprint, "maxcor": 5},
)
C:\Users\ancollet\AppData\Local\Temp\ipykernel_32856\1020373109.py:1: DeprecationWarning: scipy.optimize: The `disp` and `iprint` options of the L-BFGS-B solver are deprecated and will be removed in SciPy 1.18.0.
opt_rosenbrock_sp = minimize(
[17]:
opt_rosenbrock_sp
[17]:
message: CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH
success: True
status: 0
fun: 3.873869631473485e-07
x: [ 9.994e-01 9.988e-01]
nit: 19
jac: [-1.603e-03 1.794e-04]
nfev: 24
njev: 24
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
[18]:
def beale2(x):
return (
(1.5 - x[0] + x[0] * x[1]) ** 2
+ (2.25 - x[0] + x[0] * x[1] ** 2) ** 2
+ (2.625 - x[0] + x[0] * x[1] ** 3) ** 2
)
def beale_grad2(x):
y1 = x[1]
y2 = y1 * y1
y3 = y2 * y1
f1 = 1.5 - x[0] + x[0] * y1
f2 = 2.25 - x[0] + x[0] * y2
f3 = 2.625 - x[0] + x[0] * y3
return np.array(
[
2 * (y1 - 1) * f1 + 2 * (y2 - 1) * f2 + 2 * (y3 - 1) * f3,
2 * x[0] * f1 + 4 * x[0] * y1 * f2 + 6 * x[0] * y2 * f3,
]
)
ftol = 1e-14
gtol = 1e-10
l = -4.5 * np.ones(2)
u = -l
bounds = np.array((l, u)).T
x0 = np.array([2.5, -1.3])
New implementation#
[19]:
opt_beale = minimize_lbfgsb(
x0=x0,
fun=beale,
jac=beale_grad,
bounds=bounds,
ftol=ftol,
gtol=gtol,
logger=solver_logger,
iprint=iprint,
)
main_logger.info("")
main_logger.info("")
theoOpt = {
"x": np.array([3, 0.5]),
"f": beale(np.array([3, 0.5])),
"df": beale_grad(np.array([3, 0.5])),
}
main_logger.info("Theoretical optimal value: ")
main_logger.info(theoOpt)
main_logger.info("Optimal value found: ")
INFO: -MAIN- :
INFO: -MAIN- :
INFO: -MAIN- :Theoretical optimal value:
INFO: -MAIN- :{'x': array([3. , 0.5]), 'f': np.float64(0.0), 'df': array([0., 0.])}
INFO: -MAIN- :Optimal value found:
[20]:
main_logger.info(opt_beale)
INFO: -MAIN- : message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FTOL
success: True
status: 0
fun: 1.4336936823116482e-21
x: [ 3.000e+00 5.000e-01]
nit: 15
jac: [-7.089e-11 1.797e-10]
nfev: 20
njev: 20
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
Scipy#
[21]:
opt_beale_sp = minimize(
beale,
x0,
jac=beale_grad,
bounds=bounds,
method="l-bfgs-b",
options={"gtol": gtol, "ftol": ftol, "iprint": iprint},
)
C:\Users\ancollet\AppData\Local\Temp\ipykernel_32856\1355717965.py:1: DeprecationWarning: scipy.optimize: The `disp` and `iprint` options of the L-BFGS-B solver are deprecated and will be removed in SciPy 1.18.0.
opt_beale_sp = minimize(
[22]:
opt_beale_sp
[22]:
message: CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH
success: True
status: 0
fun: 1.4337282180350256e-21
x: [ 3.000e+00 5.000e-01]
nit: 16
jac: [-7.089e-11 1.797e-10]
nfev: 20
njev: 20
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
Comparision of convergence
[ ]:
x0 = np.array([-50])
def func(x: NDArrayFloat) -> float:
return x + np.exp(-10 * x)
# result = minimize(func, x0=10, method='L-BFGS-B',
# options={'maxls': 5, 'disp': 1})
def jac(x: NDArrayFloat) -> NDArrayFloat:
return 1.0 - 10 * np.exp(-10 * x)
res = minimize_lbfgsb(
x0=x0,
fun=func,
jac=jac,
maxls=5,
is_check_factorizations=True,
logger=solver_logger,
iprint=100,
)
solver_logger.info(res)
# assert res.message == "RESTART_FROM_LNSRCH"
INFO:L-BFGS-B:At iterate 0 , f= 1.404e+217 , |proj g|= 1.404e+218
INFO:L-BFGS-B:
INFO:L-BFGS-B:ITERATION 1
INFO:L-BFGS-B:---------------- CAUCHY entered-------------------
INFO:L-BFGS-B:---------------- exit CAUCHY----------------------
INFO:L-BFGS-B:1 variables are free at GCP, iter = 1
INFO:L-BFGS-B:Iteration #0 (max: 50): ||x||=5.000e+01, f(x)=1.404e+217, ||jac(x)||=1.404e+218, cdt_arret=1.404e+218 (eps=1.000e-05)
INFO:L-BFGS-B: message: ABNORMAL_TERMINATION_IN_LNSRCH
success: False
status: 2
fun: 1.4035922178528375e+217
x: [-5.000e+01]
nit: 0
jac: [-1.404e+218]
nfev: 6
njev: 6
hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
[24]:
import numpy as np
from lbfgsb.types import NDArrayFloat # for type hints, numpy array of floats
from lbfgsb.benchmarks import rosenbrock, rosenbrock_grad
lb = np.array([-2, -2]) # lower bounds
ub = np.array([2, 2]) # upper bounds
bounds = np.array((l, u)).T # The number of variables to optimize is len(bounds)
x0 = np.array([-0.8, -1]) # The initial guess
[25]:
from lbfgsb import minimize_lbfgsb
res = minimize_lbfgsb(
x0=x0,
fun=rosenbrock,
jac=rosenbrock_grad,
bounds=bounds,
ftol=1e-5,
gtol=1e-5,
maxcor=10,
)
res
[25]:
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FTOL
success: True
status: 0
fun: 5.834035727637504e-07
x: [ 9.994e-01 9.989e-01]
nit: 16
jac: [-2.171e-02 1.030e-02]
nfev: 19
njev: 19
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
[26]:
res_3_fun = minimize_lbfgsb(
x0=x0,
fun=rosenbrock,
jac=rosenbrock_grad,
bounds=bounds,
ftol=1e-5,
gtol=1e-5,
maxfun=3,
)
res_3_fun
[26]:
message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
success: True
status: 1
fun: 0.8619211711864526
x: [ 6.370e-01 4.913e-01]
nit: 1
jac: [-2.250e+01 1.709e+01]
nfev: 3
njev: 3
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
[27]:
res_final = minimize_lbfgsb(
x0=res_3_fun.x,
fun=rosenbrock,
jac=rosenbrock_grad,
bounds=bounds,
ftol=1e-5,
gtol=1e-5,
)
res_final
[27]:
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FTOL
success: True
status: 0
fun: 1.349454252016592e-10
x: [ 1.000e+00 1.000e+00]
nit: 14
jac: [ 1.030e-04 -6.267e-05]
nfev: 21
njev: 21
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
[28]:
res_final = minimize_lbfgsb(
x0=res_3_fun.x,
fun=rosenbrock,
jac=rosenbrock_grad,
bounds=bounds,
ftol=1e-5,
gtol=1e-5,
checkpoint=res_3_fun,
)
res_final
[28]:
message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FTOL
success: True
status: 0
fun: 5.834035727637504e-07
x: [ 9.994e-01 9.989e-01]
nit: 16
jac: [-2.171e-02 1.030e-02]
nfev: 19
njev: 19
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
[29]:
import numpy as np
from scipy.optimize import OptimizeResult
def my_custom_callback(
xk: np.typing.NDArray[np.float64], state: OptimizeResult
) -> bool:
# if the objective function has been called 10 times or more => stop
if state.nfev >= 10:
return True
return False
res_callback = minimize_lbfgsb(
x0=x0,
fun=rosenbrock,
jac=rosenbrock_grad,
bounds=bounds,
ftol=1e-5,
gtol=1e-5,
callback=my_custom_callback,
)
res_callback
[29]:
message: STOP: USER CALLBACK
success: True
status: 2
fun: 0.08537354414890938
x: [ 7.354e-01 5.284e-01]
nit: 8
jac: [ 3.115e+00 -2.478e+00]
nfev: 10
njev: 10
hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
High dimentional rosenbrock
[30]:
x0
rosenbrock
[30]:
<function lbfgsb.benchmarks.rosenbrock(x: numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]) -> float>
[31]:
lb = np.array([-2, -2]) # lower bounds
ub = np.array([2, 2]) # upper bounds
bounds = np.array((lb, ub)).T # The number of variables to optimize is len(bounds)
x0 = np.array([-0.8, -1]) # The initial guess
[32]:
bounds.shape
[32]:
(2, 2)
[33]:
n = 1_000_000
x0_large = np.repeat(x0, n)
bounds_large = np.tile(bounds, n).reshape(-1, 2)
[34]:
import numdifftools as nd
def rosenbrock_large(x: NDArrayFloat) -> float:
return np.sum(rosenbrock(x.reshape(2, -1)))
def rosenbrock_large_grad(x: NDArrayFloat) -> NDArrayFloat:
return rosenbrock_grad(x.reshape(2, -1)).ravel()
# np.testing.assert_allclose(
# rosenbrock_large_grad(x0_large), nd.Gradient(rosenbrock_large, step=1e-5)(x0_large), atol=1e-5, rtol=1e-5
# )
rosenbrock_large(x0_large), rosenbrock_large_grad(x0_large)
[34]:
(np.float64(272199999.99999994),
array([-528.4, -528.4, -528.4, ..., -328. , -328. , -328. ],
shape=(2000000,)))
[35]:
res_large = minimize_lbfgsb(
x0=x0_large,
fun=rosenbrock_large,
jac=rosenbrock_large_grad,
bounds=bounds_large,
ftol=1e-5,
gtol=1e-5,
is_use_numba_jit=True,
)
res_large
[35]:
message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
success: True
status: 1
fun: 1.5703252811463446e-06
x: [ 1.000e+00 1.000e+00 ... 1.000e+00 1.000e+00]
nit: 19
jac: [ 4.375e-06 4.375e-06 ... -3.429e-06 -3.429e-06]
nfev: 22
njev: 22
hess_inv: <2000000x2000000 LbfgsInvHessProduct with dtype=float64>
[36]:
import cProfile
cmd = """
res_large = minimize_lbfgsb(
x0=x0_large,
fun=rosenbrock_large,
jac=rosenbrock_large_grad,
bounds=bounds_large,
ftol=1e-5,
gtol=1e-5,
is_use_numba_jit=False
)"""
cProfile.run(cmd, "bench.perf")
[37]:
res_large = minimize(
x0=x0_large,
fun=rosenbrock_large,
jac=rosenbrock_large_grad,
bounds=bounds_large,
options=dict(
ftol=1e-5,
gtol=1e-5,
),
)
res_large
[37]:
message: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
success: True
status: 0
fun: 1.570325215345258e-06
x: [ 1.000e+00 1.000e+00 ... 1.000e+00 1.000e+00]
nit: 19
jac: [ 4.375e-06 4.375e-06 ... -3.429e-06 -3.429e-06]
nfev: 22
njev: 22
hess_inv: <2000000x2000000 LbfgsInvHessProduct with dtype=float64>