RustedSciThe 0.4.6

Rust framework for symbolic and numerical computing:BVP ( Newton-Raphson frozen/damped/with collocations ), IVP( BDF, Radau, Backward Euler, LSODE, LSODA, RK45, DoPri), nonlinear equations ( Levenberg, Gavin) and more
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

use super::NR::{NR, solve_linear_system};
use crate::numerical::BVP_Damp::BVP_utils::elapsed_time;
use log::{info, warn};
use nalgebra::{DMatrix, DVector};
use std::time::Instant;
/// CLASSIC LEVENBERG-MARQUARDT ALGORITHM

/// " It can be seen that simple gradient descent and Gauss-Newton iteration are complementary in
/// Levenberg proposed an algorithm based on this observation, whose
/// update rule is a blend of the above mentioned algorithms and is given as
/// xk+1 = xk − (Hk + λI)−1∇f(xk)
/// where H is the Hessian matrix evaluated at xi. This update rule is used as follows. If the error
/// goes down following an update, it implies that our quadratic assumption on f  xЎ is working and
/// we reduce λ (usually by a factor of 10) to reduce the influence of gradient descent. On the other
/// hand, if the error goes up, we would like to follow the gradient more and so λ is increased by the
/// same factor. The Levenberg algorithm is thus -
/// 1. Do an update as directed by the rule above.
/// 2. Evaluate the error at the new parameter vector.
/// 3. If the error has increased as a result the update, then retract the step (i.e. reset the weights to
/// their previous values) and increase λ by a factor of 10 or some such significant factor. Then
/// go to (1) and try an update again.
/// 4. If the error has decreased as a result of the update, then accept the step (i.e. keep the weights
/// at their new values) and decrease λ by a factor of 10 or so. "The Levenberg-Marquardt Algorithm" by Ananth Ranganathan

///

impl NR {
    pub fn step_lm(&mut self) -> (i32, Option<DVector<f64>>) {
        let now = Instant::now();
        let bounds_vec = self.bounds_vec.clone();
        let method = self.linear_sys_method.clone().unwrap();
        let parameters = self.parameters.clone().unwrap();
        let diag_flag: bool = if parameters.get("diag").unwrap() == &1.0 {
            true
        } else {
            false
        };
        let decrease_factor = parameters.get("decrease_factor").unwrap();
        let increase_factor = parameters.get("increase_factor").unwrap();
        let max_lambda = parameters.get("max_lambda").unwrap();
        let min_lambda = parameters.get("min_lambda").unwrap();
        let mut lambda = self.dumping_factor.clone();
        info!("lambda = {}", lambda);
        let y = self.y.clone();
        let Fy = self.evaluate_function(y.clone());
        let Jy = self.evaluate_jacobian(y.clone());
        // Check convergence
        let old_norm = Fy.norm_squared();
        info!("norm = {}", old_norm);
        if old_norm < self.tolerance {
            let y_result = y.clone();
            info!("Solution found!");
            return (1, Some(y_result));
        }
        // Solve (J^T J + λI)p = -J^T F
        let JtJ = Jy.transpose() * &Jy;
        let D: DMatrix<f64> = if diag_flag {
            DMatrix::from_diagonal(&JtJ.diagonal())
        } else {
            DMatrix::identity(y.len(), y.len())
        };
        let A = JtJ + lambda * D * DMatrix::identity(y.len(), y.len());
        let b = -Jy.transpose() * &Fy;
        let step = solve_linear_system(method, &A, &b).unwrap();

        let mut y_new = y.clone() + &step;

        // Project onto bounds (if they exist)
        y_new = self.clip(&y_new, &bounds_vec);

        let Fy_new = self.evaluate_function(y_new.clone());
        // let rho = (old_norm - Fy_new.norm_squared())/old_norm;

        let actual_reduction = Fy.norm_squared() - Fy_new.norm_squared();
        let predicted_reduction = Fy.norm_squared() - (Fy.clone() + Jy * &step).norm_squared();
        let rho = if predicted_reduction.abs() > 1e-12 {
            actual_reduction / (predicted_reduction + 1e-12)
        } else {
            0.0
        };

        info!("rho is: {}", rho);
        let elapsed = now.elapsed();
        elapsed_time(elapsed);
        if rho > 0.0 {
            lambda = (lambda / decrease_factor).max(*min_lambda);
            if lambda == *min_lambda {
                warn!("lambda reached minimum!");
            }
            self.dumping_factor = lambda;
            info!("lambda decreased: {}", lambda);

            return (0, Some(y_new));
        } else {
            lambda = (lambda * increase_factor).min(*max_lambda);
            if lambda == *max_lambda {
                warn!("lambda reached maximum!");
            }
            self.dumping_factor = lambda;
            info!("lambda increased: {}", lambda);
            return (0, Some(y.clone()));
        }
    }
}