Derivative-free Broyden¶
Table of contents
Algorithm Description¶
Li and Fukushima (2000) is a derivative-free variant of Broyden’s method for solving systems of nonlinear equations.
We seek a vector of values \(x^{(*)}\) such that
where \(F : \mathbb{R}^n \to \mathbb{R}^m\) is convex and differentiable. The algorithm uses an approximation to the Jacobian.
The updating rule for Broyden’s method is described below. Let \(x^{(i)}\) denote the function input values at stage \(i\) of the algorithm.
Compute the descent direction using:
\[d^{(i)} = - B^{(i)} F(x^{(i)})\]Quasi line search step.
If
\[\| F(x^{(i)} + d^{(i)}) \| \leq \rho \| F(x^{(i)} \| - \sigma_2 \| d^{(i)} \|^2\]then set \(\lambda_i = 1.0\); otherwise set \(\eta_i = 1.0 / i^2\) and find the smallest \(k\) for which
\[\| F(x^{(i)} + \lambda_i d^{(i)}) \| \leq \| F(x^{(i)} \| - \sigma_1 \| d^{(i)} \|^2 + \eta_i \| F( x^{(i)} ) \|\]holds with \(\lambda_i = \beta^k\).
Update the candidate solution vector using:
Update the approximate inverse Jacobian matrix, \(B\), using:
\[B^{(i+1)} = B^{(i)} + \frac{1}{[y^{(i+1)}]^\top y^{(i+1)}} (s^{(i+1)} - B^{(i)} y^{(i+1)}) [y^{(i+1)}]^\top\]where
\[\begin{aligned} s^{(i)} &:= x^{(i)} - x^{(i-1)} \\ y^{(i)} &:= F(x^{(i)}) - F(x^{(i-1)}) \end{aligned}\]
The algorithm stops when at least one of the following conditions are met:
\(\| F \|\) is less than
rel_objfn_change_tol
.the relative change between \(x^{(i+1)}\) and \(x^{(i)}\) is less than
rel_sol_change_tol
;the total number of iterations exceeds
iter_max
.
Function Declarations¶
-
bool broyden_df(ColVec_t &init_out_vals, std::function<ColVec_t(const ColVec_t &vals_inp, void *opt_data)> opt_objfn, void *opt_data)¶
Derivative-free variant of Broyden’s method due to Li and Fukushima (2000)
- Parameters
init_out_vals – a column vector of initial values, which will be replaced by the solution upon successful completion of the optimization algorithm.
opt_objfn – the function to be minimized, taking three arguments:
vals_inp
a vector of inputs; andopt_data
additional data passed to the user-provided function.
opt_data – additional data passed to the user-provided function.
- Returns
a boolean value indicating successful completion of the optimization algorithm.
-
bool broyden_df(ColVec_t &init_out_vals, std::function<ColVec_t(const ColVec_t &vals_inp, void *opt_data)> opt_objfn, void *opt_data, algo_settings_t &settings)¶
Derivative-free variant of Broyden’s method due to Li and Fukushima (2000)
- Parameters
init_out_vals – a column vector of initial values, which will be replaced by the solution upon successful completion of the optimization algorithm.
opt_objfn – the function to be minimized, taking three arguments:
vals_inp
a vector of inputs; andopt_data
additional data passed to the user-provided function.
opt_data – additional data passed to the user-provided function.
settings – parameters controlling the optimization routine.
- Returns
a boolean value indicating successful completion of the optimization algorithm.
-
bool broyden_df(ColVec_t &init_out_vals, std::function<ColVec_t(const ColVec_t &vals_inp, void *opt_data)> opt_objfn, void *opt_data, std::function<Mat_t(const ColVec_t &vals_inp, void *jacob_data)> jacob_objfn, void *jacob_data)¶
Derivative-free variant of Broyden’s method due to Li and Fukushima (2000)
- Parameters
init_out_vals – a column vector of initial values, which will be replaced by the solution upon successful completion of the optimization algorithm.
opt_objfn – the function to be minimized, taking three arguments:
vals_inp
a vector of inputs; andopt_data
additional data passed to the user-provided function.
opt_data – additional data passed to the user-provided function.
jacob_objfn – a function to calculate the Jacobian matrix, taking two arguments:
vals_inp
a vector of inputs; andjacob_data
additional data passed to the Jacobian function.
jacob_data – additional data passed to the Jacobian function.
- Returns
a boolean value indicating successful completion of the optimization algorithm.
-
bool broyden_df(ColVec_t &init_out_vals, std::function<ColVec_t(const ColVec_t &vals_inp, void *opt_data)> opt_objfn, void *opt_data, std::function<Mat_t(const ColVec_t &vals_inp, void *jacob_data)> jacob_objfn, void *jacob_data, algo_settings_t &settings)¶
Derivative-free variant of Broyden’s method due to Li and Fukushima (2000)
- Parameters
init_out_vals – a column vector of initial values, which will be replaced by the solution upon successful completion of the optimization algorithm.
opt_objfn – the function to be minimized, taking three arguments:
vals_inp
a vector of inputs; andopt_data
additional data passed to the user-provided function.
opt_data – additional data passed to the user-provided function.
jacob_objfn – a function to calculate the Jacobian matrix, taking two arguments:
vals_inp
a vector of inputs; andjacob_data
additional data passed to the Jacobian function.
jacob_data – additional data passed to the Jacobian function.
settings – parameters controlling the optimization routine.
- Returns
a boolean value indicating successful completion of the optimization algorithm.
Optimization Control Parameters¶
The basic control parameters are:
fp_t rel_objfn_change_tol
: the error tolerance value controlling how small \(\| F \|\) should be before ‘convergence’ is declared.fp_t rel_sol_change_tol
: the error tolerance value controlling how small the proportional change in the solution vector should be before ‘convergence’ is declared.The relative change is computed using:
\[\left\| \dfrac{x^{(i)} - x^{(i-1)}}{ |x^{(i-1)}| + \epsilon } \right\|_1\]where \(\epsilon\) is a small number added for numerical stability.
size_t iter_max
: the maximum number of iterations/updates before the algorithm exits.bool vals_bound
: whether the search space of the algorithm is bounded. Iftrue
, thenColVec_t lower_bounds
: defines the lower bounds of the search space.ColVec_t upper_bounds
: defines the upper bounds of the search space.
In addition to these:
int print_level
: Set the level of detail for printing updates on optimization progress.Level
0
: Nothing (default).Level
1
: Print the current iteration count and error values.Level
2
: Level 1 plus the current candidate solution values, \(x^{(i+1)}\).Level
3
: Level 2 plus the direction vector, \(d^{(i)}\), and the function values, \(F(x^{(i+1)})\).Level
4
: Level 3 plus the components used to update the approximate inverse Jacobian matrix: \(s^{(i+1)}\), \(y^{(i+1)}\), and \(B^{(i+1)}\).
Examples¶
Example 1¶
Code to run this example is given below.
Armadillo (Click to show/hide)
#define OPTIM_ENABLE_ARMA_WRAPPERS
#include "optim.hpp"
inline
arma::vec
zeros_test_objfn_1(const arma::vec& vals_inp, void* opt_data)
{
double x_1 = vals_inp(0);
double x_2 = vals_inp(1);
//
arma::vec ret(2);
ret(0) = std::exp(-std::exp(-(x_1+x_2))) - x_2*(1 + std::pow(x_1,2));
ret(1) = x_1*std::cos(x_2) + x_2*std::sin(x_1) - 0.5;
//
return ret;
}
inline
arma::mat
zeros_test_jacob_1(const arma::vec& vals_inp, void* opt_data)
{
double x_1 = vals_inp(0);
double x_2 = vals_inp(1);
//
arma::mat ret(2,2);
ret(0,0) = std::exp(-std::exp(-(x_1+x_2))-(x_1+x_2)) - 2*x_1*x_1;
ret(0,1) = std::exp(-std::exp(-(x_1+x_2))-(x_1+x_2)) - x_1*x_1 - 1.0;
ret(1,0) = std::cos(x_2) + x_2*std::cos(x_1);
ret(1,1) = -x_1*std::sin(x_2) + std::cos(x_1);
//
return ret;
}
int main()
{
arma::vec x = arma::zeros(2,1); // initial values (0,0)
bool success = optim::broyden_df(x, zeros_test_objfn_1, nullptr);
if (success) {
std::cout << "broyden_df: test_1 completed successfully." << "\n";
} else {
std::cout << "broyden_df: test_1 completed unsuccessfully." << "\n";
}
arma::cout << "broyden_df: solution to test_1:\n" << x << arma::endl;
//
x = arma::zeros(2,1);
success = optim::broyden_df(x, zeros_test_objfn_1, nullptr, zeros_test_jacob_1, nullptr);
if (success) {
std::cout << "broyden_df with jacobian: test_1 completed successfully." << "\n";
} else {
std::cout << "broyden_df with jacobian: test_1 completed unsuccessfully." << "\n";
}
arma::cout << "broyden_df with jacobian: solution to test_1:\n" << x << arma::endl;
//
return 0;
}
Eigen (Click to show/hide)
#define OPTIM_ENABLE_EIGEN_WRAPPERS
#include "optim.hpp"
inline
Eigen::VectorXd
zeros_test_objfn_1(const Eigen::VectorXd& vals_inp, void* opt_data)
{
double x_1 = vals_inp(0);
double x_2 = vals_inp(1);
//
Eigen::VectorXd ret(2);
ret(0) = std::exp(-std::exp(-(x_1+x_2))) - x_2*(1 + std::pow(x_1,2));
ret(1) = x_1*std::cos(x_2) + x_2*std::sin(x_1) - 0.5;
//
return ret;
}
inline
Eigen::MatrixXd
zeros_test_jacob_1(const Eigen::VectorXd& vals_inp, void* opt_data)
{
double x_1 = vals_inp(0);
double x_2 = vals_inp(1);
//
Eigen::MatrixXd ret(2,2);
ret(0,0) = std::exp(-std::exp(-(x_1+x_2))-(x_1+x_2)) - 2*x_1*x_1;
ret(0,1) = std::exp(-std::exp(-(x_1+x_2))-(x_1+x_2)) - x_1*x_1 - 1.0;
ret(1,0) = std::cos(x_2) + x_2*std::cos(x_1);
ret(1,1) = -x_1*std::sin(x_2) + std::cos(x_1);
//
return ret;
}
int main()
{
Eigen::VectorXd x = Eigen::VectorXd::Zero(2); // initial values (0,0)
bool success = optim::broyden_df(x, zeros_test_objfn_1, nullptr);
if (success) {
std::cout << "broyden_df: test_1 completed successfully." << "\n";
} else {
std::cout << "broyden_df: test_1 completed unsuccessfully." << "\n";
}
std::cout << "broyden_df: solution to test_1:\n" << x << std::endl;
//
x = Eigen::VectorXd::Zero(2);
success = optim::broyden_df(x, zeros_test_objfn_1, nullptr, zeros_test_jacob_1, nullptr);
if (success) {
std::cout << "broyden_df with jacobian: test_1 completed successfully." << "\n";
} else {
std::cout << "broyden_df with jacobian: test_1 completed unsuccessfully." << "\n";
}
std::cout << "broyden_df with jacobian: solution to test_1:\n" << x << std::endl;
//
return 0;
}
Example 2¶
Code to run this example is given below.
Armadillo (Click to show/hide)
#define OPTIM_ENABLE_ARMA_WRAPPERS
#include "optim.hpp"
inline
arma::vec
zeros_test_objfn_2(const arma::vec& vals_inp, void* opt_data)
{
double x_1 = vals_inp(0);
double x_2 = vals_inp(1);
//
arma::vec ret(2);
ret(0) = 2*x_1 - x_2 - std::exp(-x_1);
ret(1) = - x_1 + 2*x_2 - std::exp(-x_2);
//
return ret;
}
inline
arma::mat
zeros_test_jacob_2(const arma::vec& vals_inp, void* opt_data)
{
double x_1 = vals_inp(0);
double x_2 = vals_inp(1);
//
arma::mat ret(2,2);
ret(0,0) = 2 + std::exp(-x_1);
ret(0,1) = - 1.0;
ret(1,0) = - 1.0;
ret(1,1) = 2 + std::exp(-x_2);
//
return ret;
}
int main()
{
arma::vec x = arma::zeros(2,1); // initial values (0,0)
bool success = optim::broyden_df(x, zeros_test_objfn_2, nullptr);
if (success) {
std::cout << "broyden_df: test_2 completed successfully." << "\n";
} else {
std::cout << "broyden_df: test_2 completed unsuccessfully." << "\n";
}
arma::cout << "broyden_df: solution to test_2:\n" << x << arma::endl;
//
x = arma::zeros(2,1);
success = optim::broyden_df(x, zeros_test_objfn_2, nullptr, zeros_test_jacob_2, nullptr);
if (success) {
std::cout << "broyden_df with jacobian: test_2 completed successfully." << "\n";
} else {
std::cout << "broyden_df with jacobian: test_2 completed unsuccessfully." << "\n";
}
arma::cout << "broyden_df with jacobian: solution to test_2:\n" << x << arma::endl;
//
return 0;
}
Eigen (Click to show/hide)
#define OPTIM_ENABLE_EIGEN_WRAPPERS
#include "optim.hpp"
inline
Eigen::VectorXd
zeros_test_objfn_2(const Eigen::VectorXd& vals_inp, void* opt_data)
{
double x_1 = vals_inp(0);
double x_2 = vals_inp(1);
//
Eigen::VectorXd ret(2);
ret(0) = 2*x_1 - x_2 - std::exp(-x_1);
ret(1) = - x_1 + 2*x_2 - std::exp(-x_2);
//
return ret;
}
inline
Eigen::MatrixXd
zeros_test_jacob_2(const Eigen::VectorXd& vals_inp, void* opt_data)
{
double x_1 = vals_inp(0);
double x_2 = vals_inp(1);
//
Eigen::MatrixXd ret(2,2);
ret(0,0) = 2 + std::exp(-x_1);
ret(0,1) = - 1.0;
ret(1,0) = - 1.0;
ret(1,1) = 2 + std::exp(-x_2);
//
return ret;
}
int main()
{
Eigen::VectorXd x = Eigen::VectorXd::Zero(2); // initial values (0,0)
bool success = optim::broyden_df(x, zeros_test_objfn_2, nullptr);
if (success) {
std::cout << "broyden_df: test_2 completed successfully." << "\n";
} else {
std::cout << "broyden_df: test_2 completed unsuccessfully." << "\n";
}
std::cout << "broyden_df: solution to test_2:\n" << x << std::endl;
//
x = Eigen::VectorXd::Zero(2);
success = optim::broyden_df(x, zeros_test_objfn_2, nullptr, zeros_test_jacob_2, nullptr);
if (success) {
std::cout << "broyden_df with jacobian: test_2 completed successfully." << "\n";
} else {
std::cout << "broyden_df with jacobian: test_2 completed unsuccessfully." << "\n";
}
std::cout << "broyden_df with jacobian: solution to test_2:\n" << x << std::endl;
//
return 0;
}