QZ-decomposition #1067
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,74 +1,72 @@ | ||
| use na::dimension::Const; | ||
| use na::{DMatrix, OMatrix}; | ||
| use nl::GeneralizedEigen; | ||
| use num_complex::Complex; | ||
| use simba::scalar::ComplexField; | ||
|
|
||
| use crate::proptest::*; | ||
| use proptest::{prop_assert, prop_compose, proptest}; | ||
|
|
||
| prop_compose! { | ||
| fn f64_dynamic_dim_squares() | ||
| (n in PROPTEST_MATRIX_DIM) | ||
| (a in matrix(PROPTEST_F64,n,n), b in matrix(PROPTEST_F64,n,n)) -> (DMatrix<f64>, DMatrix<f64>){ | ||
| (a,b) | ||
| }} | ||
|
|
||
| proptest! { | ||
| #[test] | ||
| fn ge((a,b) in f64_dynamic_dim_squares()){ | ||
|
|
||
| let a_c = a.clone().map(|x| Complex::new(x, 0.0)); | ||
| let b_c = b.clone().map(|x| Complex::new(x, 0.0)); | ||
| let n = a.shape_generic().0; | ||
|
|
||
| let ge = GeneralizedEigen::new(a.clone(), b.clone()); | ||
| let (vsl,vsr) = ge.clone().eigenvectors(); | ||
|
|
||
|
|
||
| for (i,(alpha,beta)) in ge.raw_eigenvalues().iter().enumerate() { | ||
| let l_a = a_c.clone() * Complex::new(*beta, 0.0); | ||
| let l_b = b_c.clone() * *alpha; | ||
|
|
||
| prop_assert!( | ||
| relative_eq!( | ||
| ((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()), | ||
| OMatrix::zeros_generic(n, Const::<1>), | ||
| epsilon = 1.0e-5)); | ||
|
There was a problem hiding this comment. Relative comparisons tend to break down near zero - can't this make this test relatively brittle? There was a problem hiding this comment. I see what you mean. I believe there is the absdiff catch as shown here https://docs.rs/approx/0.4.0/src/approx/relative_eq.rs.html#68 to tackle the concern. Hope I'm not preaching to the choir. If I am, apologies There was a problem hiding this comment. I didn't know that It is a little unfortunate then that With all that said, it might be fine in this case, if you know that the "scale" (for example, the largest entry) of your numbers are expected to be somewhat close to There was a problem hiding this comment.
The eigenvectors seem to be guaranteed by LAPACK in their documentation
So I think scale wise we should be good |
||
|
|
||
| prop_assert!( | ||
| relative_eq!( | ||
| (vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()), | ||
| OMatrix::zeros_generic(Const::<1>, n), | ||
| epsilon = 1.0e-5)) | ||
| }; | ||
| } | ||
|
|
||
| #[test] | ||
| fn ge_static(a in matrix4(), b in matrix4()) { | ||
|
|
||
| let ge = GeneralizedEigen::new(a.clone(), b.clone()); | ||
| let a_c =a.clone().map(|x| Complex::new(x, 0.0)); | ||
| let b_c = b.clone().map(|x| Complex::new(x, 0.0)); | ||
| let (vsl,vsr) = ge.eigenvectors(); | ||
| let eigenvalues = ge.raw_eigenvalues(); | ||
|
|
||
| for (i,(alpha,beta)) in eigenvalues.iter().enumerate() { | ||
| let l_a = a_c.clone() * Complex::new(*beta, 0.0); | ||
| let l_b = b_c.clone() * *alpha; | ||
|
|
||
| prop_assert!( | ||
| relative_eq!( | ||
| ((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()), | ||
| OMatrix::zeros_generic(Const::<4>, Const::<1>), | ||
| epsilon = 1.0e-5)); | ||
| prop_assert!( | ||
| relative_eq!((vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()), | ||
| OMatrix::zeros_generic(Const::<1>, Const::<4>), | ||
| epsilon = 1.0e-5)) | ||
| } | ||
| } | ||
|
|
||
| } | ||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,74 +1,34 @@ | ||
| use na::DMatrix; | ||
| use nl::QZ; | ||
|
|
||
| use crate::proptest::*; | ||
| use proptest::{prop_assert, prop_compose, proptest}; | ||
|
|
||
| prop_compose! { | ||
| fn f64_dynamic_dim_squares() | ||
| (n in PROPTEST_MATRIX_DIM) | ||
| (a in matrix(PROPTEST_F64,n,n), b in matrix(PROPTEST_F64,n,n)) -> (DMatrix<f64>, DMatrix<f64>){ | ||
| (a,b) | ||
| }} | ||
|
|
||
| proptest! { | ||
| #[test] | ||
| fn qz((a,b) in f64_dynamic_dim_squares()) { | ||
|
|
||
| let qz = QZ::new(a.clone(), b.clone()); | ||
| let (vsl,s,t,vsr) = qz.clone().unpack(); | ||
|
|
||
| prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7)); | ||
| prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7)); | ||
|
|
||
| } | ||
|
|
||
| #[test] | ||
| fn qz_static(a in matrix4(), b in matrix4()) { | ||
| let qz = QZ::new(a.clone(), b.clone()); | ||
| let (vsl,s,t,vsr) = qz.unpack(); | ||
|
|
||
| prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7)); | ||
| prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7)); | ||
| } | ||
| } |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I missed this the first time around, but at first glance this logic does not sound to me. What does it do? I think this needs a solid comment and/or documentation.
For example, imagine scaling a matrix by
1e20, then all its eigenvalues get scaled by1e20and this logic will behave entirely differently. It seems to me at the very least that we'd need to look relative to the maximum absolute eigenvalue, but even then it's not bullet proof.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The scaling point is something I've never considered, it's a pretty solid point.
But in this scenario, LAPACK uses the eigenvalues to signal how to assemble eigenvectors. If the eigenvalue is complex and if the subsequent eigenvalue is the former's complex conjugate, the eigenvector to be assembled has it's real and imaginary parts in two subsequent columns.
Does the scaling edge case apply here seeing that it's more of an assembling problem by signalling than a numeric one? In the documentation, the way to check for realness is by looking at a_j values, this gets mapped to the imaginary value of an eigenvalue. While it's interesting for me to look deeper in general and perhaps into LAPACK code, was wondering if putting the burden of getting right for the edge case is more on LAPACK to get it right?
In either case, putting in a comment here does make senseIs this not apt to explain what is being done? It is possible I've been in the trenches long enough to not know if it doesn't map well https://github.com/dimforge/nalgebra/pull/1067/files#diff-a056a4c6b4fb629735a73fcb65e0d5223386cb47dedb6ab80d5a1b905ebd9d25R198There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think for me it's not clear why we need to do any kind of "post processing" at all. I think this would be nice to have a comment explaining this. And moreover, I would be very surprised if LAPACK expects you to use approximate comparison if it's signaling something through values - are you sure exact comparison is not more appropriate? I just don't have enough context to know what is the appropriate thing to do here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is fair. I think there is some level of post processing expected but not to the extent I have done which is more a result of not trusting LAPACK which may be unwarranted. The changes reflect what I think LAPACK expects us to do minimally