Closed termoshtt closed 5 years ago
If it helps, this is an instance problem that always fails. Copied from the cholesky test, but without randomness:
use ndarray::{arr1,arr2};
#[test]
fn cholesky_solve_fail() {
let a: Array2<f64> = arr2(&[[2.7518138512052546, 0.0, 0.0],
[-3.0030534432813334, 0.5503895871967794, 0.0],
[0.0, 0.0, 0.4472135954999579]]);
let x: Array1<f64> = arr1(&[-0.3633966736383765772089, -9.2503560214093027980198, 0.0]);
let b = a.dot(&x);
println!("a = \n{:#?}", a);
println!("x = \n{:#?}", x);
println!("b = \n{:#?}", b);
println!("Solution = \n{:#?}", &a.solvec(&b).unwrap());
assert_close_l2!(&a.solvec(&b).unwrap(), &x, 1e-9);
assert_close_l2!(&a.solvec_into(b.clone()).unwrap(), &x, 1e-9);
assert_close_l2!(&a.solvec_inplace(&mut b.clone()).unwrap(), &x, 1e-9);
assert_close_l2!(&a.factorizec(UPLO::Upper).unwrap().solvec(&b).unwrap(), &x, 1e-9);
assert_close_l2!(&a.factorizec(UPLO::Lower).unwrap().solvec(&b).unwrap(), &x, 1e-9);
assert_close_l2!(&a.factorizec(UPLO::Upper).unwrap().solvec_into(b.clone()).unwrap(), &x, 1e-9);
assert_close_l2!(&a.factorizec(UPLO::Lower).unwrap().solvec_into(b.clone()).unwrap(), &x, 1e-9);
assert_close_l2!(&a.factorizec(UPLO::Upper).unwrap().solvec_inplace(&mut b.clone()).unwrap(), &x, 1e-9);
assert_close_l2!(&a.factorizec(UPLO::Lower).unwrap().solvec_inplace(&mut b.clone()).unwrap(), &x, 1e-9);
}
Output:
---- cholesky_solve_fail stdout ----
a =
[[2.7518138512052546, 0.0, 0.0],
[-3.0030534432813334, 0.5503895871967794, 0.0],
[0.0, 0.0, 0.4472135954999579]] shape=[3, 3], strides=[3, 1], layout=C (0x1)
x =
[-0.3633966736383766, -9.250356021409303, 0.0] shape=[3], strides=[1], layout=C | F (0x3)
b =
[-1.0, -4.0, 0.0] shape=[3], strides=[1], layout=C | F (0x3)
Solution =
[-0.36339667363837663, -7.26757935296819, 0.0] shape=[3], strides=[1], layout=C | F (0x3)
thread 'cholesky_solve_fail' panicked at 'called `Result::unwrap()` on an `Err` value: 0.21418077766308938', libcore/result.rs:945:5
note: Run with `RUST_BACKTRACE=1` for a backtrace.
First (and third) component are OK though.
@davenza Thanks for taking a look at this. This particular example fails because a
is asymmetric. (Solving with Cholesky decomposition requires a Hermitian (or real symmetric) positive definite coefficient matrix.) Did you get that a
matrix from running the cholesky_solve
test (which calls random_hpd
to generate a random Hermitian positive-definite matrix), or did you create it yourself? In other words, is random_hpd
the real problem here?
This does bring up another question: In this particular example, a bad value is returned since a
is asymmetric. Should we explicitly check if a
is Hermitian and positive definite, and if not return an Err
?
@jturner314 Thanks, I didn't notice the symmetric requirement. I was trying to do a different thing, and I think I misunderstood the documentation.
I had a covariance matrix (not created with random_hpd), let's call U
:
[7.572479471685095, -8.263844061131207, 0.0] [-8.263844061131206, 9.321258680898515, 0.0] [0.0, 0.0, 0.2]
I found the lower cholesky decomposition U = LL^T
. L
is the matrix that I used in the code above:
[2.7518138512052546, 0.0, 0.0] [-3.0030534432813334, 0.5503895871967794, 0.0] [0.0, 0.0, 0.4472135954999579]
Then, I tried to solve Lx = [1 4 0] using solvec.
However, related with the problem: while I was doing some cargo test runs, I sometimes had errors. I could try to reproduce the errors again, printing a failing matrix if it is interesting to you.
Fwiw, if you need both the triangular factor (L
matrix) and a solution (with .solvec()
), you can factor once with .factorizec()
or .factorizec_into()
to get a CholeskyFactorized
instance. You can then solve, take the inverse, take the determinant, get the lower/upper triangular factor, etc., using the CholeskyFactorized
instance without having to factor the matrix multiple times.
Regarding the test failures: I think they're due to floating-point errors. (In the cases I've observed the tests failing, the errors have always been pretty small. If you see a case where the error is large, please post it here.) One idea to minimize issues associated with floating-point errors in the tests is to throw out randomly-generated matrices that are ill-conditioned (and generate new ones that are well-conditioned), but I haven't had a chance to try this.
In two of my executions of the cholesky_solve() test, the test failed with these matrices:
It seems that both failures happen with the f32 type and due to floating-point errors, as you said.
One idea to minimize issues associated with floating-point errors in the tests is to throw out randomly-generated matrices that are ill-conditioned
It will be a right way. Since the implementations of LAPACK are trusted, we do not have to check the singular cases.
Sometimes this test fails https://app.wercker.com/termoshtt/ndarray-linalg/runs/test-openblas/59f195128e888f0001cf3ef5?step=59f19584791109000146172c