From a643ebf037c66d8aa8eb70ba37c2473ba3469aa5 Mon Sep 17 00:00:00 2001 From: Khizir Siddiqui <14289068+khizirsiddiqui@users.noreply.github.com> Date: Wed, 18 Mar 2026 07:42:15 +0100 Subject: [PATCH] adds L-BFGS-B (L-BFGS Bounded) --- include/ensmallen.hpp | 1 + include/ensmallen_bits/lbfgsb/lbfgsb.hpp | 312 +++++++++++++++ include/ensmallen_bits/lbfgsb/lbfgsb_impl.hpp | 360 ++++++++++++++++++ tests/CMakeLists.txt | 1 + tests/lbfgsb_test.cpp | 101 +++++ 5 files changed, 775 insertions(+) create mode 100644 include/ensmallen_bits/lbfgsb/lbfgsb.hpp create mode 100644 include/ensmallen_bits/lbfgsb/lbfgsb_impl.hpp create mode 100644 tests/lbfgsb_test.cpp diff --git a/include/ensmallen.hpp b/include/ensmallen.hpp index 176407d1d..075346978 100644 --- a/include/ensmallen.hpp +++ b/include/ensmallen.hpp @@ -133,6 +133,7 @@ #include "ensmallen_bits/iqn/iqn.hpp" #include "ensmallen_bits/katyusha/katyusha.hpp" #include "ensmallen_bits/lbfgs/lbfgs.hpp" +#include "ensmallen_bits/lbfgsb/lbfgsb.hpp" #include "ensmallen_bits/lookahead/lookahead.hpp" #include "ensmallen_bits/agemoea/agemoea.hpp" #include "ensmallen_bits/moead/moead.hpp" diff --git a/include/ensmallen_bits/lbfgsb/lbfgsb.hpp b/include/ensmallen_bits/lbfgsb/lbfgsb.hpp new file mode 100644 index 000000000..9f3c0f151 --- /dev/null +++ b/include/ensmallen_bits/lbfgsb/lbfgsb.hpp @@ -0,0 +1,312 @@ +/** + * @file lbfgsb.hpp + * @author Khizir Siddiqui + * + * The generic L-BFGS-B optimizer. + * + * ensmallen is free software; you may redistribute it and/or modify it under + * the terms of the 3-clause BSD license. You should have received a copy of + * the 3-clause BSD license along with ensmallen. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ +#ifndef ENSMALLEN_LBFGSB_LBFGSB_HPP +#define ENSMALLEN_LBFGSB_LBFGSB_HPP + +#include + +namespace ens { + +/** + * The L-BFGS-B optimizer, which is an extension of L-BFGS designed to handle + * simple bound constraints. It uses gradient projection to identify the active + * of constraints, and then performs minimization using the L-BFGS approximation + * to the Hessian. + * + * L_BFGS_B can optimize differentiable functions with box constraints. + * For more details, see the documentation on function types included with this + * distribution or on the ensmallen website. + * + * The algorithm is based on the paper: + * "A Limited Memory Algorithm for Bound Constrained Optimization" + * by R.H. Byrd and P. Lu and J. Nocedal + */ +class L_BFGS_B +{ + public: + /** + * Initialize the L-BFGS-B object. There are many parameters that can be set + * for the optimization, but default values are given for each of them. + * + * @param numBasis Number of memory points to be stored (default 5). + * @param lowerBound Lower bound for the coordinates (can be a scalar or matrix). + * @param upperBound Upper bound for the coordinates (can be a scalar or matrix). + * @param maxIterations Maximum number of iterations for the optimization + * (0 means no limit and may run indefinitely). + * @param armijoConstant Controls the accuracy of the line search routine for + * determining the Armijo condition. + * @param wolfe Parameter for detecting the Wolfe condition. + * @param minGradientNorm Minimum gradient norm required to continue the + * optimization. + * @param factr Minimum relative function value decrease to continue + * the optimization. + * @param maxLineSearchTrials The maximum number of trials for the line search + * (before giving up). + * @param minStep The minimum step of the line search. + * @param maxStep The maximum step of the line search. + */ + L_BFGS_B(const size_t numBasis = 10, /* same default as scipy */ + const arma::mat& lowerBound = arma::mat(), + const arma::mat& upperBound = arma::mat(), + const size_t maxIterations = 10000, /* many but not infinite */ + const double armijoConstant = 1e-4, + const double wolfe = 0.9, + const double minGradientNorm = 1e-6, + const double factr = 1e-15, + const size_t maxLineSearchTrials = 50, + const double minStep = 1e-20, + const double maxStep = 1e20); + + /** + * Initialize the L-BFGS-B object with scalar bounds. + */ + L_BFGS_B(const size_t numBasis, + const double lowerBound, + const double upperBound, + const size_t maxIterations = 10000, + const double armijoConstant = 1e-4, + const double wolfe = 0.9, + const double minGradientNorm = 1e-6, + const double factr = 1e-15, + const size_t maxLineSearchTrials = 50, + const double minStep = 1e-20, + const double maxStep = 1e20); + + /** + * Use L-BFGS-B to optimize the given function, starting at the given iterate + * point and finding the minimum. The maximum number of iterations is set in + * the constructor (or with MaxIterations()). Alternately, another overload + * is provided which takes a maximum number of iterations as a parameter. The + * given starting point will be modified to store the finishing point of the + * algorithm, and the final objective value is returned. + * + * @tparam FunctionType Type of the function to be optimized. + * @tparam MatType Type of matrix to optimize with. + * @tparam GradType Type of matrix to use to represent function gradients. + * @tparam CallbackTypes Types of callback functions. + * @param function Function to optimize; must have Evaluate() and Gradient(). + * @param iterate Starting point (will be modified). + * @param callbacks Callback functions. + * @return Objective value of the final point. + */ + template + typename std::enable_if::value, + typename MatType::elem_type>::type + Optimize(FunctionType& function, + MatType& iterate, + CallbackTypes&&... callbacks); + + //! Forward the MatType as GradType. + template + typename MatType::elem_type Optimize(SeparableFunctionType& function, + MatType& iterate, + CallbackTypes&&... callbacks) + { + return Optimize(function, iterate, + std::forward(callbacks)...); + } + + //! Get the memory size. + size_t NumBasis() const { return numBasis; } + //! Modify the memory size. + size_t& NumBasis() { return numBasis; } + + //! Get the lower bound. + const arma::mat& LowerBound() const { return lowerBound; } + //! Modify the lower bound. + arma::mat& LowerBound() { return lowerBound; } + + //! Get the upper bound. + const arma::mat& UpperBound() const { return upperBound; } + //! Modify the upper bound. + arma::mat& UpperBound() { return upperBound; } + + //! Get the maximum number of iterations. + size_t MaxIterations() const { return maxIterations; } + //! Modify the maximum number of iterations. + size_t& MaxIterations() { return maxIterations; } + + //! Get the Armijo condition constant. + double ArmijoConstant() const { return armijoConstant; } + //! Modify the Armijo condition constant. + double& ArmijoConstant() { return armijoConstant; } + + //! Get the Wolfe parameter. + double Wolfe() const { return wolfe; } + //! Modify the Wolfe parameter. + double& Wolfe() { return wolfe; } + + //! Get the minimum gradient norm. + double MinGradientNorm() const { return minGradientNorm; } + //! Modify the minimum gradient norm. + double& MinGradientNorm() { return minGradientNorm; } + + //! Get the factr value. + double Factr() const { return factr; } + //! Modify the factr value. + double& Factr() { return factr; } + + //! Get the maximum number of line search trials. + size_t MaxLineSearchTrials() const { return maxLineSearchTrials; } + //! Modify the maximum number of line search trials. + size_t& MaxLineSearchTrials() { return maxLineSearchTrials; } + + //! Return the minimum line search step size. + double MinStep() const { return minStep; } + //! Modify the minimum line search step size. + double& MinStep() { return minStep; } + + //! Return the maximum line search step size. + double MaxStep() const { return maxStep; } + //! Modify the maximum line search step size. + double& MaxStep() { return maxStep; } + + private: + //! Size of memory for this L-BFGS-B optimizer. + size_t numBasis; + //! Lower bound for coordinates. + arma::mat lowerBound; + //! Upper bound for coordinates. + arma::mat upperBound; + //! Maximum number of iterations. + size_t maxIterations; + //! Parameter for determining the Armijo condition. + double armijoConstant; + //! Parameter for detecting the Wolfe condition. + double wolfe; + //! Minimum gradient norm required to continue the optimization. + double minGradientNorm; + //! Minimum relative function value decrease to continue the optimization. + double factr; + //! Maximum number of trials for the line search. + size_t maxLineSearchTrials; + //! Minimum step of the line search. + double minStep; + //! Maximum step of the line search. + double maxStep; + //! Controls early termination of the optimization process. + bool terminate; + //! Flag indicating whether scalar bounds were provided. + bool usingScalarBounds; + + /** + * Project the given point onto the bounds. + */ + template + void ProjectPoint(MatType& iterate); + + /** + * Find the generalized Cauchy point. + * + * @param iterate The current point. + * @param gradient The gradient at the current point. + * @param theta The scaling factor from the L-BFGS matrix. + * @param W The W matrix from the L-BFGS approximation. + * @param M The M matrix from the L-BFGS approximation. + * @param cauchyPoint Vector to store the resulting Cauchy point. + * @param c Vector to store the c vector used in subspace minimization. + * @param activeSet Boolean vector specifying if a coordinate is at bound. + */ + template + void GeneralizedCauchyPoint(const MatType& iterate, + const MatType& gradient, + const typename MatType::elem_type theta, + const arma::mat& W, + const arma::mat& M, + MatType& cauchyPoint, + arma::vec& c, + arma::uvec& activeSet); + + /** + * Perform subspace minimization over the free variables. + * + * @param iterate The current point. + * @param cauchyPoint The generalized Cauchy point. + * @param gradient The gradient at the current point. + * @param theta The scaling factor from the L-BFGS matrix. + * @param W The W matrix from the L-BFGS approximation. + * @param M The M matrix from the L-BFGS approximation. + * @param c The c vector from the Cauchy point computation. + * @param activeSet Boolean vector specifying if a coordinate is at bound. + * @param searchDirection Vector to store the resulting search direction. + */ + template + void SubspaceMinimization(const MatType& iterate, + const MatType& cauchyPoint, + const MatType& gradient, + const typename MatType::elem_type theta, + const arma::mat& W, + const arma::mat& M, + const arma::vec& c, + const arma::uvec& activeSet, + MatType& searchDirection); + + /** + * Reconstruct the L-BFGS memory representation into W and M matrices. + * + * @param iterationNum The iteration number. + * @param s Differences between the iterate and old iterate matrix. + * @param y Differences between the gradient and the old gradient matrix. + * @param theta Scaling factor to use. + * @param W The resulting W matrix. + * @param M The resulting M matrix. + */ + template + void ComputeWM(const size_t iterationNum, + const CubeType& s, + const CubeType& y, + const typename CubeType::elem_type theta, + arma::mat& W, + arma::mat& M); + + /** + * Perform a projected back-tracking line search along the search direction + * to calculate a step size satisfying the Wolfe conditions. + * + * @param function Function to optimize. + * @param functionValue Value of the function at the initial point. + * @param iterate The initial point to begin the line search from. + * @param gradient The gradient at the initial point. + * @param searchDirection A vector specifying the search direction. + * @param finalStepSize The resulting step size (0 if no step). + * @param callbacks Callback functions. + * + * @return false if no step size is suitable, true otherwise. + */ + template + bool LineSearch(FunctionType& function, + ElemType& functionValue, + MatType& iterate, + GradType& gradient, + MatType& newIterateTmp, + const GradType& searchDirection, + ElemType& finalStepSize, + CallbackTypes&... callbacks); +}; + +using LBFGSB = L_BFGS_B; + +} // namespace ens + +#include "lbfgsb_impl.hpp" + +#endif // ENSMALLEN_LBFGSB_LBFGSB_HPP diff --git a/include/ensmallen_bits/lbfgsb/lbfgsb_impl.hpp b/include/ensmallen_bits/lbfgsb/lbfgsb_impl.hpp new file mode 100644 index 000000000..42ae1c4a3 --- /dev/null +++ b/include/ensmallen_bits/lbfgsb/lbfgsb_impl.hpp @@ -0,0 +1,360 @@ +/** + * @file lbfgsb_impl.hpp + * @author Khizir Siddiqui + * + * The implementation of the L_BFGS_B optimizer. + * + * ensmallen is free software; you may redistribute it and/or modify it under + * the terms of the 3-clause BSD license. You should have received a copy of + * the 3-clause BSD license along with ensmallen. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ +#ifndef ENSMALLEN_LBFGSB_LBFGSB_IMPL_HPP +#define ENSMALLEN_LBFGSB_LBFGSB_IMPL_HPP + +// In case it hasn't been included yet. +#include "lbfgsb.hpp" + +#include + +namespace ens { + +inline L_BFGS_B::L_BFGS_B(const size_t numBasis, + const arma::mat& lowerBound, + const arma::mat& upperBound, + const size_t maxIterations, + const double armijoConstant, + const double wolfe, + const double minGradientNorm, + const double factr, + const size_t maxLineSearchTrials, + const double minStep, + const double maxStep) : + numBasis(numBasis), + lowerBound(lowerBound), + upperBound(upperBound), + maxIterations(maxIterations), + armijoConstant(armijoConstant), + wolfe(wolfe), + minGradientNorm(minGradientNorm), + factr(factr), + maxLineSearchTrials(maxLineSearchTrials), + minStep(minStep), + maxStep(maxStep), + terminate(false), + usingScalarBounds(false) +{ + // Nothing to do. +} + +inline L_BFGS_B::L_BFGS_B(const size_t numBasis, + const double lowerBound, + const double upperBound, + const size_t maxIterations, + const double armijoConstant, + const double wolfe, + const double minGradientNorm, + const double factr, + const size_t maxLineSearchTrials, + const double minStep, + const double maxStep) : + numBasis(numBasis), + lowerBound(arma::vec(1).fill(lowerBound)), + upperBound(arma::vec(1).fill(upperBound)), + maxIterations(maxIterations), + armijoConstant(armijoConstant), + wolfe(wolfe), + minGradientNorm(minGradientNorm), + factr(factr), + maxLineSearchTrials(maxLineSearchTrials), + minStep(minStep), + maxStep(maxStep), + terminate(false), + usingScalarBounds(true) +{ + // Nothing to do. +} + +template +void L_BFGS_B::ProjectPoint(MatType& iterate) +{ + if (usingScalarBounds) + { + iterate.clamp(lowerBound(0), upperBound(0)); + } + else if (!lowerBound.is_empty() && !upperBound.is_empty()) + { + for (size_t i = 0; i < iterate.n_elem; ++i) + { + if (iterate[i] < lowerBound[i]) + iterate[i] = lowerBound[i]; + else if (iterate[i] > upperBound[i]) + iterate[i] = upperBound[i]; + } + } +} + +template +typename std::enable_if::value, + typename MatType::elem_type>::type +L_BFGS_B::Optimize(FunctionType& function, + MatType& iterateIn, + CallbackTypes&&... callbacks) +{ + typedef typename MatType::elem_type ElemType; + typedef typename MatTypeTraits::BaseMatType BaseMatType; + typedef typename MatTypeTraits::BaseMatType BaseGradType; + typedef Function FullFunctionType; + + FullFunctionType& f = static_cast(function); + traits::CheckFunctionTypeAPI(); + RequireFloatingPointType(); + RequireFloatingPointType(); + RequireSameInternalTypes(); + + BaseMatType& iterate = (BaseMatType&) iterateIn; + const size_t rows = iterate.n_rows; + const size_t cols = iterate.n_cols; + + // Verify bound matrix sizes if not using scalar bounds. + if (!usingScalarBounds) + { + if (!lowerBound.is_empty() && (lowerBound.n_rows != rows || lowerBound.n_cols != cols)) + throw std::invalid_argument("L_BFGS_B: lowerBound matrix size does not match iterate matrix size."); + if (!upperBound.is_empty() && (upperBound.n_rows != rows || upperBound.n_cols != cols)) + throw std::invalid_argument("L_BFGS_B: upperBound matrix size does not match iterate matrix size."); + } + + ProjectPoint(iterate); + + BaseMatType newIterateTmp(rows, cols); + typedef typename ForwardType::bcube BaseCubeType; + BaseCubeType s(rows, cols, numBasis); + BaseCubeType y(rows, cols, numBasis); + + BaseMatType oldIterate(rows, cols); + oldIterate.zeros(); + bool optimizeUntilConvergence = (maxIterations == 0); + + BaseGradType gradient(rows, cols); + gradient.zeros(); + BaseGradType oldGradient(rows, cols); + oldGradient.zeros(); + BaseGradType searchDirection(rows, cols); + searchDirection.zeros(); + + typedef typename ForwardType::bmat DenseMatType; + + ElemType functionValue = f.EvaluateWithGradient(iterate, gradient); + terminate |= Callback::EvaluateWithGradient(*this, f, iterate, + functionValue, gradient, callbacks...); + + ElemType prevFunctionValue; + + // L-BFGS-B specific variables. + arma::mat W, M; + size_t memorySlotsUsed = 0; + + Callback::BeginOptimization(*this, f, iterate, callbacks...); + for (size_t itNum = 0; (optimizeUntilConvergence || (itNum != maxIterations)) + && !terminate; ++itNum) + { + prevFunctionValue = functionValue; + + // Check convergence: projected gradient norm + // L-BFGS-B uses the infinity norm of the projected gradient. + ElemType projectedGradNorm = 0.0; + for (size_t i = 0; i < iterate.n_elem; ++i) + { + ElemType lb = usingScalarBounds ? lowerBound(0) : + (!lowerBound.is_empty() ? lowerBound[i] : -std::numeric_limits::infinity()); + ElemType ub = usingScalarBounds ? upperBound(0) : + (!upperBound.is_empty() ? upperBound[i] : std::numeric_limits::infinity()); + + ElemType g = gradient[i]; + ElemType pgrad; + if (iterate[i] <= lb + 1e-12 && g > 0) + pgrad = 0; + else if (iterate[i] >= ub - 1e-12 && g < 0) + pgrad = 0; + else + pgrad = g; + + projectedGradNorm = std::max(projectedGradNorm, std::abs(pgrad)); + } + + if (projectedGradNorm < minGradientNorm) + { + Info << "L-BFGS-B: projected gradient norm too small (terminating successfully)." + << std::endl; + break; + } + + if (std::isnan(functionValue)) + { + Warn << "L-BFGS-B: objective value is NaN (terminating)!" << std::endl; + break; + } + + // Determine theta (L-BFGS scaling factor). + // \theta = \frac{y_k^T y_k}{s_k^T y_k} + ElemType theta = 1.0; + if (itNum > 0) + { + int previousPos = (itNum - 1) % numBasis; + const DenseMatType& sMat = s.slice(previousPos); + const DenseMatType& yMat = y.slice(previousPos); + ElemType sy = arma::dot(sMat, yMat); + ElemType yy = arma::dot(yMat, yMat); + if (sy > 1e-10 * yy) + theta = yy / sy; + } + + // 1. Cauchy Point Computation + // The generalized Cauchy point x^c is defined as the first local minimizer + // of the quadratic model along the projected gradient direction P(x - t * g, l, u). + // Here we use an active-set projection approach over the free variables. + searchDirection = gradient; + arma::uvec isFree(iterate.n_elem, arma::fill::ones); + + for (size_t i = 0; i < iterate.n_elem; ++i) + { + ElemType lb = usingScalarBounds ? lowerBound(0) : + (!lowerBound.is_empty() ? lowerBound[i] : -std::numeric_limits::infinity()); + ElemType ub = usingScalarBounds ? upperBound(0) : + (!upperBound.is_empty() ? upperBound[i] : std::numeric_limits::infinity()); + + if ((iterate[i] <= lb + 1e-12 && gradient[i] > 0) || + (iterate[i] >= ub - 1e-12 && gradient[i] < 0)) + { + searchDirection[i] = 0; + isFree[i] = 0; + } + } + + // 2. Reconstruct SearchDirection using the Subspace Minimization step. + // L-BFGS two-loop recursion over the subspace of free variables to compute H_k \nabla f(x_k). + if (memorySlotsUsed > 0) + { + arma::vec rho(numBasis, arma::fill::zeros); + arma::vec alpha(numBasis, arma::fill::zeros); + size_t limit = (numBasis > itNum) ? 0 : (itNum - numBasis); + + // First loop: \alpha_i = \rho_i s_i^T q, q = q - \alpha_i y_i + for (size_t i = itNum; i != limit; i--) + { + int pos = (i + (numBasis - 1)) % numBasis; + // Zero out bound variables from s and y computations to keep them in subspace + DenseMatType sFree = s.slice(pos) % isFree; + DenseMatType yFree = y.slice(pos) % isFree; + + ElemType sy = arma::dot(yFree, sFree); + rho[itNum - i] = (sy != 0) ? (1.0 / sy) : 1.0; + alpha[itNum - i] = rho[itNum - i] * arma::dot(sFree, searchDirection); + searchDirection -= alpha[itNum - i] * yFree; + } + + searchDirection /= theta; + + // Second loop: \beta_i = \rho_i y_i^T r, r = r + (\alpha_i - \beta_i) s_i + for (size_t i = limit; i < itNum; i++) + { + int pos = i % numBasis; + DenseMatType sFree = s.slice(pos) % isFree; + DenseMatType yFree = y.slice(pos) % isFree; + + ElemType beta = rho[itNum - i - 1] * arma::dot(yFree, searchDirection); + searchDirection += (alpha[itNum - i - 1] - beta) * sFree; + } + } + + // Fallback: if search direction isn't a descent direction, use steepest descent. + if (arma::dot(searchDirection, gradient) <= 0) + searchDirection = gradient % isFree; + + searchDirection *= -1; + + oldIterate = iterate; + oldGradient = gradient; + + // 4. Projected Backtracking Line Search + // Find step size \alpha that satisfies the Armijo condition for the projected point: + // f(P(x_k + \alpha d_k, l, u)) \le f(x_k) + c_1 \nabla f(x_k)^T (P(x_k + \alpha d_k, l, u) - x_k) + ElemType stepSize = 1.0; + if (stepSize > ElemType(maxStep)) stepSize = ElemType(maxStep); + if (stepSize < ElemType(minStep)) stepSize = ElemType(minStep); + + ElemType initialSearchDirectionDotGradient = arma::dot(gradient, searchDirection); + ElemType initialFunctionValue = functionValue; + + size_t numIterations = 0; + const ElemType dec = 0.5; + + while (true) + { + newIterateTmp = iterate + stepSize * searchDirection; + ProjectPoint(newIterateTmp); + + functionValue = f.EvaluateWithGradient(newIterateTmp, gradient); + + if (std::isnan(functionValue)) + { + Warn << "L-BFGS-B: objective value is NaN!" << std::endl; + break; + } + + terminate |= Callback::EvaluateWithGradient(*this, f, newIterateTmp, + functionValue, gradient, callbacks...); + + // Armijo condition check + // directional derivative = \nabla f(x_k)^T (x_{k+1} - x_k) + // For projected search, this is + ElemType actualMoveDotGrad = arma::dot(oldGradient, newIterateTmp - iterate); + + if (functionValue <= initialFunctionValue + armijoConstant * actualMoveDotGrad) + break; + + if (stepSize < minStep || numIterations >= maxLineSearchTrials) + break; + + stepSize *= dec; + numIterations++; + } + + if (arma::norm(newIterateTmp - iterate, "inf") == 0) + { + Info << "L-BFGS-B: step size is effectively zero (terminating successfully)." << std::endl; + break; + } + + iterate = newIterateTmp; + + const ElemType denom = std::max(ElemType(1), + std::max(std::abs(prevFunctionValue), std::abs(functionValue))); + if ((prevFunctionValue - functionValue) / denom <= factr) + { + Info << "L-BFGS-B: function value stable (terminating successfully)." << std::endl; + break; + } + + // 5. Update L-BFGS basis matrices + // s_k = x_{k+1} - x_k + // y_k = \nabla f(x_{k+1}) - \nabla f(x_k) + int overwritePos = itNum % numBasis; + s.slice(overwritePos) = iterate - oldIterate; + y.slice(overwritePos) = gradient - oldGradient; + memorySlotsUsed = std::min(numBasis, memorySlotsUsed + 1); + + terminate |= Callback::StepTaken(*this, f, iterate, callbacks...); + } + + Callback::EndOptimization(*this, f, iterate, callbacks...); + return functionValue; +} + +} // namespace ens + +#endif // ENSMALLEN_LBFGSB_LBFGSB_IMPL_HPP diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b8acfb1c6..46706588a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -31,6 +31,7 @@ set(ENSMALLEN_TESTS_SOURCES indicators_test.cpp katyusha_test.cpp lbfgs_test.cpp + lbfgsb_test.cpp line_search_test.cpp lookahead_test.cpp lrsdp_test.cpp diff --git a/tests/lbfgsb_test.cpp b/tests/lbfgsb_test.cpp new file mode 100644 index 000000000..150765656 --- /dev/null +++ b/tests/lbfgsb_test.cpp @@ -0,0 +1,101 @@ +/** + * @file lbfgsb_test.cpp + * + * ensmallen is free software; you may redistribute it and/or modify it under + * the terms of the 3-clause BSD license. You should have received a copy of + * the 3-clause BSD license along with ensmallen. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ +#if defined(ENS_USE_COOT) + #include + #include +#endif +#include +#include "catch.hpp" +#include "test_function_tools.hpp" + +using namespace ens; +using namespace ens::test; + +TEMPLATE_TEST_CASE("LBFGSB_Unbounded_RosenbrockFunction", "[LBFGSB]", + ENS_FULLPREC_TEST_TYPES, ENS_SPARSE_TEST_TYPES) +{ + typedef typename TestType::elem_type ElemType; + + L_BFGS_B lbfgsb; + lbfgsb.MaxIterations() = 10000; + + // No bounds set, should behave like L-BFGS + RosenbrockFunction f; + + TestType coords = f.GetInitialPoint(); + lbfgsb.Optimize(f, coords); + + ElemType finalValue = f.Evaluate(coords); + + REQUIRE(finalValue == Approx(0.0).margin(Tolerances::Obj)); + REQUIRE(coords(0) == Approx(1.0).epsilon(Tolerances::Coord)); + REQUIRE(coords(1) == Approx(1.0).epsilon(Tolerances::Coord)); +} + +TEMPLATE_TEST_CASE("LBFGSB_Bounded_RosenbrockFunction", "[LBFGSB]", + ENS_FULLPREC_TEST_TYPES) +{ + typedef typename TestType::elem_type ElemType; + + // Set lower bound to [2.0, 2.0] + arma::mat lowerBound(2, 1); + lowerBound(0, 0) = 2.0; + lowerBound(1, 0) = 2.0; + + // Set upper bound to [10.0, 10.0] + arma::mat upperBound(2, 1); + upperBound(0, 0) = 10.0; + upperBound(1, 0) = 10.0; + + L_BFGS_B lbfgsb(10, lowerBound, upperBound); + lbfgsb.MaxIterations() = 10000; + + RosenbrockFunction f; + + // Initial point inside bounds + TestType coords(2, 1); + coords(0, 0) = 2.5; + coords(1, 0) = 3.0; + + lbfgsb.Optimize(f, coords); + + // The true constrained minimum of Rosenbrock for x >= 2, y >= 2 is at: + // x = 2, y = 4 (since it minimizes (1-x)^2 + 100(y-x^2)^2). + // The function value there is (1-2)^2 + 100(4-4)^2 = 1.0. + ElemType finalValue = f.Evaluate(coords); + + REQUIRE(finalValue == Approx(1.0).margin(1e-4)); + REQUIRE(coords(0) == Approx(2.0).margin(1e-3)); + REQUIRE(coords(1) == Approx(4.0).margin(1e-3)); +} + +TEMPLATE_TEST_CASE("LBFGSB_ScalarBounded_RosenbrockFunction", "[LBFGSB]", + ENS_FULLPREC_TEST_TYPES) +{ + typedef typename TestType::elem_type ElemType; + + // Set scalar bounds 2.0 <= x <= 10.0 + L_BFGS_B lbfgsb(10, 2.0, 10.0); + lbfgsb.MaxIterations() = 10000; + + RosenbrockFunction f; + + // Initial point inside bounds + TestType coords(2, 1); + coords(0, 0) = 2.5; + coords(1, 0) = 3.0; + + lbfgsb.Optimize(f, coords); + + ElemType finalValue = f.Evaluate(coords); + + REQUIRE(finalValue == Approx(1.0).margin(1e-4)); + REQUIRE(coords(0) == Approx(2.0).margin(1e-3)); + REQUIRE(coords(1) == Approx(4.0).margin(1e-3)); +}