gomez 0.3.1

A pure Rust framework and implementation of (derivative-free) methods for solving nonlinear (bound-constrained) systems of equations.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

//! Analysis on initial guesses for Newton-based algorithms.
//!
//! Based on \[1\]. The analysis provide insights for focusing on which
//! variables need to be improved regarding initial guess for Newton-based
//! algorithms (e.g., [`TrustRegion`](crate::solver::trust_region)).
//!
//! \[1\] [On the Choice of Initial Guesses for the Newton-Raphson
//! Algorithm](https://arxiv.org/abs/1911.12433)

use std::marker::PhantomData;

use nalgebra::{
    allocator::Allocator, convert, storage::StorageMut, ComplexField, DefaultAllocator, DimMin,
    DimName, IsContiguous, OVector, Vector, U1,
};
use thiserror::Error;

use crate::{
    core::{Domain, Problem, ProblemError, System},
    derivatives::{Jacobian, JacobianError, EPSILON_SQRT},
};

/// Error returned from [`InitialGuessAnalysis`] solver.
#[derive(Debug, Error)]
pub enum InitialGuessAnalysisError {
    /// Error that occurred when evaluating the system.
    #[error("{0}")]
    System(#[from] ProblemError),
    /// Error that occurred when computing the Jacobian matrix.
    #[error("{0}")]
    Jacobian(#[from] JacobianError),
}

/// Initial guesses analyzer. See [module](self) documentation for more details.
pub struct InitialGuessAnalysis<F: Problem> {
    nonlinear: Vec<usize>,
    ty: PhantomData<F>,
}

impl<F: System> InitialGuessAnalysis<F>
where
    DefaultAllocator: Allocator<F::Scalar, F::Dim>,
    DefaultAllocator: Allocator<F::Scalar, F::Dim, F::Dim>,
    F::Dim: DimMin<F::Dim, Output = F::Dim>,
    DefaultAllocator: Allocator<F::Scalar, <F::Dim as DimMin<F::Dim>>::Output>,
{
    /// Analyze the system in given point.
    pub fn analyze<Sx, Sfx>(
        f: &F,
        dom: &Domain<F::Scalar>,
        x: &mut Vector<F::Scalar, F::Dim, Sx>,
        fx: &mut Vector<F::Scalar, F::Dim, Sfx>,
    ) -> Result<Self, InitialGuessAnalysisError>
    where
        Sx: StorageMut<F::Scalar, F::Dim> + IsContiguous,
        Sfx: StorageMut<F::Scalar, F::Dim>,
    {
        let scale_iter = dom.vars().iter().map(|var| var.scale());
        let scale = OVector::from_iterator_generic(f.dim(), U1::name(), scale_iter);

        // Compute F'(x) in the initial point.
        f.eval(x, fx)?;
        let jac1 = Jacobian::new(f, x, &scale, fx)?;

        // Compute Newton step.
        let mut p = fx.clone_owned();
        p.neg_mut();

        let qr = jac1.clone_owned().qr();
        qr.solve_mut(&mut p);

        // Do Newton step.
        p *= convert::<_, F::Scalar>(0.001);
        *x += p;

        // Compute F'(x) after one Newton step.
        f.eval(x, fx)?;
        let jac2 = Jacobian::new(f, x, &scale, fx)?;

        // Linear variables have no effect on the Jacobian matrix. They can be
        // recognized by observing no change in corresponding columns (i.e.,
        // columns are constant) in two different points. Being constant is
        // checked with tolerance of sqrt(eps) because a change of such
        // magnitude is caused by finite difference approach and does not
        // correspond to the analytic reality. We are interested only in
        // nonlinear variables that have influence on Jacobian matrix.
        let nonlinear = jac1
            .column_iter()
            .zip(jac2.column_iter())
            .enumerate()
            .filter(|(_, (c1, c2))| {
                c1.iter()
                    .zip(c2.iter())
                    .any(|(a, b)| (*a - *b).abs() > convert(EPSILON_SQRT))
            })
            .map(|(col, _)| col)
            .collect();

        Ok(Self {
            nonlinear,
            ty: PhantomData,
        })
    }

    /// Returns indices of variables that have influence on the Jacobian matrix
    /// of the system.
    ///
    /// Linear variables do not have influence on the Jacobian matrix (the
    /// corresponding columns remain constant during the solving process), thus
    /// they are not important for Newton-based algorithms. Thus, one should
    /// focus on initial guesses for the variables returned by this getter.
    pub fn nonlinear(&self) -> &[usize] {
        &self.nonlinear
    }
}