### Install Eigen on Ubuntu

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Install the Eigen library on Ubuntu using apt.

```bash
sudo apt install libeigen3-dev
```

--------------------------------

### Install Documentation Dependencies

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/UpdatingDocs.md

Install the required MkDocs and math support packages via pip.

```bash
pip3 install mkdocs
pip3 install python-markdown-math
```

--------------------------------

### Build Mathtoolbox C++ Library

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Build the mathtoolbox C++ library using CMake. Ensure Eigen is installed. The default build does not include examples.

```bash
git clone https://github.com/yuki-koyama/mathtoolbox.git --recursive
cd mathtoolbox
mkdir build
cd build
cmake ../
make
```

--------------------------------

### Build and Serve Documentation Locally

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/UpdatingDocs.md

Generate the static site and start a local development server.

```bash
mkdocs build
mkdocs serve
```

--------------------------------

### Install Python Prerequisites

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/README.md

System-specific commands to install build dependencies for Python bindings.

```bash
brew install cmake eigen
```

```bash
sudo apt install cmake libeigen3-dev
```

--------------------------------

### Install Mathtoolbox C++ Library

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Install the mathtoolbox C++ library to the system after building it with CMake.

```bash
make install
```

--------------------------------

### Build Mathtoolbox C++ Library with Examples

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Build the mathtoolbox C++ library with example applications enabled using CMake. This requires setting the MATHTOOLBOX_BUILD_EXAMPLES parameter to ON.

```bash
cmake ../ -DMATHTOOLBOX_BUILD_EXAMPLES=ON
make
```

--------------------------------

### Install Eigen Dependencies

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/README.md

System-specific commands to install the required Eigen library.

```bash
brew install eigen
```

```bash
sudo apt install libeigen3-dev
```

--------------------------------

### Install Eigen on macOS

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Install the Eigen library on macOS using Homebrew.

```bash
brew install eigen
```

--------------------------------

### Install pymathtoolbox Python Library

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Install the pymathtoolbox Python library from GitHub using pip. This requires CMake and Eigen to be installed.

```bash
pip install git+https://github.com/yuki-koyama/mathtoolbox
```

--------------------------------

### Configure CMake for Optional Features

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/README.md

Commands to enable example applications or Python bindings during the CMake configuration step.

```bash
cmake ../ -DMATHTOOLBOX_BUILD_EXAMPLES=ON
make
```

```bash
cmake ../ -DMATHTOOLBOX_PYTHON_BINDINGS=ON
make
```

--------------------------------

### Install CMake and Eigen on macOS for Python

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Install CMake and Eigen on macOS using Homebrew, prerequisites for using pymathtoolbox.

```bash
brew install cmake eigen
```

--------------------------------

### Add Example Tests in CMake

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/CMakeLists.txt

Conditionally adds tests for Mathtoolbox examples if MATHTOOLBOX_BUILD_EXAMPLES is set. Each test is defined with a name and the command to execute the corresponding target.

```cmake
if(MATHTOOLBOX_BUILD_EXAMPLES)
  add_test(NAME acquisition-function-test COMMAND $<TARGET_FILE:acquisition-function-test>)
  add_test(NAME bayesian-optimization-test COMMAND $<TARGET_FILE:bayesian-optimization-test>)
  add_test(NAME bfgs-test COMMAND $<TARGET_FILE:bfgs-test>)
  add_test(NAME classical-mds-test COMMAND $<TARGET_FILE:classical-mds-test>)
  add_test(NAME gaussian-process-regression-test COMMAND $<TARGET_FILE:gaussian-process-regression-test> .)
  add_test(NAME gradient-descent-test COMMAND $<TARGET_FILE:gradient-descent-test>)
  add_test(NAME l-bfgs-test COMMAND $<TARGET_FILE:l-bfgs-test>)
  add_test(NAME log-determinant-test COMMAND $<TARGET_FILE:log-determinant-test>)
  add_test(NAME matrix-inversion-test COMMAND $<TARGET_FILE:matrix-inversion-test>)
  add_test(NAME probability-distributions-test COMMAND $<TARGET_FILE:probability-distributions-test>)
  add_test(NAME rbf-interpolation-test COMMAND $<TARGET_FILE:rbf-interpolation-test>)
  add_test(NAME som-test COMMAND $<TARGET_FILE:som-test>)
endif()
```

--------------------------------

### Install CMake and Eigen on Ubuntu for Python

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Install CMake and Eigen on Ubuntu 16.04/18.04 using apt, prerequisites for using pymathtoolbox.

```bash
sudo apt install cmake libeigen3-dev
```

--------------------------------

### Configure mathtoolbox CMake build

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/CMakeLists.txt

The primary CMake configuration file for the project, defining build options, library dependencies, and installation rules.

```cmake
cmake_minimum_required(VERSION 3.16)
set(CMAKE_POLICY_VERSION_MINIMUM 3.5)

project(mathtoolbox CXX)

set(CMAKE_CXX_STANDARD 11)
if(MATHTOOLBOX_PYTHON_BINDINGS)
  set(CMAKE_POSITION_INDEPENDENT_CODE ON)
  set(CMAKE_CXX_VISIBILITY_PRESET hidden)
endif()

######################################################################
# Options
######################################################################

option(MATHTOOLBOX_BUILD_EXAMPLES  "Build example applications" OFF)
option(MATHTOOLBOX_PYTHON_BINDINGS "Build python bindings"      OFF)

######################################################################
# Prerequisites
######################################################################

find_package(Eigen3 REQUIRED)
if((NOT TARGET Eigen3::Eigen) AND (DEFINED EIGEN3_INCLUDE_DIR))
  add_library(AliasEigen3 INTERFACE)
  target_include_directories(AliasEigen3 INTERFACE ${EIGEN3_INCLUDE_DIR})
  add_library(Eigen3::Eigen ALIAS AliasEigen3)
endif()

######################################################################
# External libraries
######################################################################

if(MATHTOOLBOX_PYTHON_BINDINGS)
  set(PYBIND11_INSTALL        OFF CACHE INTERNAL "" FORCE)
  set(PYBIND11_TEST           OFF CACHE INTERNAL "" FORCE)
  set(USE_PYTHON_INCLUDE_DIR  OFF CACHE INTERNAL "" FORCE)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/external/pybind11)
endif()

######################################################################
# Core library
######################################################################

file(GLOB headers ${CMAKE_CURRENT_SOURCE_DIR}/include/mathtoolbox/*.hpp)
file(GLOB sources ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp)

add_library(mathtoolbox STATIC ${headers} ${sources})
target_link_libraries(mathtoolbox PUBLIC Eigen3::Eigen)
target_include_directories(mathtoolbox PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)

######################################################################
# Python bindings
######################################################################

set(sources_python_bindings
  ${CMAKE_CURRENT_SOURCE_DIR}/src/pybind/python-bindings.cpp
)

if(MATHTOOLBOX_PYTHON_BINDINGS)
  pybind11_add_module(pymathtoolbox ${sources_python_bindings})
  target_link_libraries(pymathtoolbox PUBLIC mathtoolbox)
endif()

######################################################################
# Installation
######################################################################

install(FILES ${headers} DESTINATION include/mathtoolbox)
install(TARGETS mathtoolbox ARCHIVE DESTINATION lib)

######################################################################
# Example demos
######################################################################

if(MATHTOOLBOX_BUILD_EXAMPLES)
  set(OTF_WITH_EIGEN ON CACHE INTERNAL "" FORCE)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/external/optimization-test-functions)

  set(TIMER_BUILD_TEST OFF CACHE INTERNAL "" FORCE)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/external/timer)

  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/acquisition-function)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/bayesian-optimization)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/bfgs)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/classical-mds)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/gaussian-process-regression)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/gradient-descent)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/l-bfgs)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/log-determinant)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/matrix-inversion)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/probability-distributions)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/rbf-interpolation)
  add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/som)
endif()
```

--------------------------------

### Calculate Acquisition Functions for Bayesian Optimization

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Computes Expected Improvement and GP-UCB values to guide sampling in Bayesian optimization.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/acquisition-functions.hpp>

int main() {
    // Surrogate model predictions (typically from Gaussian process)
    auto mu = [](const Eigen::VectorXd& x) -> double {
        return -x.squaredNorm();  // Predicted mean
    };

    auto sigma = [](const Eigen::VectorXd& x) -> double {
        return 0.5 * (1.0 - x.norm());  // Predicted std (higher uncertainty far from data)
    };

    // Current best point
    Eigen::Vector2d x_best(0.1, 0.1);

    // Query point
    Eigen::Vector2d x(0.3, 0.4);

    // Expected Improvement (EI) - most common acquisition function
    double ei = mathtoolbox::GetExpectedImprovement(x, mu, sigma, x_best);
    std::cout << "Expected Improvement at x: " << ei << std::endl;

    // GP-UCB (Upper Confidence Bound)
    double beta = 2.0;  // Exploration-exploitation trade-off (higher = more exploration)
    double ucb = mathtoolbox::GetGaussianProcessUpperConfidenceBound(x, mu, sigma, beta);
    std::cout << "GP-UCB at x: " << ucb << std::endl;

    // Derivatives for gradient-based acquisition maximization
    auto mu_deriv = [](const Eigen::VectorXd& x) -> Eigen::VectorXd {
        return -2.0 * x;
    };

    auto sigma_deriv = [](const Eigen::VectorXd& x) -> Eigen::VectorXd {
        double n = x.norm();
        return (n > 1e-10) ? -0.5 * x / n : Eigen::VectorXd::Zero(x.size());
    };

    Eigen::VectorXd ei_grad = mathtoolbox::GetExpectedImprovementDerivative(
        x, mu, sigma, x_best, mu_deriv, sigma_deriv
    );
    std::cout << "EI gradient: " << ei_grad.transpose() << std::endl;

    return 0;
}
```

--------------------------------

### Get Inverse Using Upper Left Block Inverse

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/matrix-inversion.md

Calculates the inverse of a matrix using the pre-computed inverse of its upper-left block. This is particularly useful when the upper-left block is large and its inverse is already known.

```cpp
Eigen::MatrixXd GetInverseUsingUpperLeftBlockInverse(const Eigen::MatrixXd& matrix,
                                                     const Eigen::MatrixXd& upper_left_block_inverse);
```

--------------------------------

### Build mathtoolbox with CMake

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/README.md

Standard commands to clone, configure, and build the library from source.

```bash
git clone https://github.com/yuki-koyama/mathtoolbox.git --recursive
cd mathtoolbox
mkdir build
cd build
cmake ../
make
```

```bash
make install
```

--------------------------------

### Deploy Documentation to GitHub Pages

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/UpdatingDocs.md

Build and push the documentation to the gh-pages branch.

```bash
mkdocs gh-deploy
```

--------------------------------

### Instantiate GPR Object

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/gaussian-process-regression.md

Initialize the GaussianProcessRegressor with input data, kernel type, and normalization settings.

```cpp
GaussianProcessRegressor(const Eigen::MatrixXd& X,
                         const Eigen::VectorXd& y,
                         const KernelType       kernel_type = KernelType::ArdMatern52,
                         const bool             use_data_normalization = true);
```

--------------------------------

### Add Executable and Link Libraries

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/examples/bfgs/CMakeLists.txt

Defines the build target 'bfgs-test' and links it with the 'mathtoolbox' and 'optimization-test-functions' libraries. Ensure these libraries are available in your build environment.

```cmake
add_executable(bfgs-test main.cpp)
target_link_libraries(bfgs-test mathtoolbox optimization-test-functions)
```

--------------------------------

### GaussianProcessRegressor Constructor

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/gaussian-process-regression.md

Instantiates the GPR object with input data and configuration.

```APIDOC
## Constructor: GaussianProcessRegressor

### Description
Initializes the Gaussian Process Regressor with scattered data points.

### Parameters
- **X** (Eigen::MatrixXd) - Required - Data point locations in D-dimensional space.
- **y** (Eigen::VectorXd) - Required - Associated values for each data point.
- **kernel_type** (KernelType) - Optional - Kernel function type (default: KernelType::ArdMatern52).
- **use_data_normalization** (bool) - Optional - Whether to enable automatic data normalization (default: true).
```

--------------------------------

### Perform Bounded Gradient Descent in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Executes gradient descent with backtracking line search for box-constrained optimization. Initial points outside the specified bounds are projected into the feasible region.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/gradient-descent.hpp>

int main() {
    // Simple quadratic function
    auto f = [](const Eigen::VectorXd& x) -> double {
        return x.squaredNorm();
    };

    auto g = [](const Eigen::VectorXd& x) -> Eigen::VectorXd {
        return 2.0 * x;
    };

    // Set bounds: x in [-0.5, 0.5]^3
    Eigen::VectorXd lower_bound = Eigen::VectorXd::Constant(3, -0.5);
    Eigen::VectorXd upper_bound = Eigen::VectorXd::Constant(3, 0.5);

    // Initial point (outside bounds will be projected)
    Eigen::VectorXd x_init(3);
    x_init << 0.4, -0.3, 0.2;

    // Optimization parameters
    double epsilon = 1e-06;
    double default_alpha = 1.0;  // Initial step size
    unsigned int max_iters = 1000;

    // Output variables
    Eigen::VectorXd x_star;
    unsigned int num_iters;

    // Run bounded gradient descent
    mathtoolbox::optimization::RunGradientDescent(
        x_init, f, g, lower_bound, upper_bound,
        epsilon, default_alpha, max_iters,
        x_star, num_iters
    );

    std::cout << "Converged in " << num_iters << " iterations" << std::endl;
    std::cout << "Solution: " << x_star.transpose() << std::endl;

    return 0;
}
```

--------------------------------

### Perform Bayesian Optimization

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/bayesian-optimization.md

Define a target function, instantiate the optimizer, and run iterations to find the optimal solution.

```cpp
// Define the target problem
const int  num_dims    = 5;
const auto f           = [](const Eigen::VectorXd& x) { return - x.squaredNorm(); };
const auto lower_bound = Eigen::VectorXd::Constant(num_dims, -1.0);
const auto upper_bound = Eigen::VectorXd::Constant(num_dims, +1.0);

// Instantiate an optimizer
mathtoolbox::BayesianOptimizer optimizer{f, lower_bound, upper_bound};

// Perform optimization iteration
for (int i = 0; i < 50; ++i) { optimizer.Step(); }

// Retrieve the found solution
const Eigen::VectorXd x_star = optimizer.GetCurrentSolution();
```

--------------------------------

### Perform Bayesian Optimization in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Uses a Gaussian process surrogate model to maximize an objective function over defined bounds. Requires the Eigen library and mathtoolbox headers.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/bayesian-optimization.hpp>

int main() {
    // Define the objective function (maximization problem)
    constexpr int num_dims = 5;
    auto objective_func = [](const Eigen::VectorXd& x) -> double {
        return -x.squaredNorm();  // Maximize: find x closest to origin
    };

    // Define search bounds
    Eigen::VectorXd lower_bound = Eigen::VectorXd::Constant(num_dims, -1.0);
    Eigen::VectorXd upper_bound = Eigen::VectorXd::Constant(num_dims, 1.0);

    // Create optimizer with Matern 5/2 kernel (default)
    mathtoolbox::optimization::BayesianOptimizer optimizer(
        objective_func,
        lower_bound,
        upper_bound,
        mathtoolbox::optimization::KernelType::ArdMatern52
    );

    // Run optimization iterations
    constexpr int num_iterations = 20;
    for (int i = 0; i < num_iterations; ++i) {
        // Step returns (new_point, new_value) pair
        auto [new_point, new_value] = optimizer.Step();

        // Get current best solution
        Eigen::VectorXd current_best = optimizer.GetCurrentOptimizer();
        double current_best_value = optimizer.EvaluatePoint(current_best);

        std::cout << "Iter " << i + 1
                  << ": best_value=" << current_best_value
                  << std::endl;
    }

    // Retrieve final solution
    Eigen::VectorXd solution = optimizer.GetCurrentOptimizer();
    std::cout << "Final solution: " << solution.transpose() << std::endl;

    // Access observed data
    auto [X_data, y_data] = optimizer.GetData();
    std::cout << "Total evaluations: " << y_data.size() << std::endl;

    return 0;
}
```

--------------------------------

### Perform L-BFGS Optimization in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Demonstrates large-scale optimization using the L-BFGS algorithm. Requires defining objective and gradient functions and configuring the optimization settings.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/numerical-optimization.hpp>

int main() {
    // Large-scale optimization: 1000-dimensional quadratic function
    constexpr int dims = 1000;

    // f(x) = sum(i * x_i^2) - minimize to find x = 0
    auto f = [](const Eigen::VectorXd& x) -> double {
        double sum = 0.0;
        for (int i = 0; i < x.size(); ++i) {
            sum += (i + 1) * x(i) * x(i);
        }
        return sum;
    };

    auto g = [](const Eigen::VectorXd& x) -> Eigen::VectorXd {
        Eigen::VectorXd grad(x.size());
        for (int i = 0; i < x.size(); ++i) {
            grad(i) = 2.0 * (i + 1) * x(i);
        }
        return grad;
    };

    // Use L-BFGS for large-scale problem
    mathtoolbox::optimization::Setting setting;
    setting.algorithm = mathtoolbox::optimization::Algorithm::LBfgs;
    setting.x_init = Eigen::VectorXd::Random(dims);
    setting.f = f;
    setting.g = g;
    setting.type = mathtoolbox::optimization::Type::Min;
    setting.epsilon = 1e-06;
    setting.max_num_iterations = 500;

    mathtoolbox::optimization::Result result = mathtoolbox::optimization::RunOptimization(setting);

    std::cout << "L-BFGS converged in " << result.num_iterations << " iterations" << std::endl;
    std::cout << "Final objective value: " << f(result.x_star) << std::endl;
    std::cout << "Solution norm: " << result.x_star.norm() << std::endl;

    return 0;
}
```

--------------------------------

### Include Matrix Inversion Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/matrix-inversion.md

Include the necessary header file for matrix inversion functionalities.

```cpp
#include <mathtoolbox/matrix-inversion.hpp>
```

--------------------------------

### Train a Self-Organizing Map (SOM) in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Initializes and trains a SOM on random 3D data. Requires the Eigen library and mathtoolbox/som.hpp header.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/som.hpp>

int main() {
    // Generate random 3D data (e.g., RGB color values)
    constexpr int num_data_dims = 3;
    constexpr int num_points = 100;
    Eigen::MatrixXd data = Eigen::MatrixXd::Random(num_data_dims, num_points);

    // SOM parameters
    constexpr int latent_dims = 2;        // Map dimensionality (1 or 2)
    constexpr int resolution = 10;         // Grid resolution
    constexpr double init_var = 0.20;      // Initial neighborhood variance
    constexpr double min_var = 0.01;       // Minimum variance
    constexpr double var_speed = 20.0;     // Variance decreasing speed
    constexpr bool normalize = true;       // Normalize input data

    // Create SOM
    mathtoolbox::Som som(data, latent_dims, resolution, init_var, min_var, var_speed, normalize);

    // Train for 100 iterations
    for (int i = 0; i < 100; ++i) {
        Eigen::MatrixXd prev_nodes = som.GetDataSpaceNodePositions();
        som.Step();
        Eigen::MatrixXd curr_nodes = som.GetDataSpaceNodePositions();

        double delta = (prev_nodes - curr_nodes).norm();
        if (i % 20 == 0) {
            std::cout << "Iteration " << i + 1 << ": delta=" << delta << std::endl;
        }
    }

    // Get learned manifold nodes in data space
    Eigen::MatrixXd node_positions = som.GetDataSpaceNodePositions();
    std::cout << "Node positions shape: " << node_positions.rows() << " x " << node_positions.cols() << std::endl;

    // Get data point positions in latent space
    Eigen::MatrixXd latent_positions = som.GetLatentSpaceDataPositions();
    std::cout << "Latent positions shape: " << latent_positions.rows() << " x " << latent_positions.cols() << std::endl;

    return 0;
}
```

--------------------------------

### Include Log-Determinant Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/log-determinant.md

Include the necessary header file for log-determinant calculations.

```cpp
#include <mathtoolbox/log-determinant.hpp>
```

--------------------------------

### Build Mathtoolbox C++ Library with Python Bindings

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/index.md

Build the mathtoolbox C++ library with Python bindings enabled using CMake. This requires setting the MATHTOOLBOX_PYTHON_BINDINGS parameter to ON.

```bash
cmake ../ -DMATHTOOLBOX_PYTHON_BINDINGS=ON
make
```

--------------------------------

### Perform Backtracking Line Search in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Calculates an appropriate step size using the Armijo condition. Supports both unconstrained and bounded search variants.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/backtracking-line-search.hpp>

int main() {
    // Define objective function
    auto f = [](const Eigen::VectorXd& x) -> double {
        return x.squaredNorm();
    };

    // Current point and descent direction
    Eigen::Vector2d x(1.0, 1.0);
    Eigen::Vector2d grad(2.0, 2.0);  // Gradient at x
    Eigen::Vector2d p = -grad;        // Descent direction (negative gradient)

    // Line search parameters
    double alpha_init = 1.0;  // Initial step size
    double rho = 0.5;         // Reduction factor
    double c = 1e-04;         // Armijo condition parameter

    // Find step size
    double alpha = mathtoolbox::optimization::RunBacktrackingLineSearch(
        f, grad, x, p, alpha_init, rho, c
    );

    std::cout << "Step size: " << alpha << std::endl;
    std::cout << "New point: " << (x + alpha * p).transpose() << std::endl;
    std::cout << "New value: " << f(x + alpha * p) << std::endl;

    // With bounds
    Eigen::Vector2d lower(-0.5, -0.5);
    Eigen::Vector2d upper(0.5, 0.5);

    double alpha_bounded = mathtoolbox::optimization::RunBacktrackingBoundedLineSearch(
        f, grad, x, p, lower, upper, alpha_init, rho, c
    );

    std::cout << "Bounded step size: " << alpha_bounded << std::endl;

    return 0;
}
```

--------------------------------

### Hyperparameter Configuration

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/gaussian-process-regression.md

Methods for setting hyperparameters either manually or via optimization.

```APIDOC
## Method: SetHyperparams

### Description
Manually sets the kernel and noise hyperparameters.

### Parameters
- **kernel_hyperparams** (Eigen::VectorXd) - Required - Signal variance and length-scales.
- **noise_hyperparam** (double) - Required - Noise variance.

## Method: PerformMaximumLikelihood

### Description
Determines optimal hyperparameters using the L-BFGS optimization algorithm.

### Parameters
- **kernel_hyperparams_initial** (Eigen::VectorXd) - Required - Initial guess for kernel hyperparameters.
- **noise_hyperparam_initial** (double) - Required - Initial guess for noise variance.
```

--------------------------------

### Include Bayesian Optimization Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/bayesian-optimization.md

Include the primary header file to access the Bayesian optimization functionality.

```cpp
#include <mathtoolbox/bayesian-optimization.hpp>
```

--------------------------------

### Set Hyperparameters

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/gaussian-process-regression.md

Configure hyperparameters either manually or through maximum likelihood estimation.

```cpp
void SetHyperparams(const Eigen::VectorXd& kernel_hyperparams, const double noise_hyperparam);
```

```cpp
void PerformMaximumLikelihood(const Eigen::VectorXd& kernel_hyperparams_initial,
                              const double           noise_hyperparam_initial);
```

--------------------------------

### Define Executable Target

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/examples/bayesian-optimization/CMakeLists.txt

Defines the main executable for the project and its source file. Use this to set up the primary build target.

```cmake
add_executable(bayesian-optimization-test main.cpp)
```

--------------------------------

### Perform Gaussian Process Regression in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Provides probabilistic function approximation with uncertainty quantification. Supports hyperparameter optimization via maximum likelihood estimation.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/gaussian-process-regression.hpp>

int main() {
    // Generate training data
    constexpr int num_samples = 20;
    Eigen::MatrixXd X(1, num_samples);  // 1D input
    Eigen::VectorXd y(num_samples);

    for (int i = 0; i < num_samples; ++i) {
        X(0, i) = static_cast<double>(i) / num_samples;
        y(i) = X(0, i) * std::sin(10.0 * X(0, i)) + 0.01 * (rand() / (double)RAND_MAX);
    }

    // Create regressor with Matern 5/2 kernel
    auto kernel_type = mathtoolbox::GaussianProcessRegressor::KernelType::ArdMatern52;
    // Alternative: KernelType::ArdSquaredExp

    mathtoolbox::GaussianProcessRegressor regressor(X, y, kernel_type, true);

    // Perform hyperparameter optimization via maximum likelihood
    Eigen::Vector2d initial_kernel_params(0.5, 0.5);
    double initial_noise_param = 0.01;
    regressor.PerformMaximumLikelihood(initial_kernel_params, initial_noise_param);

    // Predict at new points with uncertainty
    for (double x = 0.0; x <= 1.0; x += 0.1) {
        Eigen::VectorXd query = Eigen::VectorXd::Constant(1, x);
        double mean = regressor.PredictMean(query);
        double stdev = regressor.PredictStdev(query);

        std::cout << "x=" << x
                  << " mean=" << mean
                  << " 95% CI=[" << mean - 1.96 * stdev << ", " << mean + 1.96 * stdev << "]"
                  << std::endl;
    }

    return 0;
}
```

--------------------------------

### Include GPR Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/gaussian-process-regression.md

Include the primary header file to access Gaussian process regression functionality.

```cpp
#include <mathtoolbox/gaussian-process-regression.hpp>
```

--------------------------------

### Include Strong Wolfe Conditions Line Search Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/strong-wolfe-conditions-line-search.md

Include the necessary header file for using the strong Wolfe conditions line search functionality.

```cpp
#include <mathtoolbox/strong-wolfe-conditions-line-search.hpp>
```

--------------------------------

### Include Mathtoolbox Constants Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/constants.md

Include the necessary header file to use constants from the mathtoolbox library. Constants are defined within the mathtoolbox::constants namespace.

```cpp
#include <mathtoolbox/constants.hpp>
```

--------------------------------

### Perform RBF Interpolation, Bayesian Optimization, and Classical MDS

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Demonstrates the usage of pymathtoolbox for scattered data interpolation, black-box function optimization, and dimensionality reduction.

```python
import pymathtoolbox
import numpy as np
from scipy.spatial.distance import pdist, squareform

# Set random seed for reproducibility
pymathtoolbox.set_seed(42)

# ============ RBF Interpolation ============
# Generate noisy samples
num_samples = 50
X = np.random.rand(1, num_samples)  # 1D input
y = np.sin(X[0] * 10) + 0.1 * np.random.randn(num_samples)

# Create interpolator with Gaussian RBF kernel
kernel = pymathtoolbox.GaussianRbfKernel(theta=5.0)
interpolator = pymathtoolbox.RbfInterpolator(kernel)

# Set data and compute weights with regularization
interpolator.set_data(X, y)
interpolator.calc_weights(use_regularization=True, regularization_weight=1e-04)

# Interpolate at new points
x_test = np.linspace(0, 1, 100)
y_interp = [interpolator.calc_value(np.array([x])) for x in x_test]
print(f"Interpolated {len(y_interp)} points")

# ============ Bayesian Optimization ============
def objective(x):
    """Objective function to maximize"""
    return -np.sum(x ** 2)  # Maximum at origin

lower_bound = np.array([-1.0, -1.0])
upper_bound = np.array([1.0, 1.0])

optimizer = pymathtoolbox.BayesianOptimizer(objective, lower_bound, upper_bound)

# Run optimization
for i in range(15):
    optimizer.step()
    current_best = optimizer.get_current_optimizer()
    print(f"Iter {i+1}: best = {current_best}, value = {objective(current_best):.4f}")

# Get final result
X_observed, y_observed = optimizer.get_data()
print(f"Total evaluations: {len(y_observed)}")

# ============ Classical MDS ============
# Create high-dimensional points
points = [
    np.array([0.0, 2.0, 0.0, 3.0, 4.0]),
    np.array([1.0, 0.0, 2.0, 4.0, 3.0]),
    np.array([0.0, 1.0, 4.0, 2.0, 0.0]),
    np.array([0.0, 4.0, 1.0, 0.0, 2.0]),
    np.array([4.0, 3.0, 0.0, 1.0, 0.0]),
]

# Compute distance matrix
D = squareform(pdist(points))

# Embed into 2D
X_2d = pymathtoolbox.compute_classical_mds(D=D, target_dim=2)
print(f"\n2D Embedding shape: {X_2d.shape}")
print(X_2d)
```

--------------------------------

### Perform Strong Wolfe Conditions Line Search in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Executes a line search satisfying strong Wolfe conditions to determine an appropriate step size for quasi-Newton methods.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/strong-wolfe-conditions-line-search.hpp>

int main() {
    // Define objective and gradient
    auto f = [](const Eigen::VectorXd& x) -> double {
        return x(0) * x(0) + 4.0 * x(1) * x(1);  // Elliptic paraboloid
    };

    auto g = [](const Eigen::VectorXd& x) -> Eigen::VectorXd {
        Eigen::Vector2d grad;
        grad << 2.0 * x(0), 8.0 * x(1);
        return grad;
    };

    // Current point and search direction
    Eigen::Vector2d x(2.0, 1.0);
    Eigen::Vector2d p = -g(x);  // Negative gradient direction

    // Strong Wolfe line search parameters
    double alpha_init = 1.0;
    double alpha_max = 50.0;
    double c1 = 1e-04;  // Sufficient decrease parameter
    double c2 = 0.9;    // Curvature condition parameter

    double alpha = mathtoolbox::optimization::RunStrongWolfeConditionsLineSearch(
        f, g, x, p, alpha_init, alpha_max, c1, c2
    );

    std::cout << "Wolfe step size: " << alpha << std::endl;
    std::cout << "New point: " << (x + alpha * p).transpose() << std::endl;

    return 0;
}
```

--------------------------------

### Link Libraries to Executable

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/examples/bayesian-optimization/CMakeLists.txt

Links the specified executable to its required libraries. Ensure these libraries are available in your build environment.

```cmake
target_link_libraries(bayesian-optimization-test mathtoolbox optimization-test-functions)
```

--------------------------------

### Enable Testing in CMake

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/CMakeLists.txt

This CMake command enables the testing infrastructure for the project. It should be called before adding any tests.

```cmake
enable_testing()
```

--------------------------------

### Include Data Normalization Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/data-normalization.md

Required header file to access data normalization functionality.

```cpp
#include <mathtoolbox/data-normalization.hpp>
```

--------------------------------

### Perform BFGS Optimization in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Minimizes the Rosenbrock function using the BFGS quasi-Newton method. Requires defining the objective function and its gradient.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/numerical-optimization.hpp>

int main() {
    // Define the Rosenbrock function (minimization target)
    auto f = [](const Eigen::VectorXd& x) -> double {
        double sum = 0.0;
        for (int i = 0; i < x.size() - 1; ++i) {
            sum += 100.0 * std::pow(x(i + 1) - x(i) * x(i), 2) + std::pow(1.0 - x(i), 2);
        }
        return sum;
    };

    // Define the gradient
    auto g = [](const Eigen::VectorXd& x) -> Eigen::VectorXd {
        Eigen::VectorXd grad = Eigen::VectorXd::Zero(x.size());
        for (int i = 0; i < x.size() - 1; ++i) {
            grad(i) += -400.0 * x(i) * (x(i + 1) - x(i) * x(i)) - 2.0 * (1.0 - x(i));
            grad(i + 1) += 200.0 * (x(i + 1) - x(i) * x(i));
        }
        return grad;
    };

    // Set up optimization using the unified interface
    mathtoolbox::optimization::Setting setting;
    setting.algorithm = mathtoolbox::optimization::Algorithm::Bfgs;
    setting.x_init = Eigen::VectorXd::Random(10);  // 10-dimensional problem
    setting.f = f;
    setting.g = g;
    setting.type = mathtoolbox::optimization::Type::Min;
    setting.epsilon = 1e-05;
    setting.max_num_iterations = 1000;

    // Run optimization
    mathtoolbox::optimization::Result result = mathtoolbox::optimization::RunOptimization(setting);

    std::cout << "Converged in " << result.num_iterations << " iterations" << std::endl;
    std::cout << "Solution: " << result.x_star.transpose() << std::endl;
    std::cout << "Final value: " << f(result.x_star) << std::endl;
    // Expected: solution near (1, 1, ..., 1) with value near 0

    return 0;
}
```

--------------------------------

### Compute Kernel Functions and Derivatives in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Calculates ARD kernel values and their derivatives for Gaussian processes and hyperparameter optimization.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/kernel-functions.hpp>

int main() {
    // Two input points
    Eigen::Vector3d x_a(1.0, 2.0, 3.0);
    Eigen::Vector3d x_b(1.5, 2.5, 2.5);

    // Hyperparameters: [sigma_f, length_scale_1, length_scale_2, length_scale_3]
    // sigma_f controls overall variance, length scales control smoothness per dimension
    Eigen::Vector4d theta(1.0, 0.5, 0.5, 0.5);

    // ARD Squared Exponential kernel
    double k_se = mathtoolbox::GetArdSquaredExpKernel(x_a, x_b, theta);
    std::cout << "ARD Squared Exp kernel: " << k_se << std::endl;

    // ARD Matern 5/2 kernel (recommended for optimization)
    double k_m52 = mathtoolbox::GetArdMatern52Kernel(x_a, x_b, theta);
    std::cout << "ARD Matern 5/2 kernel: " << k_m52 << std::endl;

    // Kernel derivatives for hyperparameter optimization
    Eigen::VectorXd dk_dtheta = mathtoolbox::GetArdSquaredExpKernelThetaDerivative(x_a, x_b, theta);
    std::cout << "Derivative w.r.t. theta: " << dk_dtheta.transpose() << std::endl;

    // Derivative w.r.t. specific hyperparameter index
    double dk_dtheta_0 = mathtoolbox::GetArdSquaredExpKernelThetaIDerivative(x_a, x_b, theta, 0);
    std::cout << "Derivative w.r.t. theta[0]: " << dk_dtheta_0 << std::endl;

    // Derivative w.r.t. first argument (for gradient-based optimization)
    Eigen::VectorXd dk_dx = mathtoolbox::GetArdSquaredExpKernelFirstArgDerivative(x_a, x_b, theta);
    std::cout << "Derivative w.r.t. x_a: " << dk_dx.transpose() << std::endl;

    return 0;
}
```

--------------------------------

### Probability Distributions

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Calculates PDF, CDF, and derivatives for standard normal, general normal, and log-normal distributions. Also computes the PDF for multivariate normal distributions.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/probability-distributions.hpp>

int main() {
    double x = 0.5;

    // Standard normal distribution N(0, 1)
    double pdf = mathtoolbox::GetStandardNormalDist(x);
    double cdf = mathtoolbox::GetStandardNormalDistCdf(x);
    double pdf_deriv = mathtoolbox::GetStandardNormalDistDerivative(x);

    std::cout << "Standard Normal at x=" << x << std::endl;
    std::cout << "  PDF: " << pdf << std::endl;
    std::cout << "  CDF: " << cdf << std::endl;
    std::cout << "  PDF derivative: " << pdf_deriv << std::endl;

    // General normal distribution N(mu, sigma^2)
    double mu = 1.0;
    double sigma_2 = 0.25;  // variance

    double normal_pdf = mathtoolbox::GetNormalDist(x, mu, sigma_2);
    double normal_pdf_deriv = mathtoolbox::GetNormalDistDerivative(x, mu, sigma_2);

    std::cout << "\nNormal(1.0, 0.25) at x=" << x << std::endl;
    std::cout << "  PDF: " << normal_pdf << std::endl;
    std::cout << "  PDF derivative: " << normal_pdf_deriv << std::endl;

    // Log-normal distribution
    double lognormal_pdf = mathtoolbox::GetLogNormalDist(x, mu, sigma_2);
    double log_lognormal = mathtoolbox::GetLogOfLogNormalDist(x, mu, sigma_2);

    std::cout << "\nLogNormal(1.0, 0.25) at x=" << x << std::endl;
    std::cout << "  PDF: " << lognormal_pdf << std::endl;
    std::cout << "  Log PDF: " << log_lognormal << std::endl;

    // Multivariate normal distribution
    Eigen::Vector2d x_mv(0.5, 0.3);
    Eigen::Vector2d mu_mv(0.0, 0.0);
    Eigen::Matrix2d Sigma;
    Sigma << 1.0, 0.5, 0.5, 1.0;
    Eigen::Matrix2d Sigma_inv = Sigma.inverse();
    double Sigma_det = Sigma.determinant();

    double mv_pdf = mathtoolbox::GetNormalDist(x_mv, mu_mv, Sigma_inv, Sigma_det);
    std::cout << "\nMultivariate Normal PDF: " << mv_pdf << std::endl;

    return 0;
}

```

--------------------------------

### RbfInterpolator Class Methods

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/rbf-interpolation.md

Methods for configuring, training, and using the RbfInterpolator for scattered data interpolation.

```APIDOC
## void SetData(const Eigen::MatrixXd& X, const Eigen::VectorXd& y)

### Description
Sets the target scattered data points and their corresponding values for the interpolator.

### Parameters
#### Request Body
- **X** (Eigen::MatrixXd) - Required - Data points matrix of size m x n.
- **y** (Eigen::VectorXd) - Required - Data values vector of size m.

---

## void ComputeWeights(const bool use_regularization, const double lambda)

### Description
Calculates the weight values for the RBF interpolation. Setting use_regularization to true enables approximation mode, which is useful for noisy data.

### Parameters
#### Request Body
- **use_regularization** (bool) - Optional - If true, performs scattered data approximation (default: false).
- **lambda** (double) - Optional - Regularization parameter (default: 0.001).

---

## double CalcValue(const Eigen::VectorXd& x)

### Description
Calculates the interpolated value for a given input vector x using the computed weights.

### Parameters
#### Request Body
- **x** (Eigen::VectorXd) - Required - Input vector for which to calculate the interpolated value.
```

--------------------------------

### Include Classical MDS Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/classical-mds.md

Include the necessary header file for using the classical MDS functionality.

```cpp
#include <mathtoolbox/classical-mds.hpp>
```

--------------------------------

### Set Scattered Data Points

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/rbf-interpolation.md

Use this method to set the target scattered data points (X) and their corresponding values (y) for interpolation.

```cpp
void SetData(const Eigen::MatrixXd& X, const Eigen::VectorXd& y);
```

--------------------------------

### Include RBF Interpolation Header

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/rbf-interpolation.md

Include the necessary header file to access RBF interpolation functionality.

```cpp
#include <mathtoolbox/rbf-interpolation.hpp>
```

--------------------------------

### Predict Standard Deviation with GPR

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/gaussian-process-regression.md

Call this method to predict the standard deviation $ \sqrt{\text{Var}(f)} $ for an unknown location $ \mathbf{x} $ after a GPR object is instantiated and its hyperparameters are set.

```cpp
double PredictStdev(const Eigen::VectorXd& x) const;
```

--------------------------------

### Perform RBF Interpolation in C++

Source: https://context7.com/yuki-koyama/mathtoolbox/llms.txt

Uses radial basis functions to perform scattered data interpolation. Requires setting data points and calculating weights before querying interpolated values.

```cpp
#include <Eigen/Core>
#include <mathtoolbox/rbf-interpolation.hpp>

int main() {
    // Generate scattered data: 100 points in 2D space
    constexpr int num_samples = 100;
    Eigen::MatrixXd X = Eigen::MatrixXd::Random(2, num_samples);  // 2 x 100 matrix
    Eigen::VectorXd y(num_samples);

    // Sample function: sin(10*x) + sin(10*y)
    for (int i = 0; i < num_samples; ++i) {
        y(i) = std::sin(10.0 * X(0, i)) + std::sin(10.0 * X(1, i));
    }

    // Create interpolator with thin plate spline kernel (default)
    mathtoolbox::RbfInterpolator interpolator(mathtoolbox::ThinPlateSplineRbfKernel());

    // Alternative kernels:
    // mathtoolbox::GaussianRbfKernel(1.0)      - Gaussian with theta=1.0
    // mathtoolbox::LinearRbfKernel()           - Linear kernel
    // mathtoolbox::CubicRbfKernel()            - Cubic kernel

    // Set data points and values
    interpolator.SetData(X, y);

    // Calculate weights (with regularization for noisy data)
    bool use_regularization = true;
    double lambda = 0.001;  // Regularization parameter
    interpolator.CalcWeights(use_regularization, lambda);

    // Query interpolated value at a new point
    Eigen::Vector2d query_point(0.5, 0.3);
    double interpolated_value = interpolator.CalcValue(query_point);

    std::cout << "Interpolated value at (0.5, 0.3): " << interpolated_value << std::endl;
    // Expected output: value close to sin(5.0) + sin(3.0) ≈ -0.818

    return 0;
}
```

--------------------------------

### Predict Mean with GPR

Source: https://github.com/yuki-koyama/mathtoolbox/blob/master/docs/gaussian-process-regression.md

Call this method to predict the most likely value for an unknown location $ \mathbf{x} $ after a GPR object is instantiated and its hyperparameters are set.

```cpp
double PredictMean(const Eigen::VectorXd& x) const;
```