LBFGS++

https://lbfgspp.statr.me/
Admin
A header-only C++ library implementing L-BFGS and L-BFGS-B algorithms for unconstrained and...
Tokens:49,342
Snippets:266
Trust Score:8.7
Update:5 months ago
Show doc for...
Tokens:
Raw
### L-BFGS Minimization Example

Source: https://lbfgspp.statr.me/doc/index.html

Demonstrates how to use the LBFGS++ library to minimize a function. It shows setting up parameters, creating a solver instance, providing an initial guess, and running the minimization.

```cpp
#include <Eigen/Core>
#include <iostream>
#include <LBFGS.h>
using Eigen::VectorXd;
using namespace LBFGSpp;
class Rosenbrock
{
private:
int n;
public:
Rosenbrock(int n_) : n(n_) {}
double operator()(const VectorXd& x, VectorXd& grad)
{
double fx = 0.0;
for(int i = 0; i < n; i += 2)
{
double t1 = 1.0 - x[i];
double t2 = 10 * (x[i + 1] - x[i] * x[i]);
grad[i + 1] = 20 * t2;
grad[i] = -2.0 * (x[i] * grad[i + 1] + t1);
fx += t1 * t1 + t2 * t2;
}
return fx;
}
};

int main()
{
const int n = 10;
// Set up parameters
LBFGSParam<double> param;
param.epsilon = 1e-6;
param.max_iterations = 100;
// Create solver and function object
LBFGSSolver<double> solver(param);
Rosenbrock fun(n);
// Initial guess
VectorXd x = VectorXd::Zero(n);
// x will be overwritten to be the best point found
double fx;
int niter = solver.minimize(fun, x, fx);
std::cout << niter << " iterations" << std::endl;
std::cout << "x = \n" << x.transpose() << std::endl;
std::cout << "f(x) = " << fx << std::endl;
return 0;
}

```

--------------------------------

### Box-Constrained Optimization with LBFGSBSolver

Source: https://lbfgspp.statr.me/doc

This example shows how to use the LBFGSBSolver for the Rosenbrock function with variables constrained between 2 and 4. It initializes the solver, defines bounds, and performs the minimization.

```cpp
#include <Eigen/Core>
#include <iostream>
#include <LBFGSB.h> // Note the different header file
using Eigen::VectorXd;
using namespace LBFGSpp;
class Rosenbrock
{
private:
int n;
public:
Rosenbrock(int n_) : n(n_) {}
double operator()(const VectorXd& x, VectorXd& grad)
{
double fx = 0.0;
for(int i = 0; i < n; i += 2)
{
double t1 = 1.0 - x[i];
double t2 = 10 * (x[i + 1] - x[i] * x[i]);
grad[i + 1] = 20 * t2;
grad[i] = -2.0 * (x[i] * grad[i + 1] + t1);
fx += t1 * t1 + t2 * t2;
}
return fx;
}
};
int main()
{
const int n = 10;
// Set up parameters
LBFGSBParam<double> param; // New parameter class
param.epsilon = 1e-6;
param.max_iterations = 100;
// Create solver and function object
LBFGSBSolver<double> solver(param); // New solver class
Rosenbrock fun(n);
// Bounds
VectorXd lb = VectorXd::Constant(n, 2.0);
VectorXd ub = VectorXd::Constant(n, 4.0);
// Initial guess
VectorXd x = VectorXd::Constant(n, 3.0);
// x will be overwritten to be the best point found
double fx;
int niter = solver.minimize(fun, x, fx, lb, ub);
std::cout << niter << " iterations" << std::endl;
std::cout << "x = \n" << x.transpose() << std::endl;
std::cout << "f(x) = " << fx << std::endl;
return 0;
}
```

--------------------------------

### L-BFGS Minimization Example

Source: https://lbfgspp.statr.me/

Sets up LBFGS++ parameters, creates a solver instance, provides an initial guess, and runs the minimization for the Rosenbrock function. Requires Eigen and LBFGS++ headers.

```C++
int main()
{
    const int n = 10;
    // Set up parameters
    LBFGSParam<double> param;
    param.epsilon = 1e-6;
    param.max_iterations = 100;

    // Create solver and function object
    LBFGSSolver<double> solver(param);
    Rosenbrock fun(n);

    // Initial guess
    VectorXd x = VectorXd::Zero(n);
    // x will be overwritten to be the best point found
    double fx;
    int niter = solver.minimize(fun, x, fx);

    std::cout << niter << " iterations" << std::endl;
    std::cout << "x = \n" << x.transpose() << std::endl;
    std::cout << "f(x) = " << fx << std::endl;

    return 0;
}

```

--------------------------------

### Compilation and Execution Command

Source: https://lbfgspp.statr.me/

Example command to compile the C++ code using g++, specifying include paths for Eigen and LBFGS++. Assumes the source file is named example.cpp.

```Shell
$ g++ -I/path/to/eigen -I/path/to/lbfgspp/include -O2 example.cpp
$ ./a.out
23 iterations
x =
1 1 1 1 1 1 1 1 1 1
f(x) = 1.87948e-19

```

--------------------------------

### BFGS Matrix Update Example

Source: https://lbfgspp.statr.me/doc/BFGSMat_8h_source.html

Illustrates a typical update step within the BFGS matrix calculation, involving vector operations and matrix-vector products. This snippet is part of a larger function likely responsible for updating the BFGS approximation.

```cpp
Vector MWQtv;
apply_Mv(WQtv, MWQtv);
MWQtv.tail(m_ncorr) *= m_theta;
res.noalias() = -WP * MWQtv;
return true;
```

--------------------------------

### Initialize Subspace Minimization Variables

Source: https://lbfgspp.statr.me/doc/SubspaceMin_8h_source.html

Initializes variables for subspace minimization, including bounds, free variable sets, and initial solution vectors. This setup is crucial before entering the iterative refinement phase.

```cpp
Matrix WF = bfgs.Wb(fv_set);
Vector vecc(nfree);
bfgs.compute_FtBAb(WF, fv_set, newact_set, Wd, drt, vecc);
Vector vecl(nfree), vecu(nfree);
for (int i = 0; i < nfree; i++)
{
    const int coord = fv_set[i];
    vecl[i] = lb[coord] - x0[coord];
    vecu[i] = ub[coord] - x0[coord];
    vecc[i] += g[coord];
}
Vector vecy(nfree);
bfgs.solve_PtBP(WF, -vecc, vecy);
```

--------------------------------

### Handle Infinite Function or Gradient Values

Source: https://lbfgspp.statr.me/doc/LineSearchMoreThuente_8h_source.html

If the function or gradient values are infinite, return the midpoint between the current interval start and a given point. This is a safeguard for numerical stability.

```C++
if (!std::isfinite(ft) || !std::isfinite(gt))
    return (al + at) / Scalar(2);
```

--------------------------------

### Forming the B Matrix in L-BFGS++

Source: https://lbfgspp.statr.me/doc/BFGSMat_8h_source.html

Constructs the B matrix used in L-BFGS, starting with an initial approximation and updating it based on historical data (Y and S matrices). Requires the `Eigen` library.

```C++
inline Matrix get_Bmat() const
{
// Initial approximation theta * I
const int n = m_s.rows();
Matrix B = m_theta * Matrix::Identity(n, n);
if (m_ncorr < 1)
return B;

// Construct W matrix, W = [Y, theta * S]
// Y = [y0, y1, ..., yc]
// S = [s0, s1, ..., sc]
// We first set W = [Y, S], since later we still need Y and S matrices
// After computing Minv, we rescale the S part in W
Matrix W(n, 2 * m_ncorr);
// r = m_ptr - 1 points to the most recent element,
// (r + 1) % m_ncorr points to the oldest element
int j = m_ptr % m_ncorr;
for (int i = 0; i < m_ncorr; i++)
{
W.col(i).noalias() = m_y.col(j);
W.col(m_ncorr + i).noalias() = m_s.col(j);
j = (j + 1) % m_m;
}
// Now Y = W[:, :c], S = W[:, c:]

// Construct Minv matrix, Minv = [-D L' ]
// [ L theta * S'S]

// D = diag(y0's0, ..., yc'sc)
Matrix Minv(2 * m_ncorr, 2 * m_ncorr);
Minv.topLeftCorner(m_ncorr, m_ncorr).setZero();
Vector ys = W.leftCols(m_ncorr).cwiseProduct(W.rightCols(m_ncorr)).colwise().sum().transpose();
Minv.diagonal().head(m_ncorr).noalias() = -ys;
// L = [ 0 ]
// [ s[1]'y[0] 0 ]
// [ s[2]'y[0] s[2]'y[1] ]
// ...
// [s[c-1]'y[0] ... ... ... ... ... s[c-1]'y[c-2] 0]
Minv.bottomLeftCorner(m_ncorr, m_ncorr).setZero();
for (int i = 0; i < m_ncorr - 1; i++)
{
// Number of terms for this column
const int nterm = m_ncorr - i - 1;
// S[:, -nterm:]'Y[:, j]
Minv.col(i).tail(nterm).noalias() = W.rightCols(nterm).transpose() * W.col(i);
}
// The symmetric block
Minv.topRightCorner(m_ncorr, m_ncorr).noalias() = Minv.bottomLeftCorner(m_ncorr, m_ncorr).transpose();
// theta * S'S
Minv.bottomRightCorner(m_ncorr, m_ncorr).noalias() = m_theta * W.rightCols(m_ncorr).transpose() * W.rightCols(m_ncorr);

// Set the true W matrix
W.rightCols(m_ncorr).array() *= m_theta;

// Compute B = theta * I - W * M * W'
Eigen::PartialPivLU<Matrix> M_solver(Minv);
B.noalias() -= W * M_solver.solve(W.transpose());
return B;
}
```

--------------------------------

### Line Search by Moré and Thuente

Source: https://lbfgspp.statr.me/doc/classLBFGSpp_1_1LineSearchMoreThuente.html

This static method performs a line search to find an optimal step size satisfying strong Wolfe conditions. It requires a function object, solver parameters, current point, direction, and initial step, and outputs the calculated step, function value, gradient, and new point.

```cpp
#include <LineSearchMoreThuente.h>

template<typename Scalar>
template<typename Foo, typename SolverParam>
static void LBFGSpp::LineSearchMoreThuente< Scalar >::LineSearch (
	Foo &_f_,
	const SolverParam &_param_,
	const Vector &_xp_,
	const Vector &_drt_,
	const Scalar &_step_max_,
	Scalar &_step_,
	Scalar &_fx_,
	Vector &_grad_,
	Scalar &_dg_,
	Vector &_x_
)
```

--------------------------------

### LBFGSBSolver Constructor

Source: https://lbfgspp.statr.me/doc/LBFGSB_8h_source.html

Initializes the L-BFGS-B solver with specified parameters.

```cpp
LBFGSBSolver(const LBFGSBParam< Scalar > &param)
```

--------------------------------

### Get b-th Row of W Matrix (L-BFGS-B)

Source: https://lbfgspp.statr.me/doc/BFGSMat_8h_source.html

Extracts the b-th row of the W matrix, which is composed of Y and theta * S columns. Returns the row as a column vector, preserving the ordering of Y and S.

```cpp
inline Vector Wb(int b) const
{
Vector res(2 * m_ncorr);
for (int j = 0; j < m_ncorr; j++)
{
res[j] = m_y(b, j);
res[m_ncorr + j] = m_s(b, j);
}
res.tail(m_ncorr) *= m_theta;
return res;
}
```

--------------------------------

### Bazel Integration for LBFGS++

Source: https://lbfgspp.statr.me/doc/index.html

This snippet shows how to integrate LBFGS++ into a project using the Bazel build system. It includes adding the lbfgspp module dependency and referencing it in a BUILD.bazel file.

```Bazel
bazel_dep(name = "lbfgspp", version = "4.0.0")
git_override(
module_name = "lbfgspp",
commit = "c524a407fb85b74807f53de5a3ca2ddbcc164e54",
remote = "https://github.com/yixuan/LBFGSpp.git",
)
# ...
cc_library(
name = "my_lib",
srcs = ["my_lib.cc"],
hdrs = ["my_lib.h"],
deps = ["@lbfgspp"],
)
```

--------------------------------

### Run L-BFGS Minimization

Source: https://lbfgspp.statr.me/doc

Sets up parameters, creates an L-BFGS solver, provides an initial guess, and runs the minimization function for the Rosenbrock problem. The output shows the number of iterations and the final solution.

```cpp
int main()
{
const int n = 10;
// Set up parameters
LBFGSParam<double> param;
param.epsilon = 1e-6;
param.max_iterations = 100;
// Create solver and function object
LBFGSSolver<double> solver(param);
Rosenbrock fun(n);
// Initial guess
VectorXd x = VectorXd::Zero(n);
// x will be overwritten to be the best point found
double fx;
int niter = solver.minimize(fun, x, fx);
std::cout << niter << " iterations" << std::endl;
std::cout << "x = \n" << x.transpose() << std::endl;
std::cout << "f(x) = " << fx << std::endl;
return 0;
}

```

--------------------------------

### L-BFGS Minimization with Bracketing Line Search

Source: https://lbfgspp.statr.me/

Demonstrates using LBFGS++ with a different line search algorithm (Bracketing) by specifying it as a template parameter to LBFGSSolver. Requires Eigen and LBFGS++ headers.

```C++
int main()
{
    const int n = 10;
    // Set up parameters
    LBFGSParam<double> param;
    param.epsilon = 1e-6;
    param.max_iterations = 100;

    // Create solver and function object
    LBFGSSolver<double, LineSearchBracketing> solver(param);
    Rosenbrock fun(n);

    // Initial guess
    VectorXd x = VectorXd::Zero(n);
    // x will be overwritten to be the best point found
    double fx;
    int niter = solver.minimize(fun, x, fx);

    std::cout << niter << " iterations" << std::endl;
    std::cout << "x = \n" << x.transpose() << std::endl;
    std::cout << "f(x) = " << fx << std::endl;

    return 0;
}

```

--------------------------------

### Forming the H Matrix in L-BFGS++

Source: https://lbfgspp.statr.me/doc/BFGSMat_8h_source.html

Constructs the H matrix (inverse Hessian approximation) used in L-BFGS, starting with an initial approximation and updating it based on historical data (Y and S matrices). Requires the `Eigen` library.

```C++
inline Matrix get_Hmat() const
{
// Initial approximation 1/theta * I
const int n = m_s.rows();
Matrix H = (Scalar(1) / m_theta) * Matrix::Identity(n, n);
if (m_ncorr < 1)
return H;

// Construct W matrix, W = [1/theta * Y, S]
// Y = [y0, y1, ..., yc]
// S = [s0, s1, ..., sc]
// We first set W = [Y, S], since later we still need Y and S matrices
// After computing M, we rescale the Y part in W
Matrix W(n, 2 * m_ncorr);
// p = m_ptr - 1 points to the most recent element,
// (p + 1) % m_ncorr points to the oldest element
int j = m_ptr % m_ncorr;
for (int i = 0; i < m_ncorr; i++)
{
W.col(i).noalias() = m_y.col(j);
W.col(m_ncorr + i).noalias() = m_s.col(j);
j = (j + 1) % m_m;
}
// Now Y = W[:, :c], S = W[:, c:]

// Construct M matrix, M = [ 0 -inv(R) ]
// [ -inv(R)' inv(R)'(D + 1/theta * Y'Y)inv(R) ]
// D = diag(y0's0, ..., yc'sc)
Matrix M(2 * m_ncorr, 2 * m_ncorr);
// First use M[:c, :c] to store R
// R = [s[0]'y[0] s[0]'y[1] ... s[0]'y[c-1] ]
// [ 0 s[1]'y[1] ... s[1]'y[c-1] ]
// ...
// [ 0 0 ... s[c-1]'y[c-1] ]
for (int i = 0; i < m_ncorr; i++)
{
M.col(i).head(i + 1).noalias() = W.middleCols(m_ncorr, i + 1).transpose() * W.col(i);
}
// Compute inv(R)
Matrix Rinv = M.topLeftCorner(m_ncorr, m_ncorr).template triangularView<Eigen::Upper>().solve(Matrix::Identity(m_ncorr, m_ncorr));
// Zero out the top left block
M.topLeftCorner(m_ncorr, m_ncorr).setZero();
// Set the top right block
M.topRightCorner(m_ncorr, m_ncorr).noalias() = -Rinv;
// The symmetric block
M.bottomLeftCorner(m_ncorr, m_ncorr).noalias() = -Rinv.transpose();
// 1/theta * Y'Y
Matrix block = (Scalar(1) / m_theta) * W.leftCols(m_ncorr).transpose() * W.leftCols(m_ncorr);
// D + 1/theta * Y'Y
Vector ys = W.leftCols(m_ncorr).cwiseProduct(W.rightCols(m_ncorr)).colwise().sum().transpose();
block.diagonal().array() += ys.array();
// The bottom right block
M.bottomRightCorner(m_ncorr, m_ncorr).noalias() = Rinv.transpose() * block * Rinv;

// Set the true W matrix
W.leftCols(m_ncorr).array() *= (Scalar(1) / m_theta);

// Compute H = 1/theta * I + W * M * W'
H.noalias() += W * M * W.transpose();
return H;
}
```

--------------------------------

### Minimize Function with Bound Constraints

Source: https://lbfgspp.statr.me/doc/LBFGSB_8h_source.html

Performs minimization of a function 'f' starting from initial guess 'x', subject to lower ('lb') and upper ('ub') bound constraints. It returns an integer status code.

```C++
template <typename Foo>
inline int minimize(Foo& f, Vector& x, Scalar& fx, const Vector& lb, const Vector& ub)
{
 using std::abs;

 // Dimension of the vector
 const int n = x.size();
 if (lb.size() != n || ub.size() != n)
 throw std::invalid_argument("'lb' and 'ub' must have the same size as 'x'");

 // Check whether the initial vector is within the bounds
 // If not, project to the feasible set
 force_bounds(x, lb, ub);

 // Initialization
 reset(n);

 // The length of lag for objective function value to test convergence
 const int fpast = m_param.past;

 // Evaluate function and compute gradient
 fx = f(x, m_grad);
 m_projgnorm = proj_grad_norm(x, m_grad, lb, ub);
 if (fpast > 0)
 m_fx[0] = fx;

 // std::cout << "x0 = " << x.transpose() << std::endl;
 // std::cout << "f(x0) = " << fx << ", ||proj_grad|| = " << m_projgnorm << std::endl << std::endl;

 // Early exit if the initial x is already a minimizer
 if (m_projgnorm <= m_param.epsilon || m_projgnorm <= m_param.epsilon_rel * x.norm())
 {
 return 1;
 }

 // Compute generalized Cauchy point
 Vector xcp(n), vecc;
 IndexSet newact_set, fv_set;
 Cauchy<Scalar>::get_cauchy_point(m_bfgs, x, m_grad, lb, ub, xcp, vecc, newact_set, fv_set);

 /* Vector gcp(n);
 Scalar fcp = f(xcp, gcp);
 Scalar projgcpnorm = proj_grad_norm(xcp, gcp, lb, ub);
 std::cout << "xcp = " << xcp.transpose() << std::endl;
 std::cout << "f(xcp) = " << fcp << ", ||proj_grad|| = " << projgcpnorm << std::endl << std::endl; */

 // Initial direction
 m_drt.noalias() = xcp - x;
 m_drt.normalize();
 // Tolerance for s'y >= eps * (y'y)

```

--------------------------------

### Step Size Selection with Gradient Magnitude Comparison

Source: https://lbfgspp.statr.me/doc/LineSearchMoreThuente_8h_source.html

Selects a step size between cubic and quadratic interpolations based on gradient magnitudes. This case applies when function values decrease but gradient signs are consistent.

```C++
const Scalar deltal = Scalar(1.1), deltau = Scalar(0.66);
if (abs(gt) < abs(gl))
{
    const Scalar res = (ac_exists &&
                        (ac - at) * (at - al) > Scalar(0) &&
                        abs(ac - at) < abs(as - at)) ?
                       ac :
                       as;
    return (at > al) ?
           (std::min)(at + deltau * (au - at), res) :
           (std::max)(at + deltau * (au - at), res);
}
```

--------------------------------

### Initialize Iteration Variables for Subspace Minimization

Source: https://lbfgspp.statr.me/doc/SubspaceMin_8h_source.html

Sets up data structures and variables for the iterative subspace minimization process, including fallback solutions and dual variables.

```cpp
Vector yfallback = vecy;
Vector lambda = Vector::Zero(nfree), mu = Vector::Zero(nfree);
IndexSet L_set, U_set, P_set, yL_set, yU_set, yP_set;
L_set.reserve(nfree / 3);
yL_set.reserve(nfree / 3);
U_set.reserve(nfree / 3);
yU_set.reserve(nfree / 3);
P_set.reserve(nfree);
yP_set.reserve(nfree);
```

--------------------------------

### LBFGSBSolver Constructor

Source: https://lbfgspp.statr.me/doc/classLBFGSpp_1_1LBFGSBSolver.html

Constructor for the L-BFGS-B solver. Initializes the solver with specified parameters.

```APIDOC
## LBFGSBSolver(const LBFGSBParam< Scalar > &param)

### Description
Constructor for the L-BFGS-B solver.

### Parameters
* **param** (const LBFGSBParam< Scalar > &) - An object of LBFGSParam to store parameters for the algorithm.
```

--------------------------------

### Line Search Strong Wolfe Conditions Explanation

Source: https://lbfgspp.statr.me/doc/LineSearchMoreThuente_8h_source.html

Explains the mathematical formulation and practical implementation of the strong Wolfe conditions for line search, including interval management and exit criteria.

```C++
template <typename Foo, typename SolverParam>
static void LineSearch(Foo& f, const SolverParam& param,
                       const Vector& xp, const Vector& drt, const Scalar& step_max,
                       Scalar& step, Scalar& fx, Vector& grad, Scalar& dg, Vector& x)
{
    using std::abs;
    // The strong Wolfe conditions are:
    // f(x + alpha * d) <= f(x) + alpha * mu * g^T d,
    // |g(x + alpha * d)^T d| <= eta * |g^T d|,
    // where mu and eta are two constants, and we typically require 0 < mu < eta < 1.
    // 
    // Let phi(alpha) = f(x + alpha * d), and then the conditions reduce to
    // phi(alpha) <= phi(0) + mu * phi'(0) * alpha,
    // |phi'(alpha)| <= eta * |phi'(0)|.
    // 
    // [1] defines a set
    // T(mu) = {alpha > 0: phi(alpha) <= phi(0) + mu * phi'(0) * alpha,
    // |phi'(alpha)| <= mu * |phi'(0)| }.
    // Clearly, if we have found some alpha that belongs to T(mu), and if mu < eta,
    // then this alpha must satisfy the strong Wolfe conditions.
    // 
    // In a practical implementation, we also impose bounds on alpha:
    // alpha in Ia = [alpha_min, alpha_max]
    // The lower bound alpha_min is used to terminate the iteration.
    // The upper bound is typically used for constrained problems (e.g., L-BFGS-B).
    // 
    // The overall framework of the algorithm is to generate a sequence of iterates alpha_k
    // and a sequence of intervals I_k such that:
    // 1. alpha_k is in intersect(I_k, Ia).
    // 2. I_k eventually is a subset of Ia.
    // 3. The length of I_k is shrinking.
    // In each iteration, we test whether alpha_k satisfies the strong Wolfe conditions,
    // and will exit the line search when it meets. Other possible exiting conditions are:
    // 1. alpha reaches alpha_max.
    // 2. alpha reaches alpha_min.
    // 3. The maximum number of iterations is attained.
    // 
    // To achieve this goal, we need to make sure that we can generate an interval I such that
    // intersect(I, T(mu)) is not empty, and that for the updated interval I+, intersect(I+, T(mu))
    // is also not empty. Theorem 2.1 of [1] gives conditions for this, which first defines
    // an auxiliary function:
    // psi(alpha) = phi(alpha) - phi(0) - mu * phi'(0) * alpha.
    // Theorem 2.1 of [1] states that for an interval I = [alpha_l, alpha_u]
    // (alpha_u can be smaller than alpha_l), if
    // psi(alpha_l) <= psi(alpha_u), (*)
    // psi(alpha_l) <= 0, (*)
    // psi'(alpha_l) * (alpha_u - alpha_l) < 0, (*)
}
```

--------------------------------

### Line Search Step Selection and Interval Update Logic

Source: https://lbfgspp.statr.me/doc/LineSearchMoreThuente_8h_source.html

This snippet details the core logic for selecting a new step size and updating the bracketing interval within the line search. It handles different cases (Case 1, 2, 3) based on function and gradient values and applies safeguarding rules.

```C++
const Scalar ft = f_is_psi ? psit : fx;
const Scalar gt = f_is_psi ? dpsit : dg;

// Check and update the status of use_step_min_safeguard
// We impose a safeguarding rule that guarantees testing
// step_min if psi(alpha_k) > 0 or psi'(alpha_k) >= 0
// holds from the beginning
if (use_step_min_safeguard && (psit <= Scalar(0) && dpsit < Scalar(0)))
{
  use_step_min_safeguard = false;
}

// Update new step
Scalar new_step;
const bool in_case_2 = (psit <= psiI_lo) && (dpsit * (I_lo - step) > Scalar(0));
if (in_case_2)
{
  // For Case 2, we apply the safeguarding rule
  // newat = min(at + delta * (at - al), amax), delta in [1.1, 4]
  new_step = (std::min)(step_max, step + delta_max * (step - I_lo));

  // We can also consider the following scheme:
  // First let step_selection() decide a value, and then project to the range above
  //
  // new_step = step_selection(I_lo, I_hi, step, fI_lo, fI_hi, ft, gI_lo, gI_hi, gt);
  // const Scalar delta2 = Scalar(4)
  // const Scalar t1 = step + delta * (step - I_lo);
  // const Scalar t2 = step + delta2 * (step - I_lo);
  // const Scalar tl = std::min(t1, t2), tu = std::max(t1, t2);
  // new_step = std::min(tu, std::max(tl, new_step));
  // new_step = std::min(step_max, new_step);
}
else
{
  // For Case 1 and Case 3, use information of f and g to select new step
  new_step = step_selection(I_lo, I_hi, step, fI_lo, fI_hi, ft, gI_lo, gI_hi, gt);
  // Force new step in [step_min, step_max]
  new_step = (std::max)(new_step, step_min);
  new_step = (std::min)(new_step, step_max);

  // Apply safeguarding rule related to step_min when necessary:
  // step+ in [alpha_min, max{delta_min * step, alpha_min}]
  //
  // If use_step_min_safeguard = true, then new_step cannot be obtained
  // from Case 2, since in Case 2 we have
  // psi(alpha_k) <= 0 and psi'(alpha_k) < 0
  if (use_step_min_safeguard)
  {
    const Scalar lower = step_min;
    const Scalar upper = (std::max)(step_min, delta_min * step);
    new_step = (std::max)(new_step, lower);
    new_step = (std::min)(new_step, upper);
  }
}

// Update bracketing interval
if (psit > psiI_lo)
{
  // Case 1: psi(step) > psi(I_lo)
  I_hi = step;
  fI_hi = ft;
  gI_hi = gt;

  // std::cout << "** Case 1" << std::endl;
}
else if (in_case_2)
{
  // Case 2: psi(step) <= psi(I_lo), psi'(step)(I_lo - step) > 0
  I_lo = step;
  fI_lo = ft;
  gI_lo = gt;
  psiI_lo = psit;
  // Move x and grad to x_lo and grad_lo, respectively
  x_lo.swap(x);
  grad_lo.swap(grad);
  fx_lo = fx;
  dg_lo = dg;

  // std::cout << "** Case 2" << std::endl;
}
else
{
  // Case 3: psi(step) <= psi(I_lo), psi'(step)(I_lo - step) <= 0
  I_hi = I_lo;
  fI_hi = fI_lo;
  gI_hi = gI_lo;

  I_lo = step;
  fI_lo = ft;
  gI_lo = gt;
  psiI_lo = psit;
  // Move x and grad to x_lo and grad_lo, respectively
  x_lo.swap(x);
  grad_lo.swap(grad);
  fx_lo = fx;
  dg_lo = dg;

  // std::cout << "** Case 3" << std::endl;
}

// Check and update the status of bracketed
// bracketed is true if we have entered Case 1 or Case 3,
// and I is contained in [step_min, step_max]
if ((!bracketed) && (!in_case_2))
{
  const Scalar I_left = (std::min)(I_lo, I_hi);
  const Scalar I_right = (std::max)(I_lo, I_hi);
  bracketed = (I_left >= step_min && I_right <= step_max);
}

// If bracketed, enforce sufficient interval shrink; if not shrinking enough, use bisection
if (bracketed)
{
  I_width_prev = I_width;
  I_width = abs(I_hi - I_lo);
  // Test interval shrinkage
  if (I_width_prev < Inf && I_width > shrink * I_width_prev)
  {
    I_shrink_fail_count += 1;
  }
  else
  {
    I_shrink_fail_count = 0;
  }
  // If interval fails to shrink enough twice, select new_step using bisection
  if (I_shrink_fail_count >= 2)
  {
    new_step = (I_lo + I_hi) / Scalar(2);
    I_shrink_fail_count = 0;
  }
}

// Set the new_step
step = new_step;

// std::cout << "[I+_lo, I+_hi] = [" << I_lo << ", " << I_hi << "], step+ = " << step << std::endl << std::endl;

```

--------------------------------

### LBFGSBSolver Constructor

Source: https://lbfgspp.statr.me/doc/LBFGSB_8h_source.html

Constructs an LBFGSBSolver instance with specified parameters. The parameters are checked for validity upon construction.

```C++
LBFGSBSolver(const LBFGSBParam<Scalar>& param) : 
 m_param(param)
{
 m_param.check_param();
}
```

--------------------------------

### LBFGSSolver Constructor

Source: https://lbfgspp.statr.me/doc/classLBFGSpp_1_1LBFGSSolver.html

Constructor for the L-BFGS solver. It initializes the solver with specified parameters.

```APIDOC
## LBFGSSolver(const LBFGSParam< Scalar > &param)

### Description
Constructor for the L-BFGS solver.

### Parameters
* **param** (LBFGSParam< Scalar >) - Required - An object of LBFGSParam to store parameters for the algorithm.
```

--------------------------------

### LBFGSBParam Constructor

Source: https://lbfgspp.statr.me/doc/Param_8h_source.html

Initializes LBFGSBParam with default values for L-BFGS-B algorithm parameters. Includes settings for history size, convergence criteria, and sub-minimization.

```cpp
LBFGSBParam()
{
    // clang-format off
    m = 6;
    epsilon = Scalar(1e-5);
    epsilon_rel = Scalar(1e-5);
    past = 1;
    delta = Scalar(1e-10);
    max_iterations = 0;
    max_submin = 10;
    max_linesearch = 20;
    min_step = Scalar(1e-20);
    max_step = Scalar(1e+20);
    ftol = Scalar(1e-4);
    wolfe = Scalar(0.9);
    // clang-format on
}
```

--------------------------------

### Box-Constrained Optimization with LBFGSBSolver

Source: https://lbfgspp.statr.me/

Use LBFGSBSolver for optimization problems where parameters have simple bounds. Set lower and upper bounds using VectorXd. Infinite bounds can be specified using std::numeric_limits<double>::infinity().

```cpp
#include <Eigen/Core>
#include <iostream>
#include <LBFGSB.h>  // Note the different header file

using Eigen::VectorXd;
using namespace LBFGSpp;

class Rosenbrock
{
private:
    int n;
public:
    Rosenbrock(int n_) : n(n_) {}
    double operator()(const VectorXd& x, VectorXd& grad)
    {
        double fx = 0.0;
        for(int i = 0; i < n; i += 2)
        {
            double t1 = 1.0 - x[i];
            double t2 = 10 * (x[i + 1] - x[i] * x[i]);
            grad[i + 1] = 20 * t2;
            grad[i]     = -2.0 * (x[i] * grad[i + 1] + t1);
            fx += t1 * t1 + t2 * t2;
        }
        return fx;
    }
};

int main()
{
    const int n = 10;
    // Set up parameters
    LBFGSBParam<double> param;  // New parameter class
    param.epsilon = 1e-6;
    param.max_iterations = 100;

    // Create solver and function object
    LBFGSBSolver<double> solver(param);  // New solver class
    Rosenbrock fun(n);

    // Bounds
    VectorXd lb = VectorXd::Constant(n, 2.0);
    VectorXd ub = VectorXd::Constant(n, 4.0);

    // Initial guess
    VectorXd x = VectorXd::Constant(n, 3.0);

    // x will be overwritten to be the best point found
    double fx;
    int niter = solver.minimize(fun, x, fx, lb, ub);

    std::cout << niter << " iterations" << std::endl;
    std::cout << "x = \n" << x.transpose() << std::endl;
    std::cout << "f(x) = " << fx << std::endl;

    return 0;
}

```

--------------------------------

### Cauchy Step Calculation Logic

Source: https://lbfgspp.statr.me/doc/Cauchy_8h_source.html

This snippet outlines the core logic for calculating the Cauchy step. It handles cases where all coordinates are at the boundary and initializes variables for interval-based calculations. It also computes initial values for f', f'', and the theoretical minimum step size.

```C++
const int nord = ord.size();
const int nfree = fv_set.size();
if ((nfree < 1) && (nord < 1))
{
/* std::cout << "** All coordinates at boundary **\n";
std::cout << "\n========================= Leaving GCP search =========================\n\n"; */
return;
}

// First interval: [il=0, iu=brk[ord[0]]]
// In case ord is empty, we take iu=Inf

// p = W'd, c = 0
Vector vecp;
bfgs.apply_Wtv(vecd, vecp);
// f' = -d'd
Scalar fp = -vecd.squaredNorm();
// f'' = -theta * f' - p'Mp
Vector cache;fgs.apply_Mv(vecp, cache); // cache = Mp
Scalar fpp = -bfgs.theta() * fp - vecp.dot(cache);

// Theoretical step size to move
Scalar deltatmin = -fp / fpp;

// Limit on the current interval
Scalar il = Scalar(0);
// We have excluded the case that max(brk) <= 0
int b = 0;
Scalar iu = (nord < 1) ? inf : brk[ord[b]];
Scalar deltat = iu - il;

/* int iter = 0;
std::cout << "** Iter " << iter << " **\n";
std::cout << " fp = " << fp << ", fpp = " << fpp << ", deltatmin = " << deltatmin << std::endl;
std::cout << " il = " << il << ", iu = " << iu << ", deltat = " << deltat << std::endl; */
```

--------------------------------

### LineSearch Bracketing Implementation

Source: https://lbfgspp.statr.me/doc/LineSearchBracketing_8h_source.html

Implements a bracketing line search algorithm. It is used to find an appropriate step size in the search direction that satisfies certain conditions, such as Armijo or Wolfe conditions. Ensure that the initial step is positive and the search direction is a descent direction.

```C++
#include <Eigen/Core>
#include <stdexcept> // std::runtime_error
#include <cmath>
#include "Param.h"

namespace LBFGSpp {

template <typename Scalar>
class LineSearchBracketing
{
private:
 using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;

public:
 template <typename Foo>
 static void LineSearch(Foo& f, const LBFGSParam<Scalar>& param,
 const Vector& xp, const Vector& drt, const Scalar& step_max,
 Scalar& step, Scalar& fx, Vector& grad, Scalar& dg, Vector& x)
 {
 // Check the value of step
 if (step <= Scalar(0))
 throw std::invalid_argument("'step' must be positive");

 // Save the function value at the current x
 const Scalar fx_init = fx;
 // Projection of gradient on the search direction
 const Scalar dg_init = grad.dot(drt);
 // Make sure d points to a descent direction
 if (dg_init > 0)
 throw std::logic_error("the moving direction increases the objective function value");

 const Scalar test_decr = param.ftol * dg_init;

 // Upper and lower end of the current line search range
 Scalar step_lo = 0,
 step_hi = std::numeric_limits<Scalar>::infinity();

 int iter;
 for (iter = 0; iter < param.max_linesearch; iter++)
 {
 // x_{k+1} = x_k + step * d_k
 x.noalias() = xp + step * drt;
 // Evaluate this candidate
 fx = f(x, grad);

 if (fx > fx_init + step * test_decr || (!std::isfinite(fx)))
 {
 step_hi = step;
 }
 else
 {
 dg = grad.dot(drt);

 // Armijo condition is met
 if (param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_ARMIJO)
 break;

 if (dg < param.wolfe * dg_init)
 {
 step_lo = step;
 }
 else
 {
 // Regular Wolfe condition is met
 if (param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_WOLFE)
 break;

 if (dg > -param.wolfe * dg_init)
 {
 step_hi = step;
 }
 else
 {
 // Strong Wolfe condition is met
 break;
 }
 }
 }

 if (step_lo > step_hi)
 throw std::runtime_error("the lower bound of the bracketing interval becomes larger than the upper bound");

 if (step < param.min_step)
 throw std::runtime_error("the line search step became smaller than the minimum value allowed");

 if (step > param.max_step)
 throw std::runtime_error("the line search step became larger than the maximum value allowed");

 // continue search in mid of current search range
 step = std::isinf(step_hi) ? 2 * step : step_lo / 2 + step_hi / 2;
 }

 if (iter >= param.max_linesearch)
 throw std::runtime_error("the line search routine reached the maximum number of iterations");
 }
};

} // namespace LBFGSpp

```

--------------------------------

### LBFGSSolver Constructor

Source: https://lbfgspp.statr.me/doc/classLBFGSpp_1_1LBFGSSolver-members.html

Constructs an LBFGSSolver object with the specified parameters.

```APIDOC
## LBFGSSolver Constructor

### Description
Constructs an LBFGSSolver object with the specified parameters.

### Method
Constructor

### Parameters
#### Path Parameters
None

#### Query Parameters
None

#### Request Body
None

### Request Example
```
LBFGSpp::LBFGSSolver<Scalar, LineSearch> solver(param);
```

### Response
#### Success Response
An instance of LBFGSSolver.

#### Response Example
None
```

--------------------------------

### Apply PtBQv Transformation (Matrix)

Source: https://lbfgspp.statr.me/doc/BFGSMat_8h_source.html

Computes P'BQv where Q is represented by an explicit WQ matrix. Performs direct matrix multiplication.

```C++
inline bool apply_PtBQv(const Matrix& WP, const Matrix& WQ, const Vector& v, Vector& res) const
{
    const int nP = WP.rows();
    const int nQ = WQ.rows();
    res.resize(nP);
    if (m_ncorr < 1 || nP < 1 || nQ < 1)
    {
        res.setZero();
        return false;
    }

    // Remember that W = [Y, theta * S], so we need to multiply theta to the second half
    Vector WQtv = WQ.transpose() * v;
    WQtv.tail(m_ncorr) *= m_theta;
    return true;
}
```

--------------------------------

### Nocedal-Wright Line Search Implementation

Source: https://lbfgspp.statr.me/doc/LineSearchNocedalWright_8h_source.html

Performs a line search using the Nocedal-Wright method, ensuring sufficient decrease and curvature conditions (strong Wolfe conditions). Throws exceptions for invalid parameters or non-descent directions.

```cpp
template <typename Foo>
static void LineSearch(Foo& f, const LBFGSParam<Scalar>& param,
 const Vector& xp, const Vector& drt, const Scalar& step_max,
 Scalar& step, Scalar& fx, Vector& grad, Scalar& dg, Vector& x)
{
 using std::abs;

 // Check the value of step
 if (step <= Scalar(0))
 throw std::invalid_argument("'step' must be positive");

 if (param.linesearch != LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE)
 throw std::invalid_argument("'param.linesearch' must be 'LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE' for LineSearchNocedalWright");

 // To make this implementation more similar to the other line search
 // methods in LBFGSpp, the symbol names from the literature
 // ("Numerical Optimizations") have been changed.
 //
 // Literature | LBFGSpp
 // -----------|--------
 // alpha | step
 // phi | fx
 // phi' | dg

 // The expansion rate of the step size
 const Scalar expansion = Scalar(2);

 // Save the function value at the current x
 const Scalar fx_init = fx;
 // Projection of gradient on the search direction
 const Scalar dg_init = dg;
 // Make sure d points to a descent direction
 if (dg_init > Scalar(0))
 throw std::logic_error("the moving direction increases the objective function value");

 const Scalar test_decr = param.ftol * dg_init, // Sufficient decrease
 test_curv = -param.wolfe * dg_init; // Curvature

 // Ends of the line search range (step_lo > step_hi is allowed)
 // We can also define dg_hi, but it will never be used
 Scalar step_hi, fx_hi;
 Scalar step_lo = Scalar(0), fx_lo = fx_init, dg_lo = dg_init;
 // We also need to save x and grad for step=step_lo, since we want to return the best
 // step size along the path when strong Wolfe condition is not met
 Vector x_lo = xp, grad_lo = grad;

 // STEP 1: Bracketing Phase
 // Find a range guaranteed to contain a step satisfying strong Wolfe.
 // The bracketing phase exits if one of the following conditions is satisfied:
 // (1) Current step violates the sufficient decrease condition
 // (2) Current fx >= previous fx
}
```
Tokens:
Raw