jerome/rust-numerical-solvers
1
0
Fork 0
Code Issues 2 Pull requests Projects Releases Packages Wiki Activity

Added Nelder-Mead optimization algorithm.

Browse source
This commit is contained in:
Jérôme 2023-03-27 22:46:06 +02:00
parent 6134559471
commit 6be562ae46
5 changed files with 257 additions and 17 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

14
.vscode/launch.json vendored Normal file
View file

@ -0,0 +1,14 @@
{
// Use IntelliSense to learn about possible attributes.
// Hover to view descriptions of existing attributes.
// For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
"version": "0.2.0",
"configurations": [{
"name": "Debug launch",
"type": "lldb",
"request": "launch",
"program": "${workspaceRoot}/target/debug/rust_numerical_solvers",
"args": [],
"cwd": "${workspaceRoot}",
}]
}

3
Cargo.toml
View file

@ -6,4 +6,5 @@ edition = "2021"
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[dependencies]
colored = "2.0.0"
colored = "2.0.0"
nalgebra = "0.32.2"

78
src/main.rs
View file

@ -1,7 +1,15 @@
extern crate colored;
use colored::Colorize;
extern crate nalgebra as na;
mod univariate_solvers;
mod univariate_minimizers;
mod nelder_mead;
use colored::Colorize;
fn rosenbrock(x: &na::DVector<f64>) -> f64 {
return (1.0-x[0]).powi(2) + 100.0*(x[1] - x[0].powi(2)).powi(2);
}
/// x sin(x)
/// e + ──────
@ -52,6 +60,27 @@ fn check_result(x : f64, x_true : f64, tol : f64, test_name: &str, verbose: bool
}
}
fn check_result_optim(x: &na::DVector<f64>, f_x: f64, x_true: &na::DVector<f64>, f_x_true: f64, tol_x: f64, tol_f_x: f64, test_name: &str, verbose: bool) -> u32 {
let diff_x : f64 = (x - x_true).norm();
let diff_f_x : f64 = (f_x - f_x_true).abs();
let test_name_padded: String = format!("{:<30}", test_name);
if diff_x < tol_x && diff_f_x < tol_f_x {
if verbose {
println!("{}\t: x = {}\tf(x) = {}\t{}", test_name_padded, format!("{:<20}", x), format!("{:<40}", f_x), "passed".green());
} else {
println!("{} {}", test_name_padded, "passed".green());
}
return 1;
} else {
if verbose {
println!("{}\t: x = {}\tf(x) = {}\t{} (expected {}, delta = {})", test_name_padded, format!("{:<20}", x), format!("{:<40}", f_x), "failed".red(), x_true, (x - x_true));
} else {
println!("{} {} : expected {}, got {}", test_name_padded, "failed".red(), x_true, x);
}
return 0;
}
}
fn print_test_results(num_tests_passed: u32, num_tests_total: u32) {
let ratio_str:String = format!("{}/{} ({} %)", num_tests_passed, num_tests_total, ((num_tests_passed as f64)/(num_tests_total as f64)*100.0).round());
if num_tests_passed == num_tests_total {
@ -63,23 +92,23 @@ fn print_test_results(num_tests_passed: u32, num_tests_total: u32) {
fn test_univariate_solvers(verbose: bool) {
println!("Testing univariate numerical solvers.");
let x0 : f64 = 1.0;
let tol : f64 = 1e-10;
let x0 : f64 = 1.0;
let tol : f64 = 1e-10;
let max_iter : u32 = 100;
let dx_num : f64 = 1e-6;
let dx_num : f64 = 1e-6;
let mut num_tests_passed : u32 = 0;
let num_tests_total : u32 = 7;
let num_tests_total : u32 = 7;
let x_mathematica: f64 = -3.26650043678562449167148755288;// 30 digits of precision
let x_mathematica: f64 = -3.26650043678562449167148755288;// 30 digits of precision
let x_mathematica_2: f64 = -6.27133405258685307845641527902;// 30 digits of precision
// let x_mathematica_3: f64 = -9.42553801930504142668603949182;// 30 digits of precision
let x_newton: f64 = univariate_solvers::newton_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), x0, tol, max_iter);
let x_newton: f64 = univariate_solvers::newton_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), x0, tol, max_iter);
let x_newton_num: f64 = univariate_solvers::newton_solve_num(&(fct as fn(f64) -> f64), x0, tol, dx_num, max_iter);
let x_halley: f64 = univariate_solvers::halley_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), &(ddfct as fn(f64) -> f64), x0, tol, max_iter, false).unwrap();
let x_halley: f64 = univariate_solvers::halley_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), &(ddfct as fn(f64) -> f64), x0, tol, max_iter, false).unwrap();
let x_halley_num: f64 = univariate_solvers::halley_solve_num(&(fct as fn(f64) -> f64), x0, tol, dx_num, max_iter, false).unwrap();
let x_bisection : f64 = univariate_solvers::bisection_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol).unwrap();
let x_secant : f64 = univariate_solvers::secant_solve(&(fct as fn(f64) -> f64), -1.0, 1.0, tol, max_iter);
let x_ridder : f64 = univariate_solvers::ridder_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol, max_iter).unwrap();
let x_secant : f64 = univariate_solvers::secant_solve(&(fct as fn(f64) -> f64), -1.0, 1.0, tol, max_iter);
let x_ridder : f64 = univariate_solvers::ridder_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol, max_iter).unwrap();
num_tests_passed += check_result(x_newton, x_mathematica, tol, "Newton's method", verbose);
num_tests_passed += check_result(x_newton_num, x_mathematica, tol, "Newton's method (num)", verbose);
num_tests_passed += check_result(x_halley, x_mathematica_2, tol, "Halley's method", verbose);
@ -92,22 +121,41 @@ fn test_univariate_solvers(verbose: bool) {
}
fn test_univariate_optimizers(verbose: bool) {
let tol : f64 = 1e-10;
let tol : f64 = 1e-10;
let max_iter : u32 = 100;
let dx_num : f64 = 1e-6;
let dx_num : f64 = 1e-6;
let mut num_tests_passed : u32 = 0;
let num_tests_total : u32 = 1;
let num_tests_total : u32 = 1;
let x_mathematica: f64 = -4.54295618675514754103476876324;// 30 digits of precision
let y_mathematica: f64 = -0.206327079359226884630654987440;// 30 digits of precision
let x_mathematica: f64 = -4.54295618675514754103476876324;// 30 digits of precision
let y_mathematica: f64 = -0.206327079359226884630654987440;// 30 digits of precision
let x_golden_section: f64 = univariate_minimizers::golden_section_minimize(&(fct as fn(f64) -> f64), -7.0, -1.0, tol);
num_tests_passed += check_result(x_golden_section, x_mathematica, tol*1e2, "Golden section search", verbose);
print_test_results(num_tests_passed, num_tests_total);
}
fn test_multivariate_optimizers(verbose: bool) {
let tol : f64 = 1e-10;
let max_iter : u32 = 1000;
let dx_num : f64 = 1e-6;
let mut num_tests_passed : u32 = 0;
let num_tests_total : u32 = 1;
let tol_x: f64 = 1e-4;
let tol_f_x: f64 = 1e-6;
let x_true: na::DVector<f64> = na::DVector::from_vec(vec![1.,1.]);
let f_x_true: f64 = rosenbrock(&x_true);
let sol_nelder_mead: (na::DVector<f64>, f64) = nelder_mead::nelder_mead(&(rosenbrock as fn(&na::DVector<f64>) -> f64), &na::DVector::from_vec(vec![2.0,-1.0]), 0.1, tol, max_iter, false);
// num_tests_passed += check_result(x_nelder_mead, x_true, tol*1e2, "Nelder-Mead", verbose);
num_tests_passed += check_result_optim(&sol_nelder_mead.0, sol_nelder_mead.1, &x_true, f_x_true, tol_x, tol_f_x, "Nelder-Mead", verbose);
print_test_results(num_tests_passed, num_tests_total);
}
fn main() {
println!("Testing Rust numerical solvers.");
let verbose : bool = true;
test_univariate_solvers(verbose);
test_univariate_optimizers(verbose);
test_multivariate_optimizers(verbose);
}

177
src/nelder_mead.rs Normal file
View file

@ -0,0 +1,177 @@
extern crate nalgebra as na;
// let mut tuple_list2: Vec<(u16, u16)> = vec![(1, 5), (0, 17), (8, 2)];
// tuple_list2.sort_by(|a, b| a.1.cmp(&b.1));
/// Computes the mean of the points.
/// @param points The points.
/// @return The mean of the points.
fn mean_of_points(points: &Vec<(na::DVector<f64>, f64)>) -> f64 {
let mut sum: f64 = 0.0;
for point in points {
sum += point.1;
}
return sum / points.len() as f64;
}
/// Computes the standard deviation of the points.
/// @param points The points.
/// @return The standard deviation of the points.
fn standard_deviation_of_points(points: &Vec<(na::DVector<f64>, f64)>) -> f64 {
let mean: f64 = mean_of_points(points);
let mut sum: f64 = 0.0;
for value in points {
sum += (value.1 - mean).powi(2);
}
return f64::sqrt(sum / points.len() as f64);
}
/// Computes the average distance between each point in the simplex.
/// @param simplex The simplex.
/// @return The average distance between each point in the simplex.
fn compute_simplex_size(simplex: &Vec<(na::DVector<f64>, f64)>) -> f64 {
let mut sum: f64 = 0.0;
for i in 0..simplex.len() {
for j in 0..simplex.len() {
if &i != &j {
sum += (&simplex[i].0 - &simplex[j].0).norm();
}
}
}
return sum / (simplex.len() as f64 * simplex.len() as f64);
}
/// Implements the Nelder-Mead gradient-less optimization algorithm.
/// @param f The function to optimize.
/// @param x0 The starting point of the algorithm.
/// @param simplex_size The maximum size of the simplex at initialization.
pub fn nelder_mead<F>(f: F, x0: &na::DVector<f64>, simplex_size: f64, tol: f64, max_iter: u32, verbose: bool) -> (na::DVector<f64>, f64)
where F : Fn(&na::DVector<f64>) -> f64
{
// Parameters
let alpha: f64 = 1.0; // Reflection coefficient
let gamma: f64 = 2.0; // Expansion coefficient
let rho: f64 = 0.5; // Contraction coefficient
let sigma: f64 = 0.5; // Shrink coefficient
// Create simplex around x0
let mut simplex: Vec<(na::DVector<f64>, f64)> = Vec::new();
simplex.push((x0.clone(), f(&x0)));// Initial point of the simplex = starting point of the algorithm
for i in 0..x0.len() {
simplex.push((x0.clone(), 0.0));
simplex[i+1].0[i] += simplex_size; // Initialise each point in the simplex to be simplex_size away from x0 along aech direction
simplex[i+1].1 = f(&simplex[i+1].0);// Evaluate objective function at each point in the simplex
}
let N: &usize = &simplex.len(); // Number of points in the simplex
if verbose {
println!("Initial simplex:");
for v in &simplex {
println!("{:?}", v);
}
}
for iter in 0..max_iter {
// Sort the simplex by function value
simplex.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
if verbose {
println!("\nIteration {}\nSorted simplex:", iter);
for v in &simplex {
println!("{} -> {}", v.0, v.1);
}
}
// Termination condition 1 : convergence of the simplex values
let std_fvals: f64 = standard_deviation_of_points(&simplex);
if std_fvals < tol {
if verbose {
println!("Converged on function values after {} iterations", iter);
}
return simplex[0].clone();
}
// Termination condition 2 : convergence of the simplex size
let average_simplex_size: f64 = compute_simplex_size(&simplex);
if average_simplex_size < tol {
if verbose {
println!("Converged on simplex size after {} iterations", iter);
}
return simplex[0].clone();
}
// Compute the centroid of the simplex (excluding the worst point)
let mut centroid: na::DVector<f64> = na::DVector::zeros(x0.len());
for i in 0..(N-1) {
centroid += &simplex[i].0;
}
centroid /= (N-1) as f64;
if verbose {
println!("Centroid: {}", centroid);
}
// Compute the reflection point : xr = x0 + alpha*(x0 - x[N+1])
let reflection: na::DVector<f64> = &centroid + alpha*(&centroid - &simplex[N-1].0);
let f_reflection: f64 = f(&reflection);
if verbose {
println!("Reflection: {} -> {}", reflection, f_reflection);
}
// If the reflection point is better than the best point, replace the worst point with the reflection point
if &simplex[0].1 <= &f_reflection && &f_reflection < &simplex[N-2].1 {
simplex[N-1] = (reflection, f_reflection);
continue;// Go to next iteration
}
// If the reflection point is better than the best current point, try an expansion
if &f_reflection < &simplex[0].1 {
let expansion: na::DVector<f64> = &centroid + gamma*(&reflection - &centroid);
let f_expansion: f64 = f(&expansion);
if verbose {
println!("Expansion: {} -> {}", expansion, f_expansion);
}
if &f_expansion <= &f_reflection {
simplex[N-1] = (expansion, f_expansion);
} else {
simplex[N-1] = (reflection, f_reflection);
}
continue;// Go to next iteration
}
// If the reflection point is worse than the second worst point, try a contraction
if &f_reflection >= &simplex[N-2].1 {
let contraction: na::DVector<f64> = &centroid + rho*(&simplex[N-1].0 - &centroid);
let f_contraction: f64 = f(&contraction);
if verbose {
println!("Contraction: {} -> {}", contraction, f_contraction);
}
if &f_contraction < &simplex[N-1].1 {
simplex[N-1] = (contraction, f_contraction);
continue;// Go to next iteration
} else {
// Shrink the simplex
for i in 1..*N {
simplex[i].0 = &simplex[0].0 + sigma*(&simplex[i].0 - &simplex[0].0);
simplex[i].1 = f(&simplex[i].0);
}
if verbose {
println!("Shrinking the whole simplex");
}
}
}
}
if verbose {
println!("Maximum number of iterations reached");
}
return simplex[0].clone();
}

2
src/univariate_minimizers.rs
View file

@ -8,7 +8,7 @@
pub fn golden_section_minimize<F>(f : F, mut a: f64, mut b: f64, tol: f64) -> f64
where F : Fn(f64) -> f64
{
let invphi: f64 = (f64::sqrt(5.0) - 1.0) / 2.0; // 1 / phi
let invphi: f64 = (f64::sqrt(5.0) - 1.0) / 2.0; // 1 / phi
let invphi2: f64 = (3.0 - f64::sqrt(5.0)) / 2.0; // 1 / phi^2
a = f64::min(a, b);