# Minimizing the Rosenbrock function

MUQ contsins methods for solving both constrained and unconstrained optimization problems. Here we demonstrate some of the unconstrained optimization capabilities in MUQ by minimizing the Rosenbrock function.

Let $x\in\mathbb{R}^2$ denote the two dimensional decision variable and let $f(x)$ denote the objective function, which is given by the well known Rosenbrock function

$$f(x) = \left(1 - x_1\right)^2 + 100(x_2-x_1^2)^2$$

This function has a global minimum at $x=[x_1,x_2] = [1,1]$.

## Optimization in MUQ

In MUQ, optimization problems are defined as children of the abstract muq::Optimization::OptProbBase class. Thus, to define the Rosenbrock problem, we need inherit from this class and implement the objective function. An instance of this class can then be passed to an optimization algorithm (defined through children of the muq::Optimization::OptAlgBase class).

## Include the necessary header files

#include <Eigen/Dense>

#include "MUQ/Optimization/Problems/OptProbBase.h"
#include "MUQ/Optimization/Algorithms/OptAlgBase.h"



## Defining the objective

Here we define a class, called RoseFunc that inherits from the optimization base class muq::Optimization::OptProbBase.

### Constructor

The constructor of muq::Optimization::OptProbBase accepts the number of decision variables (2 in this case).


class RoseFunc : public muq::Optimization::OptProbBase {
public:

RoseFunc() : muq::Optimization::OptProbBase(2) {}


### Objective function, i.e. eval

All user-defined objective functions must define the objective function by implementing the eval function. Here, we create the objective function for the Rosenbrock function.

  virtual double eval(const Eigen::VectorXd& xc) override
{
return pow(1 - xc[0], 2) + 100 * pow(xc[1] - xc[0] * xc[0], 2);
}



Some optimizers can take advantage of gradient and Hessian information. This information is provied by implementing the grad and applyInvHess functions. Note that both of these functions are optional and will be replaced by finite difference approximations if they are not provided.

$$\nabla f(x) = \left[ \begin{array}{c} -2(1-x_1) - 400x_1\left(x_2-x_1^2\right) \\ 200\left(x_2-x_1^2\right) \end{array}\right]^T .$$

The grad function computes updates the gradient, which is passed by reference, and returns the objective function value. In some cases, computing these quantities at the same time can be more computationally efficient.


{

gradient[0] = -400 * (xc[1] - xc[0] * xc[0]) * xc[0] - 2 * (1 - xc[0]);
gradient[1] = 200 * (xc[1] - xc[0] * xc[0]);

return eval(xc);
}



The Hessian matrix of the Rosenbrock function is given by

$$H(x) = \left[ \begin{array}{cc} 2 - 400x_2 + 1200x_1^2 & -400x_1 \\ -400x_1 & 200\end{array}\right].$$

Notice that the applyInvHess function applies the inverse Hessian to a matrix and does not return the actual Hessian matrix. This allows for flexibility in how the inverse action is computed (e.g., adjoint methods with iterative solvers). However, for some users it may be more convenient to simply return the Hessian or inverse Hessian. In these situations, users can instead overload the getHess and getInvHess functions. Note that the getHess or getInvHess functions will not be used if the applyInvHess function is implemented.

  virtual Eigen::VectorXd applyInvHess(const Eigen::VectorXd& xc, const Eigen::VectorXd& vecIn)
{
Eigen::Matrix<double, 2, 2> Hess = Eigen::Matrix<double, 2, 2>::Zero(2, 2);

Hess(0, 0) = 1200 * pow(xc[0], 2.0) - 400 * xc[1] + 2;
Hess(0, 1) = -400 * xc[0];
Hess(1, 0) = -400 * xc[0];
Hess(1, 1) = 200;

return Hess.lu().solve(vecIn);
}

}; // End of RoseFunc class



## Solving the problem

Now that we've defined our problem in the RoseFunc class, we can set up an optimizer and minimize the objective function.

### Create an instance of the objective

We begin the main function here by creating an instance of the Rosenbrock optimization function defined above.



int main()
{

// create an instance of the optimization problem
auto prob = std::make_shared<RoseFunc>();



### Set up the optimizer

Now we create an optimizer. MUQ uses the factory method muq::Optimization::OptAlgBase::Create to construct an optimizer based on the problem and other optimization parameters defined in a boost::property_tree::ptree.

The optimization algorithm, set by the Opt.Method parameter, can be one of

• SD_Line, which will create a steepest descent solver
• BFGS_Line, which will create a BFGS solver
• Newton, which will create a Newton solver

If MUQ was compiled with NLOPT. It is also possible to set Opt.Method to NLOPT. The specific optimization algorithm is then specified by the Opt.NLOPT.Method key. See the parameter list here for more options.


// set the initial condition
Eigen::VectorXd x0(2);
x0 << -1, 3;

boost::property_tree::ptree params;

// set some of the optimization parameters
params.put("Opt.MaxIts", 10000);
params.put("Opt.ftol", 1e-8);
params.put("Opt.xtol", 1e-8);
params.put("Opt.gtol", 1e-8);
params.put("Opt.LineSearch.LineIts", 100);
params.put("Opt.StepLength", 1);
params.put("Opt.verbosity", 3);

//params.put("Opt.Method", "SD_Line");    // Use the steepest descent algorithm
//params.put("Opt.Method", "BFGS_Line");  // Use the BFGS algorithm
params.put("Opt.Method", "Newton");     // Use the Newton algorithm

// set up the solver
auto Solver = muq::Optimization::OptAlgBase::Create(prob, params);


### Run the optimizer

The OptAlgBase::solve function runs the optimization algorithm and, upon successful completion, returns the optimal point.

The termination status of the solver can be checked with the OptAlgBase::GetStatus function. Note that MUQ's termination codes are identical to NLOPT: positive numbers indicate successful termination and negative numbers indicate that an error occured.


// solve the optimization problem
Eigen::VectorXd xOpt = Solver->solve(x0);

// Get the termination status
int optStat = Solver->GetStatus();

std::cout << "Optimal solution = " << xOpt.transpose() << std::endl;

return 0;
} // End of "int main()"


## Build the executable

cd build; cmake ../ > BuildLog.txt; make; cd ../

Scanning dependencies of target RosenbrockOpt
[100%] Building CXX object CMakeFiles/RosenbrockOpt.dir/RosenbrockOpt.cpp.o
[100%] Built target RosenbrockOpt



## Run the executable

build/RosenbrockOpt

Using optimization method: Newton
Iteration: 1   fval: 404
Iteration: 2   fval: 4.02008
Iteration: 3   fval: 3.46927
Iteration: 4   fval: 2.66982
Iteration: 5   fval: 2.25626
Iteration: 6   fval: 1.71473
Iteration: 7   fval: 1.40064
Iteration: 8   fval: 0.943349
Iteration: 9   fval: 0.754829
Iteration: 10   fval: 0.488025
Iteration: 11   fval: 0.362813
Iteration: 12   fval: 0.223237
Iteration: 13   fval: 0.133336
Iteration: 14   fval: 0.0661269
Iteration: 15   fval: 0.0353271
Iteration: 16   fval: 0.0098764
Iteration: 17   fval: 0.00387419
Iteration: 18   fval: 0.000141491
Iteration: 19   fval: 1.87429e-06
Iteration: 20   fval: 1.63594e-19
Terminating with status: 3
Optimal solution = 1 1



## Completed code:

### RosenbrockOpt.cpp

#include <Eigen/Dense>

#include "MUQ/Optimization/Problems/OptProbBase.h"
#include "MUQ/Optimization/Algorithms/OptAlgBase.h"

class RoseFunc : public muq::Optimization::OptProbBase {
public:

RoseFunc() : muq::Optimization::OptProbBase(2) {}

virtual double eval(const Eigen::VectorXd& xc) override
{
return pow(1 - xc[0], 2) + 100 * pow(xc[1] - xc[0] * xc[0], 2);
}

{

gradient[0] = -400 * (xc[1] - xc[0] * xc[0]) * xc[0] - 2 * (1 - xc[0]);
gradient[1] = 200 * (xc[1] - xc[0] * xc[0]);

return eval(xc);
}

virtual Eigen::VectorXd applyInvHess(const Eigen::VectorXd& xc, const Eigen::VectorXd& vecIn)
{
Eigen::Matrix<double, 2, 2> Hess = Eigen::Matrix<double, 2, 2>::Zero(2, 2);

Hess(0, 0) = 1200 * pow(xc[0], 2.0) - 400 * xc[1] + 2;
Hess(0, 1) = -400 * xc[0];
Hess(1, 0) = -400 * xc[0];
Hess(1, 1) = 200;

return Hess.lu().solve(vecIn);
}

}; // End of RoseFunc class

int main()
{

// create an instance of the optimization problem
auto prob = std::make_shared<RoseFunc>();

// set the initial condition
Eigen::VectorXd x0(2);
x0 << -1, 3;

boost::property_tree::ptree params;

// set some of the optimization parameters
params.put("Opt.MaxIts", 10000);
params.put("Opt.ftol", 1e-8);
params.put("Opt.xtol", 1e-8);
params.put("Opt.gtol", 1e-8);
params.put("Opt.LineSearch.LineIts", 100);
params.put("Opt.StepLength", 1);
params.put("Opt.verbosity", 3);

//params.put("Opt.Method", "SD_Line");    // Use the steepest descent algorithm
//params.put("Opt.Method", "BFGS_Line");  // Use the BFGS algorithm
params.put("Opt.Method", "Newton");     // Use the Newton algorithm

// set up the solver
auto Solver = muq::Optimization::OptAlgBase::Create(prob, params);

// solve the optimization problem
Eigen::VectorXd xOpt = Solver->solve(x0);

// Get the termination status
int optStat = Solver->GetStatus();

std::cout << "Optimal solution = " << xOpt.transpose() << std::endl;

return 0;
} // End of "int main()"


Test Status
• ##### Develop Branch
OS Compiler
OSX Clang
Ubuntu Clang
Ubuntu g++ 4.7
Ubuntu g++ 4.8
Ubuntu g++ 4.9
Acknowledgments
This material is based upon work supported by the National Science Foundation under Grant No. 1550487.

This work was supported by the DOE Office of Science through the QUEST SciDAC Institute.

Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.