gpmp::optim::QuasiNewton Class Reference

Class implementing Quasi-Newton optimization methods. More...

#include <quasi.hpp>

Public Member Functions
std::vector< double >	bhhh_optimize (const std::function< double(const std::vector< double > &)> &func, const std::vector< double > &initial_point, double tolerance, size_t max_iterations)
	Optimize a function using the Berndt–Hall–Hall–Hausman (BHHH) algorithm. More...
std::vector< double >	calculate_gradient (const std::function< double(const std::vector< double > &)> &func, const std::vector< double > &point, double epsilon)
	Calculate the gradient of a function at a given point. More...
std::vector< std::vector< double > >	calculate_bhhh_matrix (const std::vector< double > &gradient)
	Calculate the BHHH matrix from the gradient. More...
std::vector< double >	update_point (const std::vector< double > &current_point, const std::vector< double > &gradient, const std::vector< std::vector< double >> &bhhh_matrix)
	Update the current point using the BHHH matrix. More...
std::vector< double >	bfgs_optimize (const std::function< double(const std::vector< double > &)> &func, const std::vector< double > &initial_point, double tolerance, size_t max_iterations)
	Optimize a function using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. More...
std::vector< double >	calculate_search_direction (const std::vector< double > &gradient, const std::vector< std::vector< double >> &hessian_inverse)
	Calculate the search direction using the BFGS method. More...
double	line_search (const std::function< double(const std::vector< double > &)> &func, const std::vector< double > &current_point, const std::vector< double > &search_direction)
	Perform line search to find an appropriate step size. More...
std::vector< double >	update_point (const std::vector< double > &current_point, const std::vector< double > &search_direction, double step_size)
	Update the current point using the line search and step size. More...
std::vector< double >	calculate_gradient_difference (const std::vector< double > &next_point, const std::vector< double > &current_point, const std::vector< double > &gradient)
	Calculate the gradient difference between two points. More...
std::vector< std::vector< double > >	update_hessian_inverse (const std::vector< std::vector< double >> &hessian_inverse, const std::vector< double > &gradient_difference, const std::vector< double > &search_direction)
	Update the inverse of the Hessian matrix using the BFGS method. More...
double	dot_product (const std::vector< double > &a, const std::vector< double > &b)
	Calculate the dot product of two vectors. More...
std::vector< double >	vector_subtraction (const std::vector< double > &a, const std::vector< double > &b) const
	Subtract two vectors element-wise. More...
std::tuple< std::vector< double >, double >	lbfgs_optimize (const std::function< double(const std::vector< double > &)> &f, const std::vector< double > &initial_point, double tolerance=1e-4, size_t max_iterations=100, size_t memory_size=5)
	Optimize a function using the Limited-Memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm. More...

Detailed Description

Class implementing Quasi-Newton optimization methods.

Definition at line 48 of file quasi.hpp.

Member Function Documentation

bfgs_optimize()

std::vector< double > gpmp::optim::QuasiNewton::bfgs_optimize	(	const std::function< double(const std::vector< double > &)> &	func,
		const std::vector< double > &	initial_point,
		double	tolerance,
		size_t	max_iterations
	)

Optimize a function using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.

Parameters

func	The objective function to minimize
initial_point	The initial guess for the optimal parameters
tolerance	The tolerance for stopping criterion
max_iterations	The maximum number of iterations

Returns

The vector of optimal parameters

Definition at line 145 of file quasi.cpp.

                           {
    std::vector<double> current_point = initial_point;
    size_t n = initial_point.size();
 
    // Initialize Hessian approximation as the identity matrix
    std::vector<std::vector<double>> hessian_inverse(
        n,
        std::vector<double>(n, 0.0));
    for (size_t i = 0; i < n; ++i) {
        hessian_inverse[i][i] = 1.0;
    }
 
    for (size_t iteration = 0; iteration < max_iterations; ++iteration) {
        // Calculate the gradient
        std::vector<double> gradient =
            calculate_gradient(func, current_point, 1e-6);
 
        // Check convergence
        double gradient_norm = 0.0;
        for (size_t i = 0; i < n; ++i) {
            gradient_norm += gradient[i] * gradient[i];
        }
        gradient_norm = std::sqrt(gradient_norm);
 
        if (gradient_norm < tolerance) {
            std::cout << "Converged after " << iteration << " iterations."
                      << std::endl;
            return current_point;
        }
 
        // Calculate search direction
        std::vector<double> search_direction =
            calculate_search_direction(gradient, hessian_inverse);
 
        // Perform line search to find the step size
        double step_size = line_search(func, current_point, search_direction);
 
        // Update the current point
        std::vector<double> next_point =
            update_point(current_point, search_direction, step_size);
 
        // Update the Hessian approximation
        std::vector<double> gradient_difference =
            calculate_gradient_difference(next_point, current_point, gradient);
        hessian_inverse = update_hessian_inverse(hessian_inverse,
                                                 gradient_difference,
                                                 search_direction);
 
        // Move to the next iteration
        current_point = next_point;
    }
 
    std::cout << "Reached maximum iterations without convergence." << std::endl;
    return current_point;
}

bhhh_optimize()

std::vector< double > gpmp::optim::QuasiNewton::bhhh_optimize	(	const std::function< double(const std::vector< double > &)> &	func,
		const std::vector< double > &	initial_point,
		double	tolerance,
		size_t	max_iterations
	)

Optimize a function using the Berndt–Hall–Hall–Hausman (BHHH) algorithm.

Parameters

func	The objective function to minimize
initial_point	The initial guess for the optimal parameters
tolerance	The tolerance for stopping criterion
max_iterations	The maximum number of iterations

Returns

The vector of optimal parameters

Definition at line 58 of file quasi.cpp.

                           {
    std::vector<double> current_point = initial_point;
    size_t n = initial_point.size();
 
    for (size_t iteration = 0; iteration < max_iterations; ++iteration) {
        // Calculate the gradient
        std::vector<double> gradient =
            calculate_gradient(func, current_point, 1e-6);
 
        // Check convergence
        double gradient_norm = 0.0;
        for (size_t i = 0; i < n; ++i) {
            gradient_norm += gradient[i] * gradient[i];
        }
        gradient_norm = std::sqrt(gradient_norm);
 
        if (gradient_norm < tolerance) {
            std::cout << "Converged after " << iteration << " iterations."
                      << std::endl;
            return current_point;
        }
 
        // Calculate the BHHH matrix
        std::vector<std::vector<double>> bhhh_matrix =
            calculate_bhhh_matrix(gradient);
 
        // Update the current point
        current_point = update_point(current_point, gradient, bhhh_matrix);
    }
 
    std::cout << "Reached maximum iterations without convergence." << std::endl;
    return current_point;
}

calculate_bhhh_matrix()

std::vector< std::vector< double > > gpmp::optim::QuasiNewton::calculate_bhhh_matrix

(

const std::vector< double > &

gradient

)

Calculate the BHHH matrix from the gradient.

Parameters

gradient

The gradient vector

Returns

The BHHH matrix

Definition at line 117 of file quasi.cpp.

                                       {
    size_t n = gradient.size();
    std::vector<std::vector<double>> bhhh_matrix(n, std::vector<double>(n));
 
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            bhhh_matrix[i][j] = gradient[i] * gradient[j];
        }
    }
 
    return bhhh_matrix;
}

calculate_gradient()

std::vector< double > gpmp::optim::QuasiNewton::calculate_gradient	(	const std::function< double(const std::vector< double > &)> &	func,
		const std::vector< double > &	point,
		double	epsilon
	)

Calculate the gradient of a function at a given point.

Parameters

func	The objective function
point	The point at which to calculate the gradient
epsilon	The perturbation for finite differences

Returns

The gradient vector

Definition at line 96 of file quasi.cpp.

                    {
    size_t n = point.size();
    std::vector<double> gradient(n);
 
    for (size_t i = 0; i < n; ++i) {
        std::vector<double> perturbed_point = point;
        perturbed_point[i] += epsilon;
 
        double perturbed_value = func(perturbed_point);
        double original_value = func(point);
 
        gradient[i] = (perturbed_value - original_value) / epsilon;
    }
 
    return gradient;
}

calculate_gradient_difference()

std::vector< double > gpmp::optim::QuasiNewton::calculate_gradient_difference	(	const std::vector< double > &	next_point,
		const std::vector< double > &	current_point,
		const std::vector< double > &	gradient
	)

Calculate the gradient difference between two points.

Parameters

next_point	The next point
current_point	The current point
gradient	The gradient vector

Returns

The gradient difference vector

Definition at line 271 of file quasi.cpp.

                                       {
    size_t n = next_point.size();
    std::vector<double> gradient_difference(n);
 
    for (size_t i = 0; i < n; ++i) {
        gradient_difference[i] =
            gradient[i] * (next_point[i] - current_point[i]);
    }
 
    return gradient_difference;
}

calculate_search_direction()

std::vector< double > gpmp::optim::QuasiNewton::calculate_search_direction	(	const std::vector< double > &	gradient,
		const std::vector< std::vector< double >> &	hessian_inverse
	)

Calculate the search direction using the BFGS method.

Parameters

gradient	The gradient vector
hessian_inverse	The inverse of the Hessian matrix

Returns

The search direction vector

Definition at line 205 of file quasi.cpp.

                                                         {
    size_t n = gradient.size();
    std::vector<double> search_direction(n);
 
    for (size_t i = 0; i < n; ++i) {
        search_direction[i] = 0.0;
        for (size_t j = 0; j < n; ++j) {
            search_direction[i] -= hessian_inverse[i][j] * gradient[j];
        }
    }
 
    return search_direction;
}

dot_product()

double gpmp::optim::QuasiNewton::dot_product	(	const std::vector< double > &	a,
		const std::vector< double > &	b
	)

Calculate the dot product of two vectors.

Parameters

a	The first vector
b	The second vector

Returns

The dot product value

Definition at line 311 of file quasi.cpp.

                                                                         {
    // Ensure vectors have the same size
    if (a.size() != b.size()) {
        throw std::invalid_argument(
            "Vectors must have the same size for dot product.");
    }
 
    double result = 0.0;
    for (size_t i = 0; i < a.size(); ++i) {
        result += a[i] * b[i];
    }
    return result;
}

lbfgs_optimize()

std::tuple< std::vector< double >, double > gpmp::optim::QuasiNewton::lbfgs_optimize	(	const std::function< double(const std::vector< double > &)> &	f,
		const std::vector< double > &	initial_point,
		double	tolerance = 1e-4,
		size_t	max_iterations = 100,
		size_t	memory_size = 5
	)

Optimize a function using the Limited-Memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm.

Parameters

func	The objective function to optimize
initial_point	The initial guess for the optimal point
tolerance	The tolerance for stopping criterion
max_iterations	The maximum number of iterations
memory_size	The size of the limited-memory history (s and y vectors)

Returns

The optimized point that minimizes the given objective function

Definition at line 327 of file quasi.cpp.

                        {
 
    const double eps = 1e-8;
 
    size_t n = initial_point.size();
    std::vector<double> x = initial_point;
    std::vector<double> g(n); // Gradient vector
    std::vector<std::vector<double>> s(memory_size,
                                       std::vector<double>(n)); // s vectors
    std::vector<std::vector<double>> y(memory_size,
                                       std::vector<double>(n)); // y vectors
    std::vector<double> rho(memory_size);                       // rho values
 
    // Evaluate the objective function and gradient at initial_point
    double fx = f(x);
    // Calculate gradient at initial_point
    // Gradient calculation logic to be implemented
    // Assign gradient to 'g'
 
    for (size_t iter = 0; iter < max_iterations; ++iter) {
        // Check for convergence
        double norm_grad = 0.0;
        for (size_t i = 0; i < n; ++i) {
            norm_grad += g[i] * g[i];
        }
        norm_grad = sqrt(norm_grad);
        if (norm_grad < tolerance) {
            break;
        }
 
        // Compute search direction (use initial guess)
        std::vector<double> d = g;
 
        // L-BFGS two-loop recursion
        size_t start = std::min(iter, memory_size);
        // for (size_t i = start - 1; i >= 0; --i) {
        for (size_t i = start; i > 0; --i) {
 
            rho[i] = 1.0 /
                     inner_product(s[i].begin(), s[i].end(), y[i].begin(), 0.0);
            double alpha =
                rho[i] *
                inner_product(s[i].begin(), s[i].end(), d.begin(), 0.0);
            for (size_t j = 0; j < n; ++j) {
                d[j] -= alpha * y[i][j];
            }
        }
 
        // Perform scaling
        for (size_t i = 0; i < n; ++i) {
            d[i] *= rho[i];
        }
 
        // Compute gradient of the objective function along the search direction
        // Gradient calculation logic to be implemented
        // Assign gradient to 'dg'
        double dg = inner_product(d.begin(), d.end(), g.begin(), 0.0);
 
        // Limit curvature
        if (dg > 0) {
            break;
        }
 
        // Line search
        double step_size = 1.0;
        std::vector<double> x_new = x;
        for (size_t i = 0; i < n; ++i) {
            x_new[i] += step_size * d[i];
        }
 
        double fx_new = f(x_new);
        if (fx_new < fx + eps * step_size * dg) {
            // Update x
            x = x_new;
            fx = fx_new;
 
            // Evaluate gradient at new point
            // Gradient calculation logic to be implemented
            // Assign gradient to 'g'
 
            // Update s and y
            for (size_t i = 0; i < n; ++i) {
                s[iter % memory_size][i] = x_new[i] - x[i];
                y[iter % memory_size][i] = g[i] - d[i];
            }
        }
    }
 
    return std::make_tuple(x, fx);
}

line_search()

double gpmp::optim::QuasiNewton::line_search	(	const std::function< double(const std::vector< double > &)> &	func,
		const std::vector< double > &	current_point,
		const std::vector< double > &	search_direction
	)

Perform line search to find an appropriate step size.

Parameters

func	The objective function
current_point	The current point
search_direction	The search direction vector

Returns

The optimal step size

Definition at line 221 of file quasi.cpp.

                                               {
    const double alpha = 0.001;      // Step size multiplier
    const double beta = 0.5;         // Factor for reducing the step size
    const int maxIterations = 100;   // Maximum number of iterations
    const double minStepSize = 1e-6; // Minimum step size
 
    double step_size = 1.0; // Initial step size
    std::vector<double> updated_point = current_point;
 
    // Evaluate the objective function at the current point
    double f_current = func(current_point);
 
    // Calculate the directional derivative (gradient dot search_direction)
    double directional_derivative =
        dot_product(calculate_gradient(func, current_point, 1e-6),
                    search_direction);
 
    int iteration = 0;
    while (step_size > minStepSize && iteration < maxIterations) {
        updated_point =
            update_point(current_point, search_direction, step_size);
        double f_updated = func(updated_point);
        if (f_updated <=
            f_current + alpha * step_size * directional_derivative) {
            break; // Stop if Armijo condition is satisfied
        }
        step_size *= beta; // Reduce the step size
        ++iteration;
    }
 
    return step_size;
}

References gpmp::linalg::dot_product().

update_hessian_inverse()

std::vector< std::vector< double > > gpmp::optim::QuasiNewton::update_hessian_inverse	(	const std::vector< std::vector< double >> &	hessian_inverse,
		const std::vector< double > &	gradient_difference,
		const std::vector< double > &	search_direction
	)

Update the inverse of the Hessian matrix using the BFGS method.

Parameters

hessian_inverse	The current inverse of the Hessian matrix
gradient_difference	The gradient difference vector
search_direction	The search direction vector

Returns

The updated inverse of the Hessian matrix

Definition at line 287 of file quasi.cpp.

                                               {
 
    size_t n = hessian_inverse.size();
    std::vector<std::vector<double>> updated_hessian_inverse(
        n,
        std::vector<double>(n));
 
    // Update Hessian using BFGS update formula
    double rho = dot_product(gradient_difference, search_direction);
 
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            updated_hessian_inverse[i][j] =
                hessian_inverse[i][j] +
                rho * gradient_difference[i] * gradient_difference[j];
        }
    }
 
    return updated_hessian_inverse;
}

References gpmp::linalg::dot_product().

update_point() [1/2]

std::vector< double > gpmp::optim::QuasiNewton::update_point	(	const std::vector< double > &	current_point,
		const std::vector< double > &	gradient,
		const std::vector< std::vector< double >> &	bhhh_matrix
	)

Update the current point using the BHHH matrix.

Parameters

current_point	The current point
gradient	The gradient vector
bhhh_matrix	The BHHH matrix

Returns

The updated point

Definition at line 131 of file quasi.cpp.

                                                     {
    size_t n = current_point.size();
    std::vector<double> updated_point(n);
 
    for (size_t i = 0; i < n; ++i) {
        updated_point[i] = current_point[i] - gradient[i] / bhhh_matrix[i][i];
    }
 
    return updated_point;
}

update_point() [2/2]

std::vector< double > gpmp::optim::QuasiNewton::update_point	(	const std::vector< double > &	current_point,
		const std::vector< double > &	search_direction,
		double	step_size
	)

Update the current point using the line search and step size.

Parameters

current_point	The current point
search_direction	The search direction vector
step_size	The step size from the line search

Returns

The updated point

Definition at line 257 of file quasi.cpp.

                      {
    size_t n = current_point.size();
    std::vector<double> updated_point(n);
 
    for (size_t i = 0; i < n; ++i) {
        updated_point[i] = current_point[i] + step_size * search_direction[i];
    }
 
    return updated_point;
}

vector_subtraction()

std::vector< double > gpmp::optim::QuasiNewton::vector_subtraction	(	const std::vector< double > &	a,
		const std::vector< double > &	b
	)		const

Subtract two vectors element-wise.

Parameters

a	The first vector
b	The second vector

Returns

The result of subtracting each element of vector b from vector a

Definition at line 40 of file quasi.cpp.

                                      {
    if (a.size() != b.size()) {
        throw std::invalid_argument(
            "Error: Vector dimensions do not match for subtraction.");
    }
 
    std::vector<double> result;
    result.reserve(a.size());
 
    for (size_t i = 0; i < a.size(); ++i) {
        result.push_back(a[i] - b[i]);
    }
 
    return result;
}

---

The documentation for this class was generated from the following files:

• include/openGPMP/optim/quasi.hpp
• modules/optim/quasi.cpp