Skip to main content

gmac/morph/
linear_algebra.rs

1use hologram::linear_algebra::lu_linear_solver;
2
3use crate::error::Result;
4
5/// Solves a system of linear equations using least squares.
6///
7/// # Arguments
8/// * `mat`: The design matrix.
9/// * `rhs`: A vector of n-dimensional points as the right-hand side.
10///
11/// # Returns
12/// * `Ok(Vec<Vec<f64>>)`: A Result containing either:
13/// A vector of transformed points (`Ok`)
14/// An error if something goes wrong (`Err`)
15pub fn least_squares_solver(mat: &[Vec<f64>], rhs: &[Vec<f64>]) -> Result<Vec<Vec<f64>>> {
16    let mat_rows = mat.len();
17    let mat_cols = mat[0].len();
18    let rhs_dim = rhs[0].len();
19
20    // Calculate A^T A
21    let mut ata = vec![vec![0.0; mat_cols]; mat_cols];
22    for i in 0..mat_cols {
23        for j in 0..mat_cols {
24            for k in 0..mat_rows {
25                ata[i][j] += mat[k][i] * mat[k][j];
26            }
27        }
28    }
29
30    // Calculate A^T b for each dimension
31    let mut atb = vec![vec![0.0; rhs_dim]; mat_cols];
32    for i in 0..mat_cols {
33        for d in 0..rhs_dim {
34            for k in 0..mat_rows {
35                atb[i][d] += mat[k][i] * rhs[k][d];
36            }
37        }
38    }
39
40    let x = lu_linear_solver(&ata, &atb)?;
41
42    Ok(x)
43}