nalgebra/src/linalg/schur.rs

582 lines
19 KiB
Rust
Raw Normal View History

#[cfg(feature = "serde-serialize")]
2018-10-22 13:00:10 +08:00
use serde::{Deserialize, Serialize};
use approx::AbsDiffEq;
use alga::general::{Complex, Real};
use num_complex::Complex as NumComplex;
use std::cmp;
use allocator::Allocator;
use base::dimension::{Dim, DimDiff, DimSub, Dynamic, U1, U2, U3};
use base::storage::Storage;
use base::{DefaultAllocator, MatrixN, SquareMatrix, Unit, Vector2, Vector3, VectorN};
use constraint::{DimEq, ShapeConstraint};
use geometry::Reflection;
use linalg::householder;
use linalg::Hessenberg;
use linalg::givens::GivensRotation;
/// Schur decomposition of a square matrix.
///
/// If this is a real matrix, this will be a Real Schur decomposition.
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
#[cfg_attr(
feature = "serde-serialize",
2018-10-22 13:00:10 +08:00
serde(bound(
serialize = "DefaultAllocator: Allocator<N, D, D>,
MatrixN<N, D>: Serialize"
2018-10-22 13:00:10 +08:00
))
)]
#[cfg_attr(
feature = "serde-serialize",
2018-10-22 13:00:10 +08:00
serde(bound(
deserialize = "DefaultAllocator: Allocator<N, D, D>,
MatrixN<N, D>: Deserialize<'de>"
2018-10-22 13:00:10 +08:00
))
)]
#[derive(Clone, Debug)]
pub struct Schur<N: Complex, D: Dim>
2018-10-22 13:00:10 +08:00
where DefaultAllocator: Allocator<N, D, D>
2018-02-02 19:26:35 +08:00
{
q: MatrixN<N, D>,
2018-02-02 19:26:35 +08:00
t: MatrixN<N, D>,
}
impl<N: Complex, D: Dim> Copy for Schur<N, D>
2018-02-02 19:26:35 +08:00
where
DefaultAllocator: Allocator<N, D, D>,
MatrixN<N, D>: Copy,
2018-10-22 13:00:10 +08:00
{}
impl<N: Complex, D: Dim> Schur<N, D>
2018-02-02 19:26:35 +08:00
where
D: DimSub<U1>, // For Hessenberg.
ShapeConstraint: DimEq<Dynamic, DimDiff<D, U1>>, // For Hessenberg.
DefaultAllocator: Allocator<N, D, DimDiff<D, U1>>
+ Allocator<N, DimDiff<D, U1>>
+ Allocator<N, D, D>
+ Allocator<N, D>,
{
/// Computes the Schur decomposition of a square matrix.
pub fn new(m: MatrixN<N, D>) -> Self {
Self::try_new(m, N::Real::default_epsilon(), 0).unwrap()
}
/// Attempts to compute the Schur decomposition of a square matrix.
///
/// If only eigenvalues are needed, it is more efficient to call the matrix method
/// `.eigenvalues()` instead.
///
/// # Arguments
///
2018-09-24 12:48:42 +08:00
/// * `eps` tolerance used to determine when a value converged to 0.
/// * `max_niter` maximum total number of iterations performed by the algorithm. If this
/// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
/// continues indefinitely until convergence.
pub fn try_new(m: MatrixN<N, D>, eps: N::Real, max_niter: usize) -> Option<Self> {
let mut work = unsafe { VectorN::new_uninitialized_generic(m.data.shape().0, U1) };
Self::do_decompose(m, &mut work, eps, max_niter, true).map(|(q, t)| Schur {
2018-02-02 19:26:35 +08:00
q: q.unwrap(),
t: t,
})
}
2018-02-02 19:26:35 +08:00
fn do_decompose(
mut m: MatrixN<N, D>,
work: &mut VectorN<N, D>,
eps: N::Real,
2018-02-02 19:26:35 +08:00
max_niter: usize,
compute_q: bool,
2018-10-22 13:00:10 +08:00
) -> Option<(Option<MatrixN<N, D>>, MatrixN<N, D>)>
{
2018-02-02 19:26:35 +08:00
assert!(
m.is_square(),
"Unable to compute the eigenvectors and eigenvalues of a non-square matrix."
);
let dim = m.data.shape().0;
2019-03-12 20:15:02 +08:00
// Specialization would make this easier.
if dim.value() == 0 {
let vecs = Some(MatrixN::from_element_generic(dim, dim, N::zero()));
let vals = MatrixN::from_element_generic(dim, dim, N::zero());
return Some((vecs, vals));
2018-02-02 19:26:35 +08:00
} else if dim.value() == 1 {
if compute_q {
let q = MatrixN::from_element_generic(dim, dim, N::one());
return Some((Some(q), m));
2018-02-02 19:26:35 +08:00
} else {
return Some((None, m));
}
2019-03-12 20:15:02 +08:00
} else if dim.value() == 2 {
return decompose_2x2(m, compute_q);
}
let amax_m = m.camax();
m.unscale_mut(amax_m);
let hess = Hessenberg::new_with_workspace(m, work);
let mut q;
let mut t;
if compute_q {
// FIXME: could we work without unpacking? Using only the internal representation of
// hessenberg decomposition.
let (vecs, vals) = hess.unpack();
q = Some(vecs);
t = vals;
2018-02-02 19:26:35 +08:00
} else {
q = None;
t = hess.unpack_h()
}
// Implicit double-shift QR method.
let mut niter = 0;
let (mut start, mut end) = Self::delimit_subproblem(&mut t, eps, dim.value() - 1);
while end != start {
let subdim = end - start + 1;
if subdim > 2 {
let m = end - 1;
let n = end;
let h11 = t[(start + 0, start + 0)];
let h12 = t[(start + 0, start + 1)];
let h21 = t[(start + 1, start + 0)];
let h22 = t[(start + 1, start + 1)];
let h32 = t[(start + 2, start + 1)];
let hnn = t[(n, n)];
let hmm = t[(m, m)];
let hnm = t[(n, m)];
let hmn = t[(m, n)];
let tra = hnn + hmm;
let det = hnn * hmm - hnm * hmn;
2018-02-02 19:26:35 +08:00
let mut axis = Vector3::new(
h11 * h11 + h12 * h21 - tra * h11 + det,
h21 * (h11 + h22 - tra),
h21 * h32,
);
2018-02-02 19:26:35 +08:00
for k in start..n - 1 {
let (norm, not_zero) = householder::reflection_axis_mut(&mut axis);
if not_zero {
if k > start {
t[(k + 0, k - 1)] = norm;
t[(k + 1, k - 1)] = N::zero();
t[(k + 2, k - 1)] = N::zero();
}
let refl = Reflection::new(Unit::new_unchecked(axis), N::zero());
{
let krows = cmp::min(k + 4, end + 1);
let mut work = work.rows_mut(0, krows);
2018-10-22 13:00:10 +08:00
refl.reflect(
&mut t
.generic_slice_mut((k, k), (U3, Dynamic::new(dim.value() - k))),
);
2018-02-02 19:26:35 +08:00
refl.reflect_rows(
&mut t.generic_slice_mut((0, k), (Dynamic::new(krows), U3)),
&mut work,
);
}
if let Some(ref mut q) = q {
refl.reflect_rows(&mut q.generic_slice_mut((0, k), (dim, U3)), work);
}
}
axis.x = t[(k + 1, k)];
axis.y = t[(k + 2, k)];
if k < n - 2 {
axis.z = t[(k + 3, k)];
}
}
let mut axis = Vector2::new(axis.x, axis.y);
let (norm, not_zero) = householder::reflection_axis_mut(&mut axis);
if not_zero {
let refl = Reflection::new(Unit::new_unchecked(axis), N::zero());
t[(m, m - 1)] = norm;
t[(n, m - 1)] = N::zero();
{
let mut work = work.rows_mut(0, end + 1);
2018-10-22 13:00:10 +08:00
refl.reflect(
&mut t.generic_slice_mut((m, m), (U2, Dynamic::new(dim.value() - m))),
);
2018-02-02 19:26:35 +08:00
refl.reflect_rows(
&mut t.generic_slice_mut((0, m), (Dynamic::new(end + 1), U2)),
&mut work,
);
}
if let Some(ref mut q) = q {
refl.reflect_rows(&mut q.generic_slice_mut((0, m), (dim, U2)), work);
}
}
2018-02-02 19:26:35 +08:00
} else {
// Decouple the 2x2 block if it has real eigenvalues.
if let Some(rot) = compute_2x2_basis(&t.fixed_slice::<U2, U2>(start, start)) {
let inv_rot = rot.inverse();
2018-02-02 19:26:35 +08:00
inv_rot.rotate(&mut t.generic_slice_mut(
(start, start),
(U2, Dynamic::new(dim.value() - start)),
));
2018-10-22 13:00:10 +08:00
rot.rotate_rows(
&mut t.generic_slice_mut((0, start), (Dynamic::new(end + 1), U2)),
);
t[(end, start)] = N::zero();
if let Some(ref mut q) = q {
rot.rotate_rows(&mut q.generic_slice_mut((0, start), (dim, U2)));
}
}
// Check if we reached the beginning of the matrix.
if end > 2 {
end -= 2;
2018-02-02 19:26:35 +08:00
} else {
break;
}
}
let sub = Self::delimit_subproblem(&mut t, eps, end);
start = sub.0;
2018-02-02 19:26:35 +08:00
end = sub.1;
niter += 1;
if niter == max_niter {
return None;
}
}
t.scale_mut(amax_m);
Some((q, t))
}
/// Computes the eigenvalues of the decomposed matrix.
fn do_eigenvalues(t: &MatrixN<N, D>, out: &mut VectorN<N, D>) -> bool {
let dim = t.nrows();
let mut m = 0;
while m < dim - 1 {
let n = m + 1;
if t[(n, m)].is_zero() {
out[m] = t[(m, m)];
m += 1;
2018-02-02 19:26:35 +08:00
} else {
// Complex eigenvalue.
return false;
}
}
if m == dim - 1 {
out[m] = t[(m, m)];
}
true
}
/// Computes the complex eigenvalues of the decomposed matrix.
fn do_complex_eigenvalues(t: &MatrixN<N, D>, out: &mut VectorN<NumComplex<N>, D>)
where N: Real,
DefaultAllocator: Allocator<NumComplex<N>, D> {
let dim = t.nrows();
let mut m = 0;
while m < dim - 1 {
let n = m + 1;
if t[(n, m)].is_zero() {
out[m] = NumComplex::new(t[(m, m)], N::zero());
m += 1;
2018-02-02 19:26:35 +08:00
} else {
// Solve the 2x2 eigenvalue subproblem.
let hmm = t[(m, m)];
let hnm = t[(n, m)];
let hmn = t[(m, n)];
let hnn = t[(n, n)];
2018-02-02 19:26:35 +08:00
let tra = hnn + hmm;
let det = hnn * hmm - hnm * hmn;
let discr = tra * tra * ::convert(0.25) - det;
// All 2x2 blocks have negative discriminant because we already decoupled those
// with positive eigenvalues..
let sqrt_discr = NumComplex::new(N::zero(), (-discr).sqrt());
out[m] = NumComplex::new(tra * ::convert(0.5), N::zero()) + sqrt_discr;
out[m + 1] = NumComplex::new(tra * ::convert(0.5), N::zero()) - sqrt_discr;
m += 2;
}
}
if m == dim - 1 {
out[m] = NumComplex::new(t[(m, m)], N::zero());
}
}
fn delimit_subproblem(t: &mut MatrixN<N, D>, eps: N::Real, end: usize) -> (usize, usize)
2018-02-02 19:26:35 +08:00
where
D: DimSub<U1>,
DefaultAllocator: Allocator<N, DimDiff<D, U1>>,
{
let mut n = end;
while n > 0 {
let m = n - 1;
if t[(n, m)].modulus() <= eps * (t[(n, n)].modulus() + t[(m, m)].modulus()) {
t[(n, m)] = N::zero();
2018-02-02 19:26:35 +08:00
} else {
break;
}
n -= 1;
}
if n == 0 {
return (0, 0);
}
let mut new_start = n - 1;
while new_start > 0 {
let m = new_start - 1;
let off_diag = t[(new_start, m)];
2018-02-02 19:26:35 +08:00
if off_diag.is_zero()
|| off_diag.modulus() <= eps * (t[(new_start, new_start)].modulus() + t[(m, m)].modulus())
2018-02-02 19:26:35 +08:00
{
t[(new_start, m)] = N::zero();
break;
}
new_start -= 1;
}
(new_start, n)
}
/// Retrieves the unitary matrix `Q` and the upper-quasitriangular matrix `T` such that the
/// decomposed matrix equals `Q * T * Q.transpose()`.
pub fn unpack(self) -> (MatrixN<N, D>, MatrixN<N, D>) {
(self.q, self.t)
}
/// Computes the real eigenvalues of the decomposed matrix.
///
/// Return `None` if some eigenvalues are complex.
pub fn eigenvalues(&self) -> Option<VectorN<N, D>> {
let mut out = unsafe { VectorN::new_uninitialized_generic(self.t.data.shape().0, U1) };
if Self::do_eigenvalues(&self.t, &mut out) {
Some(out)
2018-02-02 19:26:35 +08:00
} else {
None
}
}
/// Computes the complex eigenvalues of the decomposed matrix.
pub fn complex_eigenvalues(&self) -> VectorN<NumComplex<N>, D>
where N: Real,
DefaultAllocator: Allocator<NumComplex<N>, D> {
let mut out = unsafe { VectorN::new_uninitialized_generic(self.t.data.shape().0, U1) };
Self::do_complex_eigenvalues(&self.t, &mut out);
out
}
}
fn decompose_2x2<N: Complex, D: Dim>(
2018-02-02 19:26:35 +08:00
mut m: MatrixN<N, D>,
compute_q: bool,
) -> Option<(Option<MatrixN<N, D>>, MatrixN<N, D>)>
where
DefaultAllocator: Allocator<N, D, D>,
{
let dim = m.data.shape().0;
let mut q = None;
match compute_2x2_basis(&m.fixed_slice::<U2, U2>(0, 0)) {
Some(rot) => {
let mut m = m.fixed_slice_mut::<U2, U2>(0, 0);
let inv_rot = rot.inverse();
inv_rot.rotate(&mut m);
rot.rotate_rows(&mut m);
if compute_q {
// XXX: we have to build the matrix manually because
// rot.to_rotation_matrix().unwrap() causes an ICE.
2019-03-18 18:23:19 +08:00
let c = N::from_real(rot.c());
2018-02-02 19:26:35 +08:00
q = Some(MatrixN::from_column_slice_generic(
dim,
dim,
2019-03-18 18:23:19 +08:00
&[c, rot.s(), -rot.s().conjugate(), c],
2018-02-02 19:26:35 +08:00
));
}
2018-02-02 19:26:35 +08:00
}
2018-10-22 13:00:10 +08:00
None => {
if compute_q {
q = Some(MatrixN::identity_generic(dim, dim));
}
}
};
Some((q, m))
}
fn compute_2x2_eigvals<N: Complex, S: Storage<N, U2, U2>>(
2018-02-02 19:26:35 +08:00
m: &SquareMatrix<N, U2, S>,
) -> Option<(N, N)> {
// Solve the 2x2 eigenvalue subproblem.
let h00 = m[(0, 0)];
let h10 = m[(1, 0)];
let h01 = m[(0, 1)];
let h11 = m[(1, 1)];
2019-03-12 20:15:02 +08:00
// NOTE: this discriminant computation is more stable than the
// one based on the trace and determinant: 0.25 * tra * tra - det
2019-03-12 20:15:02 +08:00
// because it ensures positiveness for symmetric matrices.
2018-02-02 19:26:35 +08:00
let val = (h00 - h11) * ::convert(0.5);
let discr = h10 * h01 + val * val;
discr.try_sqrt().map(|sqrt_discr| {
2018-02-02 19:26:35 +08:00
let half_tra = (h00 + h11) * ::convert(0.5);
(half_tra + sqrt_discr, half_tra - sqrt_discr)
})
}
// Computes the 2x2 transformation that upper-triangulates a 2x2 matrix with real eigenvalues.
/// Computes the singular vectors for a 2x2 matrix.
///
/// Returns `None` if the matrix has complex eigenvalues, or is upper-triangular. In both case,
/// the basis is the identity.
fn compute_2x2_basis<N: Complex, S: Storage<N, U2, U2>>(
2018-02-02 19:26:35 +08:00
m: &SquareMatrix<N, U2, S>,
) -> Option<GivensRotation<N>> {
let h10 = m[(1, 0)];
if h10.is_zero() {
return None;
}
if let Some((eigval1, eigval2)) = compute_2x2_eigvals(m) {
2019-03-12 20:15:02 +08:00
let x1 = eigval1 - m[(1, 1)];
let x2 = eigval2 - m[(1, 1)];
// NOTE: Choose the one that yields a larger x component.
// This is necessary for numerical stability of the normalization of the complex
// number.
if x1.modulus() > x2.modulus() {
2019-03-18 18:23:19 +08:00
Some(GivensRotation::new(x1, h10).0)
2018-02-02 19:26:35 +08:00
} else {
2019-03-18 18:23:19 +08:00
Some(GivensRotation::new(x2, h10).0)
}
2018-02-02 19:26:35 +08:00
} else {
None
}
}
impl<N: Complex, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S>
2018-02-02 19:26:35 +08:00
where
D: DimSub<U1>, // For Hessenberg.
ShapeConstraint: DimEq<Dynamic, DimDiff<D, U1>>, // For Hessenberg.
DefaultAllocator: Allocator<N, D, DimDiff<D, U1>>
+ Allocator<N, DimDiff<D, U1>>
+ Allocator<N, D, D>
+ Allocator<N, D>,
{
/// Computes the Schur decomposition of a square matrix.
pub fn schur(self) -> Schur<N, D> {
Schur::new(self.into_owned())
}
/// Attempts to compute the Schur decomposition of a square matrix.
///
/// If only eigenvalues are needed, it is more efficient to call the matrix method
/// `.eigenvalues()` instead.
///
/// # Arguments
///
2018-09-24 12:48:42 +08:00
/// * `eps` tolerance used to determine when a value converged to 0.
/// * `max_niter` maximum total number of iterations performed by the algorithm. If this
/// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
/// continues indefinitely until convergence.
pub fn try_schur(self, eps: N::Real, max_niter: usize) -> Option<Schur<N, D>> {
Schur::try_new(self.into_owned(), eps, max_niter)
}
/// Computes the eigenvalues of this matrix.
pub fn eigenvalues(&self) -> Option<VectorN<N, D>> {
2018-02-02 19:26:35 +08:00
assert!(
self.is_square(),
"Unable to compute eigenvalues of a non-square matrix."
);
2018-02-02 19:26:35 +08:00
let mut work = unsafe { VectorN::new_uninitialized_generic(self.data.shape().0, U1) };
2018-09-24 12:48:42 +08:00
// Special case for 2x2 matrices.
if self.nrows() == 2 {
// FIXME: can we avoid this slicing
// (which is needed here just to transform D to U2)?
let me = self.fixed_slice::<U2, U2>(0, 0);
return match compute_2x2_eigvals(&me) {
Some((a, b)) => {
work[0] = a;
work[1] = b;
Some(work)
2018-02-02 19:26:35 +08:00
}
None => None,
};
}
// FIXME: add balancing?
let schur = Schur::do_decompose(
2018-02-02 19:26:35 +08:00
self.clone_owned(),
&mut work,
N::Real::default_epsilon(),
2018-02-02 19:26:35 +08:00
0,
false,
2018-10-22 13:00:10 +08:00
)
.unwrap();
if Schur::do_eigenvalues(&schur.1, &mut work) {
Some(work)
2018-02-02 19:26:35 +08:00
} else {
None
}
}
/// Computes the eigenvalues of this matrix.
pub fn complex_eigenvalues(&self) -> VectorN<NumComplex<N>, D>
2018-02-02 19:26:35 +08:00
// FIXME: add balancing?
where N: Real,
DefaultAllocator: Allocator<NumComplex<N>, D> {
let dim = self.data.shape().0;
let mut work = unsafe { VectorN::new_uninitialized_generic(dim, U1) };
let schur = Schur::do_decompose(
2018-02-02 19:26:35 +08:00
self.clone_owned(),
&mut work,
N::default_epsilon(),
0,
false,
2018-10-22 13:00:10 +08:00
)
.unwrap();
let mut eig = unsafe { VectorN::new_uninitialized_generic(dim, U1) };
Schur::do_complex_eigenvalues(&schur.1, &mut eig);
eig
}
}