commit e44c3a36a92a3639f53dde7d3d0be1c5fa5321a7 Author: Jerome Date: Fri Mar 17 20:29:20 2023 +0100 Initial commit. Python and C++ version using Eigen. diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4caddcf --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +cpp/build/* +python/__pycache__/* \ No newline at end of file diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt new file mode 100644 index 0000000..3d13685 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/cpp/Makefile b/cpp/Makefile new file mode 100644 index 0000000..fb1c189 --- /dev/null +++ b/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) diff --git a/cpp/include/CubicLagrangeMinimize.hpp b/cpp/include/CubicLagrangeMinimize.hpp new file mode 100644 index 0000000..941246f --- /dev/null +++ b/cpp/include/CubicLagrangeMinimize.hpp @@ -0,0 +1,24 @@ +#ifndef DEF_CubicLagrangeMinimize +#define DEF_CubicLagrangeMinimize + +#include +#include + +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 f, double a, double b, double tol=1e-9); + +#endif \ No newline at end of file diff --git a/cpp/src/CubicLagrangeMinimize.cpp b/cpp/src/CubicLagrangeMinimize.cpp new file mode 100644 index 0000000..9a5b9cf --- /dev/null +++ b/cpp/src/CubicLagrangeMinimize.cpp @@ -0,0 +1,107 @@ +#include + +/// @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 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(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; +} diff --git a/cpp/src/tests.cpp b/cpp/src/tests.cpp new file mode 100644 index 0000000..09022f8 --- /dev/null +++ b/cpp/src/tests.cpp @@ -0,0 +1,80 @@ +#include +#include +#include + +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> fcts = {fct_01, fct_02, fct_03, fct_04, fct_05, fct_06, fct_07}; + std::vector> dfcts = {dfct_01, dfct_02, dfct_03, dfct_04, dfct_05, dfct_06, dfct_07}; + std::vector mins = {-1.2, -10., -1.5, -10., -1., -2., -3.}; + std::vector 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; +} diff --git a/python/CubicLagrangeMinimize.py b/python/CubicLagrangeMinimize.py new file mode 100644 index 0000000..ed4b71b --- /dev/null +++ b/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() \ No newline at end of file