Module ceres_solver::nlls_problem

source ·
Expand description

Non-Linear Least Squares problem builder and solver.

The diagram shows the lifecycle of a NllsProblem:

       x
       │ NllsProblem::new()
       │
    ┌──▼────────┐ .solve(self, options) ┌───────────────────┐
┌──►│NllsProblem├──────────────────────►│NllsProblemSolution│
│   └──┬────────┘                       └───────────────────┘
│      │  .residual_block_builder(self)
│   ┌──▼─────────────────┐
│   │ResidualBlockBuilder│
│   └──▲─┬───────────────┘
│      │ │.set_cost(self, func, num_residuals)
│      └─┤
│      ▲ │
│      └─┤.set_loss(self, loss)
│        │
│      ▲ │
│      └─┤.set_parameters(self,
│        │
└────────┘.build_into_problem(self)

We start with NllsProblem with no residual blocks and cannot be solved. Next we should add a residual block by calling NllsProblem::residual_block_builder which is a destructive method which consumes the problem and returns a ResidualBlockBuilder which can be used to build a new residual block. Here we add mandatory cost function crate::cost::CostFunctionType and parameter blocks crate::parameter_block::ParameterBlock. We can also set optional loss function crate::loss::LossFunction. Once we are done, we call ResidualBlockBuilder::build_into_problem which returns previously consumed NllsProblem. Now we can optionally add more residual blocks repeating the process: call NllsProblem::residual_block_builder consuming NllsProblem, add what we need and rebuild the problem. The only difference that now we can re-use parameter blocks used in the previous residual blocks, adding them by their indexes. Once we are done, we can call NllsProblem::solve which consumes the problem, solves it and returns NllsProblemSolution which contains the solution and summary of the solver run. It returns an error if the problem has no residual blocks.

Examples

Multiple residual blocks with shared parameters

Let’s solve a problem of fitting a family of functions y_ij = a + b_i * exp(c_i * x_ij): all of them have the same offset a, but different scale parameters b_i and c_i, i in 0..=k-1 for k (N_CURVES bellow) different sets of data.

use ceres_solver::parameter_block::ParameterBlockOrIndex;
use ceres_solver::{CostFunctionType, NllsProblem, SolverOptions};

// Get parameters, x, y and return tuple of function value and its derivatives
fn target_function(parameters: &[f64; 3], x: f64) -> (f64, [f64; 3]) {
    let [a, b, c] = parameters;
    let y = a + b * f64::exp(c * x);
    let dy_da = 1.0;
    let dy_db = f64::exp(c * x);
    let dy_dc = b * x * f64::exp(c * x);
    (y, [dy_da, dy_db, dy_dc])
}

const N_OBS_PER_CURVE: usize = 100;
const N_CURVES: usize = 3;

// True parameters
let a_true = -2.0;
let b_true: [_; N_CURVES] = [2.0, 2.0, -1.0];
let c_true: [_; N_CURVES] = [3.0, -1.0, 3.0];

// Initial parameter guesses
let a_init = 0.0;
let b_init = 1.0;
let c_init = 1.0;

// Generate data
let x = vec![
    (0..N_OBS_PER_CURVE)
        .map(|i| (i as f64) / (N_OBS_PER_CURVE as f64))
        .collect::<Vec<_>>();
    3
];
let y: Vec<Vec<_>> = x
    .iter()
    .zip(b_true.iter().zip(c_true.iter()))
    .map(|(x, (&b, &c))| {
        x.iter()
            .map(|&x| {
                let (y, _) = target_function(&[a_true, b, c], x);
                // True value + "noise"
                y + 0.001 + f64::sin(1e6 * x)
            })
            .collect()
    })
    .collect();

// Build the problem
let mut problem = NllsProblem::new();
for (i, (x, y)) in x.into_iter().zip(y.into_iter()).enumerate() {
    let cost: CostFunctionType = Box::new(
        move |parameters: &[&[f64]],
              residuals: &mut [f64],
              mut jacobians: Option<&mut [Option<&mut [&mut [f64]]>]>| {
            assert_eq!(parameters.len(), 3);
            let a = parameters[0][0];
            let b = parameters[1][0];
            let c = parameters[2][0];
            // Number of residuls equal to the number of observations
            assert_eq!(residuals.len(), N_OBS_PER_CURVE);
            for (j, (&x, &y)) in x.iter().zip(y.iter()).enumerate() {
                let (y_model, derivatives) = target_function(&[a, b, c], x);
                residuals[j] = y - y_model;
                // jacobians can be None, then you don't need to provide them
                if let Some(jacobians) = jacobians.as_mut() {
                    // The size of the jacobians array is equal to the number of parameters,
                    // each element is Option<&mut [&mut [f64]]>
                    for (mut jacobian, &derivative) in jacobians.iter_mut().zip(&derivatives) {
                        if let Some(jacobian) = &mut jacobian {
                            // Each element in the jacobians array is slice of slices:
                            // the first index is for different residuals components,
                            // the second index is for different components of the parameter vector
                            jacobian[j][0] = -derivative;
                        }
                    }
                }
            }
            true
        },
    );
    let a_parameter: ParameterBlockOrIndex = if i == 0 {
        vec![c_init].into()
    } else {
        0.into()
    };
    problem = problem
        .residual_block_builder()
        .set_cost(cost, N_OBS_PER_CURVE)
        .add_parameter(a_parameter)
        .add_parameter(vec![b_init])
        .add_parameter(vec![c_init])
        .build_into_problem()
        .unwrap()
        .0;
}

// Solve the problem
let solution = problem.solve(&SolverOptions::default()).unwrap();
println!("Brief summary: {:?}", solution.summary);
// Getting parameter values
let a = solution.parameters[0][0];
assert!((a - a_true).abs() < 0.03);
let (b, c): (Vec<_>, Vec<_>) = solution.parameters[1..]
    .chunks(2)
    .map(|sl| (sl[0][0], sl[1][0]))
    .unzip();
for (b, &b_true) in b.iter().zip(b_true.iter()) {
    assert!((b - b_true).abs() < 0.03);
}
for (c, &c_true) in c.iter().zip(c_true.iter()) {
    assert!((c - c_true).abs() < 0.03);
}

Parameter constraints

Let’s find a minimum of the Himmelblau’s function: f(x, y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2 with boundaries x ∈ [0; 3.5], y ∈ [-1.8; 3.5] and initial guess x = 3.45, y = -1.8. This function have four global minima, all having the the same value f(x, y) = 0, one of them is within the boundaries and another one is just outside of them, near the initial guess. The solver converges to the corner of the boundary.

use ceres_solver::{CostFunctionType, NllsProblem, ParameterBlock, SolverOptions};

const LOWER_X: f64 = 0.0;
const UPPER_X: f64 = 3.5;
const LOWER_Y: f64 = -1.8;
const UPPER_Y: f64 = 3.5;

fn solve_himmelblau(initial_x: f64, initial_y: f64) -> (f64, f64) {
    let x_block = {
        let mut block = ParameterBlock::new(vec![initial_x]);
        block.set_all_lower_bounds(vec![LOWER_X]);
        block.set_all_upper_bounds(vec![UPPER_X]);
        block
    };
    let y_block = {
        let mut block = ParameterBlock::new(vec![initial_y]);
        block.set_all_lower_bounds(vec![LOWER_Y]);
        block.set_all_upper_bounds(vec![UPPER_Y]);
        block
    };

    // You can skip type annotations in the closure definition, we use them for verbosity only.
    let cost: CostFunctionType = Box::new(
        move |parameters: &[&[f64]],
              residuals: &mut [f64],
              mut jacobians: Option<&mut [Option<&mut [&mut [f64]]>]>| {
            let x = parameters[0][0];
            let y = parameters[1][0];
            // residuals have the size of your data set, in our case it is two
            residuals[0] = x.powi(2) + y - 11.0;
            residuals[1] = x + y.powi(2) - 7.0;
            // jacobians can be None, then you don't need to provide them
            if let Some(jacobians) = jacobians {
                // The size of the jacobians array is equal to the number of parameters,
                // each element is Option<&mut [&mut [f64]]>
                if let Some(d_dx) = &mut jacobians[0] {
                    // Each element in the jacobians array is slice of slices:
                    // the first index is for different residuals components,
                    // the second index is for different components of the parameter vector
                    d_dx[0][0] = 2.0 * x;
                    d_dx[1][0] = 1.0;
                }
                if let Some(d_dy) = &mut jacobians[1] {
                    d_dy[0][0] = 1.0;
                    d_dy[1][0] = 2.0 * y;
                }
            }
            true
        },
    );

    let solution = NllsProblem::new()
        .residual_block_builder() // create a builder for residual block
        .set_cost(cost, 2) // 2 is the number of residuals
        .set_parameters([x_block, y_block])
        .build_into_problem()
        .unwrap()
        .0 // build_into_problem returns a tuple (NllsProblem, ResidualBlockId)
        .solve(&SolverOptions::default()) // SolverOptions can be customized
        .unwrap(); // Err should happen only if we added no residual blocks

    // Print the full solver report
    println!("{}", solution.summary.full_report());

    (solution.parameters[0][0], solution.parameters[1][0])
}

// The solver converges to the corner of the boundary rectangle.
let (x, y) = solve_himmelblau(3.4, -1.0);
assert_eq!(UPPER_X, x);
assert_eq!(LOWER_Y, y);

// The solver converges to the global minimum inside the boundaries.
let (x, y) = solve_himmelblau(1.0, 1.0);
assert!((3.0 - x).abs() < 1e-8);
assert!((2.0 - y).abs() < 1e-8);

Structs