lstsq/
lib.rs

1#![cfg_attr(not(feature = "std"), no_std)]
2
3//! Return the least-squares solution to a linear matrix equation
4//!
5//! The crate implements the linear least squares solution to a linear matrix
6//! equation.
7//!
8//! Characteristics:
9//!
10//! * Linear algebra and types from the [`nalgebra`](https://docs.rs/nalgebra)
11//!   crate.
12//! * Maximum compatibility with the
13//!   [`numpy.linalg.lstsq`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html)
14//!   Python library function.
15//! * No standard library is required (disable the default features to disable
16//!   use of `std`) and no heap allocations. In other words, this can run on a
17//!   bare-metal microcontroller with no OS.
18//!
19//! Example:
20//!
21//! ```rust
22//! use nalgebra::{self as na, OMatrix, OVector, U2};
23//!
24//! let a = OMatrix::<f64, na::Dyn, U2>::from_row_slice(&[
25//!     1.0, 1.0,
26//!     2.0, 1.0,
27//!     3.0, 1.0,
28//!     4.0, 1.0,
29//! ]);
30//!
31//! let b = OVector::<f64, na::Dyn>::from_row_slice(&[2.5, 4.4, 6.6, 8.5]);
32//!
33//! let epsilon = 1e-14;
34//! let results = lstsq::lstsq(&a, &b, epsilon).unwrap();
35//!
36//! assert_eq!(results.solution.nrows(), 2);
37//! approx::assert_relative_eq!(results.solution[0], 2.02, epsilon = epsilon);
38//! approx::assert_relative_eq!(results.solution[1], 0.45, epsilon = epsilon);
39//! approx::assert_relative_eq!(results.residuals, 0.018, epsilon = epsilon);
40//! assert_eq!(results.rank, 2);
41//! ```
42
43use nalgebra::allocator::Allocator;
44use nalgebra::base::{OMatrix, OVector};
45use nalgebra::dimension::{Dim, DimDiff, DimMin, DimMinimum, DimSub, U1};
46use nalgebra::{DefaultAllocator, RealField};
47
48/// Results of [lstsq]
49pub struct Lstsq<R: RealField, N: Dim>
50where
51    DefaultAllocator: Allocator<N>,
52{
53    /// Least-squares solution.
54    ///
55    /// This is the variable `x` that approximatively solves the equation `a * x = b`.
56    pub solution: OVector<R, N>,
57    /// Sums of squared residuals: Squared Euclidean 2-norm.
58    pub residuals: R,
59    /// Rank of matrix `a`.
60    pub rank: usize,
61}
62
63/// Return the least-squares solution to a linear matrix equation.
64///
65/// Computes the vector x that approximatively solves the equation `a * x = b`.
66/// Usage is maximally compatible with Python's `numpy.linalg.lstsq`.
67///
68/// Arguments:
69///
70/// - `a`: "Coefficient" matrix (shape: M rows, N columns)
71/// - `b`: Ordinate or “dependent variable” values (shape: M dimensional)
72/// - `epsilon`: singular values less than this are assumed to be zero.
73///
74/// Returns:
75/// - `Result<`[Lstsq]`,&'static str>`
76///
77/// See the module level documentation for example of usage.
78pub fn lstsq<R, M, N>(
79    a: &OMatrix<R, M, N>,
80    b: &OVector<R, M>,
81    epsilon: R,
82) -> Result<Lstsq<R, N>, &'static str>
83where
84    R: RealField,
85    M: DimMin<N>,
86    N: Dim,
87    DimMinimum<M, N>: DimSub<U1>, // for Bidiagonal.
88    DefaultAllocator: Allocator<M, N>
89        + Allocator<N>
90        + Allocator<M>
91        + Allocator<DimDiff<DimMinimum<M, N>, U1>>
92        + Allocator<DimMinimum<M, N>, N>
93        + Allocator<M, DimMinimum<M, N>>
94        + Allocator<DimMinimum<M, N>>,
95{
96    // calculate solution with epsilon
97    let svd = nalgebra::linalg::SVD::new(a.clone(), true, true);
98    let solution = svd.solve(b, epsilon.clone())?;
99
100    // calculate residuals
101    let model: OVector<R, M> = a * &solution;
102    let l1: OVector<R, M> = model - b;
103    let residuals: R = l1.dot(&l1);
104
105    // calculate rank with epsilon
106    let rank = svd.rank(epsilon);
107
108    Ok(Lstsq {
109        solution,
110        residuals,
111        rank,
112    })
113}
114
115#[cfg(test)]
116mod tests {
117    use crate::lstsq;
118
119    use na::{OMatrix, OVector, RealField, U2};
120    use nalgebra as na;
121
122    fn check_residuals<R: RealField + Copy>(epsilon: R) {
123        /*
124        import numpy as np
125        A = np.array([[1.0, 1.0], [2.0, 1.0], [3.0, 1.0], [4.0, 1.0]])
126        b = np.array([2.5, 4.4, 6.6, 8.5])
127        x,residuals,rank,s = np.linalg.lstsq(A,b)
128        */
129        let a: Vec<R> = vec![1.0, 1.0, 2.0, 1.0, 3.0, 1.0, 4.0, 1.0]
130            .into_iter()
131            .map(na::convert)
132            .collect();
133
134        let a = OMatrix::<R, na::Dyn, U2>::from_row_slice(&a);
135
136        let b_data: Vec<R> = vec![2.5, 4.4, 6.6, 8.5]
137            .into_iter()
138            .map(na::convert)
139            .collect();
140        let b = OVector::<R, na::Dyn>::from_row_slice(&b_data);
141
142        let results = lstsq(&a, &b, R::default_epsilon()).unwrap();
143        assert_eq!(results.solution.nrows(), 2);
144        approx::assert_relative_eq!(results.solution[0], na::convert(2.02), epsilon = epsilon);
145        approx::assert_relative_eq!(results.solution[1], na::convert(0.45), epsilon = epsilon);
146        approx::assert_relative_eq!(results.residuals, na::convert(0.018), epsilon = epsilon);
147        assert_eq!(results.rank, 2);
148    }
149
150    #[test]
151    fn test_residuals_f64() {
152        check_residuals::<f64>(1e-14)
153    }
154
155    #[test]
156    fn test_residuals_f32() {
157        check_residuals::<f32>(1e-5)
158    }
159}
lstsq/lib.rs

lstsq/
lib.rs
