Optimization Demo

The sangi optimization module handles a wide range of problems, from unconstrained nonlinear minimization to linear programming. This page walks through three demos that illustrate representative use cases.

Demo 1: Rosenbrock minimization (BFGS)

The Rosenbrock function $f(x, y) = (1 - x)^2 + 100(y - x^2)^2$ is a standard benchmark for optimization algorithms. Its global minimum is at $(x, y) = (1, 1)$ with $f = 0$, but converging along the narrow banana-shaped valley is challenging. Here we minimize starting from $(-1, -1)$ with BFGS.

=== Rosenbrock minimization (BFGS) === initial point: (-1, -1) iterations: 34 minimum at: (1.000000, 1.000000) f(x,y) = 2.53e-20 converged: true

Demo 1 source code

// Rosenbrock minimization with BFGS
#include <math/optimization/optimization_nd.hpp>
#include <iostream>
#include <iomanip>
using namespace sangi;

void demo_rosenbrock() {
    std::cout << "=== Rosenbrock minimization (BFGS) ===\n";

    // Objective: f(x,y) = (1-x)^2 + 100(y-x^2)^2
    auto rosenbrock = [](const std::vector<double>& x) {
        double a = 1.0 - x[0];
        double b = x[1] - x[0] * x[0];
        return a * a + 100.0 * b * b;
    };

    std::vector<double> x0 = {-1.0, -1.0};
    std::cout << "initial point: (" << x0[0] << ", " << x0[1] << ")\n";

    auto result = bfgs_minimize(rosenbrock, x0);

    std::cout << "iterations: " << result.iterations << "\n";
    std::cout << std::setprecision(6) << std::fixed;
    std::cout << "minimum at: (" << result.x[0] << ", " << result.x[1] << ")\n";
    std::cout << std::scientific << std::setprecision(2);
    std::cout << "f(x,y)   = " << result.value << "\n";
    std::cout << "converged: " << (result.converged ? "true" : "false") << "\n";
}

BFGS is a quasi-Newton method that incrementally updates an approximation of the inverse Hessian to determine the search direction. The gradient is computed automatically via numerical differentiation, so the user is not required to supply derivatives. 34 iterations bring $f$ down to $\approx 10^{-20}$, which is essentially the exact solution.

Demo 2: Curve fitting (curveFit)

Generate noisy Gaussian data of the form $y = a \exp\!\bigl(-(x - b)^2 / (2c^2)\bigr)$ and estimate the parameters $a$, $b$, and $c$ with curveFit. We compare the recovered parameters against the ground truth to assess accuracy.

=== Gaussian curve fitting === true parameters: a=3.000, b=2.000, c=0.500 fitted parameters: a=2.987, b=2.003, c=0.508 residual norm: 0.142

Demo 2 source code

// Gaussian curve fitting
#include <math/optimization/optimization_nd.hpp>
#include <iostream>
#include <iomanip>
#include <random>
using namespace sangi;

void demo_curve_fit() {
    std::cout << "\n=== Gaussian curve fitting ===\n";

    // Ground-truth parameters
    double a_true = 3.0, b_true = 2.0, c_true = 0.5;

    // Generate noisy data
    std::mt19937 rng(42);
    std::normal_distribution<double> noise(0.0, 0.05);

    std::vector<double> xs, ys;
    for (int i = 0; i <= 40; ++i) {
        double x = i * 0.1;  // 0.0 to 4.0
        double y = a_true * std::exp(-(x - b_true) * (x - b_true)
                                      / (2.0 * c_true * c_true))
                   + noise(rng);
        xs.push_back(x);
        ys.push_back(y);
    }

    // Model function: params = {a, b, c}
    auto gaussian = [](double x, const std::vector<double>& p) {
        return p[0] * std::exp(-(x - p[1]) * (x - p[1])
                                / (2.0 * p[2] * p[2]));
    };

    // Initial guess
    std::vector<double> p0 = {1.0, 1.0, 1.0};
    auto result = curveFit(gaussian, xs, ys, p0);

    std::cout << std::fixed << std::setprecision(3);
    std::cout << "true parameters:    a=" << a_true
              << ", b=" << b_true << ", c=" << c_true << "\n";
    std::cout << "fitted parameters:  a=" << result.params[0]
              << ", b=" << result.params[1]
              << ", c=" << result.params[2] << "\n";
    std::cout << "residual norm: " << result.residualNorm << "\n";
}

Internally, curveFit runs nonlinear least squares with the Levenberg-Marquardt algorithm. With noise of standard deviation 0.05, the recovered parameters match the ground truth to high accuracy. The initial guess $(1, 1, 1)$ is far from the true values, yet the fit converges without trouble.

Demo 3: Linear programming (Simplex)

Solve a linear program with the simplex method.

$$\text{maximize} \quad 3x + 5y$$ $$\text{subject to} \quad x + y \le 4, \quad x + 3y \le 6, \quad x, y \ge 0$$

The simplex method visits the vertices of the feasible region defined by the constraints to find the optimum.

=== Linear programming (Simplex) === maximize: 3x + 5y subject to: x + y <= 4 x + 3y <= 6 x, y >= 0 optimal value: 16.000000 optimal solution: x=3.000000, y=1.000000 status: optimal

Demo 3 source code

// Linear programming with the simplex method
#include <math/optimization/optimization_nd.hpp>
#include <iostream>
#include <iomanip>
using namespace sangi;

void demo_simplex() {
    std::cout << "\n=== Linear programming (Simplex) ===\n";
    std::cout << "maximize: 3x + 5y\n";
    std::cout << "subject to:\n";
    std::cout << "  x +  y <= 4\n";
    std::cout << "  x + 3y <= 6\n";
    std::cout << "  x, y >= 0\n\n";

    // Objective coefficients (maximization)
    std::vector<double> c = {3.0, 5.0};

    // Inequality constraints Ax <= b
    std::vector<std::vector<double>> A = {
        {1.0, 1.0},   // x + y <= 4
        {1.0, 3.0}    // x + 3y <= 6
    };
    std::vector<double> b = {4.0, 6.0};

    auto result = simplexMaximize(c, A, b);

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "optimal value: " << result.optimalValue << "\n";
    std::cout << "optimal solution: x=" << result.solution[0]
              << ", y=" << result.solution[1] << "\n";
    std::cout << "status: " << (result.feasible ? "optimal" : "infeasible") << "\n";
}

The optimal solution is $(x, y) = (3, 1)$ with an objective value of $3 \times 3 + 5 \times 1 = 16$. This is the intersection of the two active constraints $x + y = 4$ and $x + 3y = 6$. simplexMaximize handles conversion to standard form and the pivoting operations internally — the user need not deal with slack variables explicitly.

Source code and how to build

example_optimization.cpp (full source)

// example_optimization.cpp
// Optimization demo: Rosenbrock BFGS, Gaussian curve fitting, simplex LP

#define _USE_MATH_DEFINES
#include <cmath>
#include <math/optimization/optimization_nd_impl.hpp>
#include <math/optimization/linear_programming_impl.hpp>
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>
#include <iostream>
#include <iomanip>
#include <random>
#include <functional>

using namespace sangi;

int main() {
    // --- Demo 1: Rosenbrock BFGS ---
    {
        std::cout << "=== Rosenbrock minimization (BFGS) ===" << std::endl;

        auto rosenbrock = [](const Vector<double>& x) -> double {
            double a = 1.0 - x[0];
            double b = x[1] - x[0] * x[0];
            return a * a + 100.0 * b * b;
        };
        auto rosenbrock_grad = [](const Vector<double>& x) -> Vector<double> {
            Vector<double> g(2);
            g[0] = -2.0 * (1.0 - x[0]) - 400.0 * x[0] * (x[1] - x[0] * x[0]);
            g[1] = 200.0 * (x[1] - x[0] * x[0]);
            return g;
        };

        Vector<double> x0(2);
        x0[0] = -1.0; x0[1] = -1.0;
        std::cout << "initial point: (" << x0[0] << ", " << x0[1] << ")" << std::endl;

        OptimizationOptions<double> opts;
        opts.tolerance = 1e-8;
        opts.max_iterations = 10000;

        auto result = bfgs_minimize<double, Vector<double>, Matrix<double>>(
            rosenbrock, rosenbrock_grad, x0, opts);

        std::cout << "iterations: " << result.iterations << std::endl;
        std::cout << std::fixed << std::setprecision(6);
        std::cout << "minimum at: (" << result.point[0] << ", " << result.point[1] << ")" << std::endl;
        std::cout << std::scientific << std::setprecision(2);
        std::cout << "f(x,y)    = " << result.value << std::endl;
        std::cout << "converged: " << (result.success ? "true" : "false") << std::endl;
    }

    // --- Demo 2: Gaussian curve fit ---
    {
        std::cout << "\n=== Gaussian curve fitting ===" << std::endl;

        double a_true = 3.0, b_true = 2.0, c_true = 0.5;

        std::mt19937 rng(42);
        std::normal_distribution<double> noise(0.0, 0.05);

        std::vector<double> xs, ys;
        for (int i = 0; i <= 40; ++i) {
            double x = i * 0.1;
            double y = a_true * std::exp(-(x - b_true) * (x - b_true)
                                          / (2.0 * c_true * c_true))
                       + noise(rng);
            xs.push_back(x);
            ys.push_back(y);
        }

        auto gaussian = [](double x, const Vector<double>& p) -> double {
            return p[0] * std::exp(-(x - p[1]) * (x - p[1])
                                    / (2.0 * p[2] * p[2]));
        };

        Vector<double> p0(3);
        p0[0] = 1.0; p0[1] = 1.0; p0[2] = 1.0;
        auto result = curveFit<double>(gaussian, xs, ys, p0);

        std::cout << std::fixed << std::setprecision(3);
        std::cout << "true parameters:    a=" << a_true
                  << ", b=" << b_true << ", c=" << c_true << std::endl;
        std::cout << "fitted parameters:  a=" << result.parameters[0]
                  << ", b=" << result.parameters[1]
                  << ", c=" << result.parameters[2] << std::endl;
        std::cout << "chi^2: " << result.chiSquared << std::endl;
    }

    // --- Demo 3: Linear programming (Simplex) ---
    {
        std::cout << "\n=== Linear programming (Simplex) ===" << std::endl;
        std::cout << "maximize: 3x + 5y" << std::endl;
        std::cout << "subject to:" << std::endl;
        std::cout << "  x +  y <= 4" << std::endl;
        std::cout << "  x + 3y <= 6" << std::endl;
        std::cout << "  x, y >= 0" << std::endl;

        LinearProgram<double> lp;
        lp.objective = {3.0, 5.0};
        lp.constraints.push_back(LinearConstraint<double>({1.0, 1.0}, ConstraintType::LessEqual, 4.0));
        lp.constraints.push_back(LinearConstraint<double>({1.0, 3.0}, ConstraintType::LessEqual, 6.0));

        auto result = simplexMaximize(lp);

        std::cout << std::fixed << std::setprecision(3);
        std::cout << "objective: " << result.objectiveValue << std::endl;
        std::cout << "x = " << result.solution[0] << ", y = " << result.solution[1] << std::endl;
    }

    return 0;
}

For full API details, see the Optimization API reference.

Build and run

Because of template instantiation, the user must include optimization_nd_impl.hpp and linear_programming_impl.hpp (header-only).

cd sangi
mkdir build && cd build
cmake .. -G "Visual Studio 17 2022" -A x64
cmake --build . --config Release --target example-optimization
examples\Release\example-optimization.exe

CMakeLists.txt snippet

target_include_directories(demo PRIVATE <sangi>/include)