From 714f2ac98795935acecfdf47da7339a171bd2ec7 Mon Sep 17 00:00:00 2001 From: metric-space Date: Thu, 3 Feb 2022 06:36:10 -0500 Subject: [PATCH] New code and modified tests for generalized_eigenvalues --- .../src/generalized_eigenvalues.rs | 89 ++++++++++++++++--- .../tests/linalg/generalized_eigenvalues.rs | 63 ++++++++----- 2 files changed, 118 insertions(+), 34 deletions(-) diff --git a/nalgebra-lapack/src/generalized_eigenvalues.rs b/nalgebra-lapack/src/generalized_eigenvalues.rs index 6332f2db..5c273e9b 100644 --- a/nalgebra-lapack/src/generalized_eigenvalues.rs +++ b/nalgebra-lapack/src/generalized_eigenvalues.rs @@ -4,7 +4,7 @@ use serde::{Deserialize, Serialize}; use num::Zero; use num_complex::Complex; -use simba::scalar:: RealField; +use simba::scalar::RealField; use crate::ComplexHelper; use na::allocator::Allocator; @@ -14,6 +14,19 @@ use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar}; use lapack; /// Generalized eigenvalues and generalized eigenvectors(left and right) of a pair of N*N square matrices. +/// +/// Each generalized eigenvalue (lambda) satisfies determinant(A - lambda*B) = 0 +/// +/// The right eigenvector v(j) corresponding to the eigenvalue lambda(j) +/// of (A,B) satisfies +/// +/// A * v(j) = lambda(j) * B * v(j). +/// +/// The left eigenvector u(j) corresponding to the eigenvalue lambda(j) +/// of (A,B) satisfies +/// +/// u(j)**H * A = lambda(j) * u(j)**H * B . +/// where u(j)**H is the conjugate-transpose of u(j). #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize", @@ -55,11 +68,21 @@ impl GE where DefaultAllocator: Allocator + Allocator, { - /// Attempts to compute the generalized eigenvalues (and eigenvectors) via the raw returns from LAPACK's - /// dggev and sggev routines + /// Attempts to compute the generalized eigenvalues, and left and right associated eigenvectors + /// via the raw returns from LAPACK's dggev and sggev routines /// - /// For each e in generalized eigenvalues and the associated eigenvectors e_l and e_r (left andf right) - /// it satisfies e_l*a = e*e_l*b and a*e_r = e*b*e_r + /// Each generalized eigenvalue (lambda) satisfies determinant(A - lambda*B) = 0 + /// + /// The right eigenvector v(j) corresponding to the eigenvalue lambda(j) + /// of (A,B) satisfies + /// + /// A * v(j) = lambda(j) * B * v(j). + /// + /// The left eigenvector u(j) corresponding to the eigenvalue lambda(j) + /// of (A,B) satisfies + /// + /// u(j)**H * A = lambda(j) * u(j)**H * B . + /// where u(j)**H is the conjugate-transpose of u(j). /// /// Panics if the method did not converge. pub fn new(a: OMatrix, b: OMatrix) -> Self { @@ -69,8 +92,18 @@ where /// Attempts to compute the generalized eigenvalues (and eigenvectors) via the raw returns from LAPACK's /// dggev and sggev routines /// - /// For each e in generalized eigenvalues and the associated eigenvectors e_l and e_r (left andf right) - /// it satisfies e_l*a = e*e_l*b and a*e_r = e*b*e_r + /// Each generalized eigenvalue (lambda) satisfies determinant(A - lambda*B) = 0 + /// + /// The right eigenvector v(j) corresponding to the eigenvalue lambda(j) + /// of (A,B) satisfies + /// + /// A * v(j) = lambda(j) * B * v(j). + /// + /// The left eigenvector u(j) corresponding to the eigenvalue lambda(j) + /// of (A,B) satisfies + /// + /// u(j)**H * A = lambda(j) * u(j)**H * B . + /// where u(j)**H is the conjugate-transpose of u(j). /// /// Returns `None` if the method did not converge. pub fn try_new(mut a: OMatrix, mut b: OMatrix) -> Option { @@ -147,9 +180,24 @@ where } /// Calculates the generalized eigenvectors (left and right) associated with the generalized eigenvalues + /// Outputs two matrices, the first one containing the left eigenvectors of the generalized eigenvalues + /// as columns and the second matrix contains the right eigenvectors of the generalized eigenvalues + /// as columns + /// + /// The right eigenvector v(j) corresponding to the eigenvalue lambda(j) + /// of (A,B) satisfies + /// + /// A * v(j) = lambda(j) * B * v(j). + /// + /// The left eigenvector u(j) corresponding to the eigenvalue lambda(j) + /// of (A,B) satisfies + /// + /// u(j)**H * A = lambda(j) * u(j)**H * B . + /// where u(j)**H is the conjugate-transpose of u(j). pub fn eigenvectors(self) -> (OMatrix, D, D>, OMatrix, D, D>) where - DefaultAllocator: Allocator, D, D> + Allocator, D>, + DefaultAllocator: + Allocator, D, D> + Allocator, D> + Allocator<(Complex, T), D>, { let n = self.vsl.shape().0; let mut l = self @@ -199,9 +247,10 @@ where (l, r) } - /// computes the generalized eigenvalues + /// computes the generalized eigenvalues i.e values of lambda that satisfy the following equation + /// determinant(A - lambda* B) = 0 #[must_use] - pub fn eigenvalues(&self) -> OVector, D> + fn eigenvalues(&self) -> OVector, D> where DefaultAllocator: Allocator, D>, { @@ -233,6 +282,26 @@ where out } + + /// outputs the unprocessed (almost) version of generalized eigenvalues ((alphar, alpai), beta) + /// straight from LAPACK + #[must_use] + pub fn raw_eigenvalues(&self) -> OVector<(Complex, T), D> + where + DefaultAllocator: Allocator<(Complex, T), D>, + { + let mut out = Matrix::from_element_generic( + self.vsl.shape_generic().0, + Const::<1>, + (Complex::zero(), T::RealField::zero()), + ); + + for i in 0..out.len() { + out[i] = (Complex::new(self.alphar[i], self.alphai[i]), self.beta[i]) + } + + out + } } /* diff --git a/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs b/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs index 275691c8..8da21b30 100644 --- a/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs +++ b/nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs @@ -17,21 +17,29 @@ proptest! { let a_condition_no = a.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&a))); let b_condition_no = b.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&b))); - if a_condition_no.unwrap_or(200000.0) < 10.0 && b_condition_no.unwrap_or(200000.0) < 10.0 { - let a_c =a.clone().map(|x| Complex::new(x, 0.0)); - let b_c = b.clone().map(|x| Complex::new(x, 0.0)); + if a_condition_no.unwrap_or(200000.0) < 5.0 && b_condition_no.unwrap_or(200000.0) < 5.0 { + 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 ge = GE::new(a.clone(), b.clone()); - let (vsl,vsr) = ge.clone().eigenvectors(); - let eigenvalues = ge.clone().eigenvalues(); + let ge = GE::new(a.clone(), b.clone()); + let (vsl,vsr) = ge.clone().eigenvectors(); - for i in 0..n { - let left_eigenvector = &vsl.column(i); - prop_assert!(relative_eq!((left_eigenvector.transpose()*&a_c - left_eigenvector.transpose()*&b_c*eigenvalues[i]).map(|x| x.modulus()), OMatrix::zeros_generic(Const::<1>,Dynamic::new(n)) ,epsilon = 1.0e-7)); + 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; - let right_eigenvector = &vsr.column(i); - prop_assert!(relative_eq!((&a_c*right_eigenvector - &b_c*right_eigenvector*eigenvalues[i]).map(|x| x.modulus()), OMatrix::zeros_generic(Dynamic::new(n), Const::<1>) ,epsilon = 1.0e-7)); - }; + prop_assert!( + relative_eq!( + ((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()), + OMatrix::zeros_generic(Dynamic::new(n), Const::<1>), + epsilon = 1.0e-7)); + + prop_assert!( + relative_eq!( + (vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()), + OMatrix::zeros_generic(Const::<1>, Dynamic::new(n)), + epsilon = 1.0e-7)) + }; }; } @@ -40,20 +48,27 @@ proptest! { let a_condition_no = a.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&a))); let b_condition_no = b.clone().try_inverse().and_then(|x| Some(EuclideanNorm.norm(&x)* EuclideanNorm.norm(&b))); - if a_condition_no.unwrap_or(200000.0) < 10.0 && b_condition_no.unwrap_or(200000.0) < 10.0{ - let ge = GE::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.eigenvalues(); + if a_condition_no.unwrap_or(200000.0) < 5.0 && b_condition_no.unwrap_or(200000.0) < 5.0 { + let ge = GE::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 in 0..4 { - let left_eigenvector = &vsl.column(i); - prop_assert!(relative_eq!((left_eigenvector.transpose()*&a_c - left_eigenvector.transpose()*&b_c*eigenvalues[i]).map(|x| x.modulus()), OMatrix::zeros_generic(Const::<1>,Const::<4>) ,epsilon = 1.0e-7)); + 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; - let right_eigenvector = &vsr.column(i); - prop_assert!(relative_eq!((&a_c*right_eigenvector - &b_c*right_eigenvector*eigenvalues[i]).map(|x| x.modulus()), OMatrix::zeros_generic(Const::<4>, Const::<1>) ,epsilon = 1.0e-7)); - }; + 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-7)); + 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-7)) + } }; } }