New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
SVD produces wrong singular values #1072
Comments
Thanks for reporting! Ugh, that's pretty bad. Here's a modified version which also compiles for older versions of nalgebra: use nalgebra::{DMatrix, DVector};
fn main() {
let x = DMatrix::from_row_slice(2, 2, &[-6.206610118536945f64, -3.67612186839874, -1.2755730783423473, 6.047238193479124]);
// let mut x_svd = x.try_svd(true, true, 1e-8, 0).unwrap();// This gives almost the same result as the line below.
let mut x_svd = x.svd(true, true); // This gives almost the same result as the line above.
println!("Original singular values: {}", x_svd.singular_values);
x_svd.singular_values = DVector::from_row_slice(&[1.0, 0.0]);
println!("Modified singular values: {}", x_svd.singular_values);
let y = x_svd.recompose().unwrap();
let y_svd = y.clone().svd(true, true); // This produces wrong singular values!
println!("Recomposed singular values (default eps): {}", y_svd.singular_values);
let y_svd = y.try_svd(true, true, 1e-8, 0).unwrap(); // This produces correct singular values! For small enough `eps` the wrong behavior occurs.
println!("Recomposed singular values (eps = 1e-8): {}", y_svd.singular_values);
} It seems to always have been broken, although the incorrect exact values produced seems to have changed a little between |
Update: it seems setting |
I briefly looked at the SVD implementation. I believe it's probably based off of the SVD algorithm lined out in the book Matrix Computations by Golub & Van Loan. I couldn't find any comments though... @sebcrozet can hopefully clarify that. Here's an excerpt: Note that it says that EDIT: Simply multiplying by a factor 2 seems to be sufficient in this case. I'd probably prefer something like 4 to have a bit more leeway, however. Is there a way we can test this problem more thoroughly? |
Increasing use nalgebra::proptest::matrix;
use nalgebra::DVector;
use proptest::prelude::*;
proptest! {
fn test_singular_values(x in matrix(-10f64..=10f64, 1..=5, 1..=5)) {
let mut x_svd = x.svd(true, true);
let mut true_singular_vals = DVector::zeros(x_svd.singular_values.len());
true_singular_vals[0] = 1.0;
x_svd.singular_values = true_singular_vals.clone();
let y = x_svd.recompose().unwrap();
let y_svd = y.try_svd(true, true, 4.0 * f64::EPSILON, 0).unwrap();
prop_assert!(y_svd.singular_values.relative_eq(&true_singular_vals, 1e-12, 1e-12));
}
}
fn main() {
test_singular_values();
} produces correct results. Because you mentioned the SVD algorithm from the book by Golub & Van Loan, I took a look at the GNU scientific library (written in C) that uses this algorithm to compute the SVD. It seems that in their implementation the |
|
@cfunky: I had forgotten about or entirely missed #983. I would suspect there's more to the story here than just the Given we rely on SVD for our work, this is personally very important to me to resolve, but unfortunately I also don't have the bandwidth to dig into this any time soon. Any assistance here would be much appreciated! |
Sorry for the delayed response. I'd be interested in investigating this further, as I also rely on nalgebra for work. I probably won't have a lot of time in the next couple of weeks due to upcoming deadlines, but would be happy to take a closer look afterwards. |
I'm somewhat cautiously closing this as I hope/believe it has been addressed by #1089. Feel free to re-open or make a new issue if similar issues crop up again, however! |
The SVD with the default tolerance parameter
eps
can produce significantly wrong singular values. This behavior occurs in nalgebra version 0.30.1, potentially also in other versions (I did not try this). It looks like the problem occurs when, as in the example above, an singular value is close to zero and the tolerance parametereps
ofMatrix::try_svd(...)
orMatrix::try_svd_unordered(...)
is set to a very low value, as is done internally, e.g., inMatrix::svd(...)
. Here's a minimal example reproducing this issue:Link to playground
The above code prints
The text was updated successfully, but these errors were encountered: