Unexpected results using `faer_core::solve::solve_lower_triangular_in_place_with_conj`
Closed this issue · 1 comments
DJDuque commented
First of all, thanks a lot for your work in faer
, it is impressive. It has also been a lot of fun to play around with the library.
Describe the bug
The results of solving a lower triangular system are unexpected.
To Reproduce
The following code reproduces the error. It also shows how the equivalent code with nalgebra
gets the correct solution:
# Cargo.toml
[package]
name = "test"
version = "0.1.0"
edition = "2021"
[dependencies]
faer-core = "0.9.1"
nalgebra = "0.32.2"
// main.rs
const RESPONSE: [f64; 3] = [-18.0, -60.0, -80.0];
// Toeplitz and lower triangular
fn r_element(i: usize, j: usize) -> f64 {
if i >= j {
let diff = i - j;
if diff < RESPONSE.len() {
RESPONSE[diff]
} else {
0.0
}
} else {
0.0
}
}
fn r_matrix_faer(n: usize) -> faer_core::Mat<f64> {
faer_core::Mat::with_dims(n, n, |i, j| r_element(i, j))
}
fn r_matrix_nalgebra(n: usize) -> nalgebra::DMatrix<f64> {
nalgebra::DMatrix::from_fn(n, n, |i, j| r_element(i, j))
}
// Solve for X in Y = R * X
fn solve_x_faer(r: &faer_core::Mat<f64>, y: &mut faer_core::Mat<f64>) {
faer_core::solve::solve_lower_triangular_in_place_with_conj(
r.as_ref(),
faer_core::Conj::No,
y.as_mut(),
faer_core::Parallelism::None,
);
}
// Solve for X in Y = R * X
fn solve_x_nalgebra(r: &nalgebra::DMatrix<f64>, y: &mut nalgebra::DMatrix<f64>) {
r.solve_lower_triangular_mut(y);
}
fn main() {
let (i, j) = (400, 1);
let r = r_matrix_faer(i);
let mut x = faer_core::Mat::zeros(i, j);
x.write(0, 0, 150.0);
let mut y = r.clone() * x;
solve_x_faer(&r, &mut y); // The values in `y` blow up
// Now, the same thing with nalgebra
let r = r_matrix_nalgebra(i);
let mut x = nalgebra::DMatrix::zeros(i, j);
x[(0, 0)] = 150.0;
let mut y = r.clone() * x;
solve_x_nalgebra(&r, &mut y); // The values in `y` are correct
}
Expected behavior
I would expect the same results as nalgebra
sarah-quinones commented
can't fix this, as it turns out to be due to numerical errors arising from the very high condition number of the r
matrix