Initial commit. Python and C++ version using Eigen.

.gitignore

@@ -0,0 +1,2 @@
+cpp/build/*
+python/__pycache__/*

cpp/CMakeLists.txt

@@ -0,0 +1,22 @@
+cmake_minimum_required(VERSION 3.0)
+
+project(CubicLagrangeMinimize)
+
+# Set the C++ standard to C++11
+set(CMAKE_CXX_STANDARD 11)
+
+# Find the Eigen library
+find_package(Eigen3 REQUIRED)
+include_directories(${EIGEN3_INCLUDE_DIRS})
+
+# Add the source files
+set(SOURCES
+    src/tests.cpp
+    src/CubicLagrangeMinimize.cpp
+)
+
+# Add an executable target
+add_executable(CubicLagrangeMinimize ${SOURCES})
+
+# Link against the Eigen library
+target_link_libraries(CubicLagrangeMinimize Eigen3::Eigen)

cpp/Makefile

@@ -0,0 +1,40 @@
+# Makefile for compiling C++ project
+
+# Compiler
+CXX = g++
+
+# Compiler options
+CXXFLAGS = -std=c++11 -Wall -Wextra
+
+# Linker options
+LDFLAGS =
+
+# Include path
+INC_PATH = -I./include -ID:/Users/Jerome/Documents/Ingenierie/Programmation/eigen-3.4.0
+
+# Source path
+SRC_PATH = ./src
+
+# Build path
+BUILD_PATH = ./build
+
+# Source files
+SRCS := $(wildcard $(SRC_PATH)/*.cpp)
+
+# Object files
+OBJS := $(patsubst $(SRC_PATH)/%.cpp,$(BUILD_PATH)/%.o,$(SRCS))
+
+# Executable
+EXEC = $(BUILD_PATH)/CubicLagrangeMinimize
+
+# Targets
+all: $(EXEC)
+
+$(EXEC): $(OBJS)
+	$(CXX) $(LDFLAGS) $(OBJS) -o $(EXEC)
+
+$(BUILD_PATH)/%.o: $(SRC_PATH)/%.cpp
+	$(CXX) $(CXXFLAGS) $(INC_PATH) -c $< -o $@
+
+clean:
+	rm -f $(OBJS) $(EXEC)

cpp/include/CubicLagrangeMinimize.hpp

@@ -0,0 +1,24 @@
+#ifndef DEF_CubicLagrangeMinimize
+#define DEF_CubicLagrangeMinimize
+
+#include <cmath>
+#include <Eigen/Dense>
+
+using Eigen::MatrixXd;
+using Eigen::VectorXd;
+using Eigen::Matrix4d;
+using Eigen::Vector4d;
+
+/// @brief Function to find minimum of f over interval [a, b] using cubic Lagrange polynomial interpolation.
+/// If the function is monotonic, then the minimum is one of the bounds of the interval, and the minimum is found in a single iteration.
+/// The best estimate of the minimum within the current interval is returned once the interval is smaller than the tolerance.
+/// The number of function evaluations is 2 + 2*Niter.
+/// The interval width is reduced by a factor of 3 every iteration.
+/// @param f The univariate function to minimize.
+/// @param a The lower bound of the interval.
+/// @param b The upper bound of the interval.
+/// @param tol The tolerance on the interval width.
+/// @return The best estimate of the minimum within the current interval once its width is smaller than the tolerance.
+double CubicLagrangeMinimize(std::function<double(double)> f, double a, double b, double tol=1e-9);
+
+#endif

cpp/src/CubicLagrangeMinimize.cpp

@@ -0,0 +1,107 @@
+#include <CubicLagrangeMinimize.hpp>
+
+/// @brief Linear least-squares fit of a cubic polynomial to four points.
+/// @param x1 Abscissa of first point.
+/// @param x2 Abscissa of second point.
+/// @param x3 Abscissa of third point.
+/// @param x4 Abscissa of fourth point.
+/// @param y1 Ordinate of first point.
+/// @param y2 Ordinate of second point.
+/// @param y3 Ordinate of third point.
+/// @param y4 Ordinate of fourth point.
+/// @return 4-dimensional vector of polynomial coefficients from low to high order : [a0 a1 a2 a3] -> a0 + a1*x + a2*x^2 + a3*x^3
+Vector4d polyfit4(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4) {
+    Matrix4d A(4, 4);
+    double x1x1 = x1*x1, x2x2 = x2*x2, x3x3 = x3*x3, x4x4 = x4*x4;
+    double x1x1x1 = x1x1*x1, x2x2x2 = x2x2*x2, x3x3x3 = x3x3*x3, x4x4x4 = x4x4*x4;
+    A << 1, x1, x1x1, x1x1x1,
+         1, x2, x2x2, x2x2x2,
+         1, x3, x3x3, x3x3x3,
+         1, x4, x4x4, x4x4x4;
+    VectorXd y(4);
+    y << y1, y2, y3, y4;
+    return (A.transpose() * A).colPivHouseholderQr().solve(A.transpose() * y);
+}
+
+/// @brief Cubic polynomial function.
+/// @param x    Abscissa.
+/// @param a0   Constant term.
+/// @param a1   Linear term.
+/// @param a2   Quadratic term.
+/// @param a3   Cubic term.
+/// @return a0 + a1*x + a2*x^2 + a3*x^3
+double cubic_poly(double x, double a0, double a1, double a2, double a3) {
+    return a0 + a1*x + a2*x*x + a3*x*x*x;
+}
+
+double CubicLagrangeMinimize(std::function<double(double)> f, double a, double b, double tol) {
+    // initialize interval endpoints and function values
+    double x0 = a, x1 = a*2./3. + b*1./3., x2 = a*1./3. + b*2./3., x3 = b;  // endpoints and two points in the interval
+    double f0 = f(x0), f1 = 0., f2 = 0., f3 = f(x3);                        // function values at the endpoints and two points in the interval
+    double x_prev = x0;                                                     // previous value of x_sol to track convergence of the solution
+    double a0, a1, a2, a3;                                                  // coefficients of the cubic Lagrange polynomial
+    double delta;                                                           // determinant of the quadratic equation of the Lagrange polynomial
+    double x_sol, x_sol_1, x_sol_2, y_sol, y_sol_1, y_sol_2;                // solutions of the Lagrange polynomial
+    constexpr double small_coefficient = 1e-9;                              // threshold for small coefficients to avoid ill-conditioning of the quadratic equation
+
+    unsigned int Niter = static_cast<unsigned int>(std::ceil(std::log(std::fabs(b - a)/tol)/std::log(3.)));// number of iterations to reduce interval width by a factor of 3
+
+    for(unsigned int i = 0 ; i < Niter ; ++i) {
+        // Compute function values at two points in the interval
+        x1 = x0*2./3. + x3*1./3., x2 = x0*1./3. + x3*2./3.;
+        f1 = f(x1), f2 = f(x2);
+
+        // compute Lagrange polynomial using least-squares fit to 4 points, which is equivalent to the cubic Lagrange polynomial
+        Vector4d A = polyfit4(x0, x1, x2, x3, f0, f1, f2, f3);
+        a0 = A[0]; a1 = A[1]; a2 = A[2]; a3 = A[3];
+
+        // Solve the first derivative of the Lagrange polynomial for a zero
+        if(std::fabs(a3) > small_coefficient) {
+            delta = -3*a1*a3 + a2*a2;
+
+            if(delta < 0) {
+                x_sol = (f0 < f3) ? x0 : x3;   // just choose the interval tha contains the minimum of the linear polynomial
+                y_sol = cubic_poly(x_sol, a0, a1, a2, a3);
+            } else {                             // solve for the two solutions of the quadratic equation of the first derivative of the Lagrange polynomial
+                x_sol_1 = (-a2 + std::sqrt(delta))/(3.*a3);
+                x_sol_2 = (-a2 - std::sqrt(delta))/(3.*a3);
+
+                y_sol_1 = cubic_poly(x_sol_1, a0, a1, a2, a3);
+                y_sol_2 = cubic_poly(x_sol_2, a0, a1, a2, a3);
+
+                x_sol = (y_sol_1 < y_sol_2) ? x_sol_1 : x_sol_2;
+                y_sol = (y_sol_1 < y_sol_2) ? y_sol_1 : y_sol_2;
+            }
+        }
+        else if(std::fabs(a2) > small_coefficient) { // if a3 is zero, then the Lagrange polynomial is a quadratic polynomial
+            x_sol = -a1/(2.*a2);
+            y_sol = cubic_poly(x_sol, a0, a1, a2, a3);
+        } else {                                     // if a3 and a2 are zero, then the Lagrange polynomial is a linear polynomial
+            x_sol = (f0 < f3) ? x0 : x3;             // just choose the interval tha contains the minimum of the linear polynomial
+            y_sol = cubic_poly(x_sol, a0, a1, a2, a3);
+        }
+
+        // Check convergence
+        if(std::fabs(x_sol - x_prev) < tol) { break; }
+
+        // Determine which interval contains the minimum of the cubic polynomial
+        if(x_sol < x1) {
+            x3 = x1; f3 = f1;
+        } else if(x_sol < x2) {
+            x0 = x1; f0 = f1;
+            x3 = x2; f3 = f2;
+        } else {
+            x0 = x2; f0 = f2;
+        }
+        
+        x_prev = x_sol;
+    }
+
+    // return best estimate of minimum
+    if((y_sol < f0) && (y_sol < f3))
+        return x_sol;
+    else if(f0 < f3)
+        return x0;
+    else
+        return x3;
+}

cpp/src/tests.cpp

@@ -0,0 +1,80 @@
+#include <iostream>
+#include <vector>
+#include <CubicLagrangeMinimize.hpp>
+
+using std::cout;
+using std::endl;
+
+double fct_01(double x) {
+    return x*x + std::sin(5.*x);
+}
+
+double dfct_01(double x) {
+    return 2.*x + 5.*std::cos(5.*x);
+}
+
+double fct_02(double x) {
+    return std::exp(x);
+}
+
+double dfct_02(double x) {
+    return std::exp(x);
+}
+
+double fct_03(double x) {
+    return -std::exp(-x*x);
+}
+
+double dfct_03(double x) {
+    return 2.*x*std::exp(-x*x);
+}
+
+double fct_04(double x) {
+    return std::exp(-x);
+}
+
+double dfct_04(double x) {
+    return std::exp(-x);
+}
+
+double fct_05(double x) {
+    return std::sin(x);
+}
+
+double dfct_05(double x) {
+    return std::cos(x);
+}
+
+double fct_06(double x) {
+    return x*x;
+}
+
+double dfct_06(double x) {
+    return 2.*x;
+}
+
+double fct_07(double x) {
+    return x;
+}
+
+double dfct_07(double) {
+    return 1.;
+}
+
+void test_CubicLagrangeMinimize() {
+    std::vector<std::function<double(double)>> fcts = {fct_01, fct_02, fct_03, fct_04, fct_05, fct_06, fct_07};
+    std::vector<std::function<double(double)>> dfcts = {dfct_01, dfct_02, dfct_03, dfct_04, dfct_05, dfct_06, dfct_07};
+    std::vector<double> mins = {-1.2, -10., -1.5, -10., -1., -2., -3.};
+    std::vector<double> maxs = {1.5, 10., 3., 5., 6., 3., 4.};
+    for(unsigned int i = 0 ; i < fcts.size() ; ++i) {
+        auto f = fcts[i]; auto df = dfcts[i];
+        auto x_min = CubicLagrangeMinimize(f, mins[i], maxs[i]);
+        cout << "[" << mins[i] << " " << maxs[i] << "]\tf(" << x_min << ") = " << f(x_min) << " df(x_min)/dt = " << df(x_min) << endl;
+    }
+}
+
+int main(int, char**) {
+    std::cout.precision(15);
+    test_CubicLagrangeMinimize();
+    return 0;
+}

python/CubicLagrangeMinimize.py

@@ -0,0 +1,129 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+import matplotlib.pyplot as plt
+
+def cubic_poly(x, a0, a1, a2, a3):
+    return a0 + a1*x + a2*x**2 + a3*x**3
+
+def polyfit4(x1, x2, x3, x4, y1, y2, y3, y4):
+    A = np.array([[1, x1, x1**2, x1**3], [1, x2, x2**2, x2**3], [1, x3, x3**2, x3**3], [1, x4, x4**2, x4**3]])
+    y = np.array([y1, y2, y3, y4])
+    a = np.linalg.solve(A.transpose()@A, A.transpose()@y)
+    return a
+
+def cubic_lagrange_minimize(f, a, b, tol=1e-6):
+    ''' Function to find minimum of f over interval [a, b] using cubic Lagrange polynomial interpolation.
+        If the function is monotonic, then the minimum is one of the bounds of the interval, and the minimum is found in a single iteration.
+        The best estimate of the minimum is returned.
+        The number of function evaluations is 2 + 2*Niter.
+    '''
+    # initialize interval endpoints and function values
+    x0, x1, x2, x3 = a, a*2./3. + b*1./3., a*1./3. + b*2./3., b
+    f0, f3 = f(x0), f(x3)
+    x_prev = x0
+    delta = 1# only needed to debug print
+    x_sol_1 = 1# only needed to debug print
+    x_sol_2 = 1# only needed to debug print
+    y_sol_1 = 1# only needed to debug print
+    y_sol_2 = 1# only needed to debug print
+
+    reduction_factor = 3# reduction factor of interval size at each iteration is 3 because we sample two points in the interval each iteration
+    Niter = int(np.ceil(np.log(np.abs(b-a)/tol)/np.log(reduction_factor)))
+
+    for i in range(Niter):
+        # Compute function values at two points in the interval
+        x1, x2 = x0*2./3. + x3*1./3., x0*1./3. + x3*2./3.
+        f1, f2 = f(x1), f(x2)
+
+        print('-----------------')
+        print(f'Iteration {i}')
+        print(x0, x1, x2, x3)
+
+        # compute Lagrange polynomial using least-squares fit to 4 points, which is equivalent to the cubic Lagrange polynomial
+        a0, a1, a2, a3 = polyfit4(x0, x1, x2, x3, f0, f1, f2, f3)
+        print('p : ', a0, a1, a2, a3)
+
+        # Solve the first derivative of the Lagrange polynomial for a zero
+        if np.abs(a3) > 1e-9:
+            delta = -3*a1*a3 + a2**2
+
+            if delta < 0:
+                x_sol = x0 if f0 < f3 else x3   # just choose the interval tha contains the minimum of the linear polynomial
+                y_sol = cubic_poly(x_sol, a0, a1, a2, a3)
+            else:                               # solve for the two solutions of the quadratic equation of the first derivative of the Lagrange polynomial
+                x_sol_1 = (-a2 + np.sqrt(delta))/(3*a3)
+                x_sol_2 = (-a2 - np.sqrt(delta))/(3*a3)
+
+                y_sol_1 = cubic_poly(x_sol_1, a0, a1, a2, a3)
+                y_sol_2 = cubic_poly(x_sol_2, a0, a1, a2, a3)
+
+                x_sol = x_sol_1 if y_sol_1 < y_sol_2 else x_sol_2
+                y_sol = y_sol_1 if y_sol_1 < y_sol_2 else y_sol_2
+            
+        elif np.abs(a2) > 1e-9: # if a3 is zero, then the Lagrange polynomial is a quadratic polynomial
+            x_sol = -a1/(2*a2)
+            y_sol = cubic_poly(x_sol, a0, a1, a2, a3)
+        else:   # if a3 and a2 are zero, then the Lagrange polynomial is a linear polynomial
+            x_sol = x0 if f0 < f3 else x3   # just choose the interval tha contains the minimum of the linear polynomial
+            y_sol = cubic_poly(x_sol, a0, a1, a2, a3)
+        
+        print(f'{delta}, f({x_sol_1}) = {y_sol_1} f({x_sol_2}) = {y_sol_2}\t{x_sol}')
+
+        if np.abs(x_sol - x_prev) < tol:
+            break
+
+        # Determine which interval contains the minimum of the cubic polynomial
+        if x_sol < x1:
+            x3, f3 = x1, f1
+        elif x_sol < x2:
+            x0, f0 = x1, f1
+            x3, f3 = x2, f2
+        else:
+            x0, f0 = x2, f2
+        
+        x_prev = x_sol
+    
+    # return best estimate of minimum
+    if y_sol < f0 and y_sol < f3:
+        return x_sol
+    elif f0 < f3:
+        return x0
+    else:
+        return x3
+
+def f(x):
+    return x**2 + np.sin(5*x)
+    # return (x-1)**2
+    # return -x
+    # return np.exp(x)
+    # return -np.exp(-x**2)
+    # return np.exp(-x**2)
+    # return np.ones_like(x)
+
+if __name__ == '__main__':
+    def test_cubic_lagrange_polynomial():
+        x = np.linspace(0, 16, 100)
+        x1, x2, x3, x4 = 1, 3, 7, 15
+        y1, y2, y3, y4 = 1, -2, 3.5, 5
+
+        a0, a1, a2, a3 = polyfit4(x1, x2, x3, x4, y1, y2, y3, y4)
+        print(a0, a1, a2, a3)
+        y = cubic_poly(x, a0, a1, a2, a3)
+
+        # plt.plot(x, cubic_lagrange(x, x1, x2, x3, x4, y1, y2, y3, y4), 'b', label='cubic Lagrange ChatGPT one shot')
+        plt.plot(x, y, 'g', label='cubic Lagrange through coeffs polyfit style baby')
+        plt.plot([x1, x2, x3, x4], [y1, y2, y3, y4], 'ro')
+        plt.grid()
+        plt.show()
+    
+    def test_cubic_lagrange_minimize():
+        x = np.linspace(-2, 2, 100)
+        x_min = cubic_lagrange_minimize(f, -1.2, 1.5, tol=1e-6)
+        print('Solution : ', x_min, f(x_min))
+        plt.plot(x, f(x), 'b')
+        plt.plot(x_min, f(x_min), 'ro')
+        plt.grid()
+        plt.show()
+    
+    test_cubic_lagrange_polynomial()
+    test_cubic_lagrange_minimize()