conspire/math/tensor/rank_2/eigen/
mod.rs

1use crate::math::Quantity;
2use crate::units::Dimensionless;
3use std::f64::consts::TAU;
4
5use super::{
6    super::{
7        Rank2, Tensor, TensorArray, TensorError,
8        rank_0::{TensorRank0, list::TensorRank0List},
9        rank_1::{TensorRank1, cross::CrossProduct},
10    },
11    TensorRank2,
12};
13use crate::{ABS_TOL, math::assert::Assert};
14
15impl<I> TensorRank2<3, I, I, Dimensionless> {
16    /// Returns the invariants of the 3x3 symmetric tensor.
17    pub fn invariants(&self) -> TensorRank0List<3> {
18        TensorRank0List::from([
19            self.trace().value(),
20            self.second_invariant().value(),
21            self.determinant(),
22        ])
23    }
24    /// Returns the eigenvalues and (row-wise) eigenvectors of the 3x3 symmetric tensor.
25    ///
26    /// Reuse this alongside [`TensorRank2::powm_from_eigen`]/[`TensorRank2::dpowm_from_eigen`]
27    /// to evaluate several exponents against the same tensor while sharing one
28    /// eigendecomposition, instead of paying for a fresh cubic solve on every
29    /// [`TensorRank2::powm`]/[`TensorRank2::dpowm`] call.
30    pub fn eigen(&self) -> Result<(TensorRank0List<3>, Self), TensorError> {
31        let eigenvalues = solve_cubic_symmetric(self.invariants())?;
32        let eigenvectors = find_orthonormal_eigenvectors(&eigenvalues, self);
33        Ok((eigenvalues, eigenvectors))
34    }
35}
36
37/// Shared by [`super::logarithm`] and [`super::power`].
38pub(super) fn solve_cubic_symmetric(
39    coefficients: TensorRank0List<3>,
40) -> Result<TensorRank0List<3>, TensorError> {
41    let c2 = coefficients[0];
42    let c1 = coefficients[1];
43    let c0 = coefficients[2];
44    let p = c1 - c2 * c2 / 3.0;
45    let q = -(2.0 * c2.powi(3) - 9.0 * c2 * c1 + 27.0 * c0) / 27.0;
46    if p.abs() < ABS_TOL {
47        let t = (-q).cbrt();
48        let lambda = t + c2 / 3.0;
49        return Ok(TensorRank0List::from([lambda; _]));
50    }
51    let discriminant = -4.0 * p * p * p - 27.0 * q * q;
52    let scale = (4.0 * p * p * p).abs().max(27.0 * q * q);
53    if discriminant.abs() <= 1e-13 * scale {
54        let r = (q / 2.0).cbrt();
55        let lambda_double = r + c2 / 3.0;
56        let lambda_simple = -2.0 * r + c2 / 3.0;
57        let lambdas = if lambda_double >= lambda_simple {
58            [lambda_double, lambda_double, lambda_simple]
59        } else {
60            [lambda_simple, lambda_double, lambda_double]
61        };
62        Ok(TensorRank0List::from(lambdas))
63    } else if discriminant > 0.0 {
64        let sqrt_term = (-p / 3.0).sqrt();
65        let cos_arg = 3.0 * q / (2.0 * p * (-p / 3.0).sqrt());
66        let cos_arg = cos_arg.clamp(-1.0, 1.0);
67        let theta = cos_arg.acos();
68        let mut lambdas = [
69            2.0 * sqrt_term * (theta / 3.0).cos() + c2 / 3.0,
70            2.0 * sqrt_term * ((theta + TAU) / 3.0).cos() + c2 / 3.0,
71            2.0 * sqrt_term * ((theta + 2.0 * TAU) / 3.0).cos() + c2 / 3.0,
72        ];
73        lambdas.iter_mut().for_each(|lambda| {
74            for _ in 0..2 {
75                let x = *lambda;
76                let f = x * x * x - c2 * x * x + c1 * x - c0;
77                let f_prime = 3.0 * x * x - 2.0 * c2 * x + c1;
78                if f_prime.abs() < ABS_TOL {
79                    break;
80                }
81                *lambda -= f / f_prime;
82            }
83        });
84        lambdas.sort_by(|a, b| b.partial_cmp(a).unwrap());
85        Ok(TensorRank0List::from(lambdas))
86    } else {
87        Err(TensorError::SymmetricMatrixComplexEigenvalues)
88    }
89}
90
91/// Shared by [`super::logarithm`] and [`super::power`]; see [`solve_cubic_symmetric`].
92pub(super) fn find_orthonormal_eigenvectors<I>(
93    eigenvalues: &TensorRank0List<3>,
94    tensor: &TensorRank2<3, I, I, Dimensionless>,
95) -> TensorRank2<3, I, I, Dimensionless> {
96    if Assert::default()
97        .eq_within_tols(eigenvalues[0], &eigenvalues[1])
98        .is_ok()
99    {
100        let mut eigenvectors = TensorRank2::zero();
101        eigenvectors[2] = eigenvector_symmetric(eigenvalues[2], tensor);
102        eigenvectors[0] = orthogonal_unit_vector(&eigenvectors[2]);
103        eigenvectors[1] = eigenvectors[2].cross(&eigenvectors[0]);
104        eigenvectors
105    } else if Assert::default()
106        .eq_within_tols(eigenvalues[1], &eigenvalues[2])
107        .is_ok()
108    {
109        let mut eigenvectors = TensorRank2::zero();
110        eigenvectors[0] = eigenvector_symmetric(eigenvalues[0], tensor);
111        eigenvectors[1] = orthogonal_unit_vector(&eigenvectors[0]);
112        eigenvectors[2] = eigenvectors[0].cross(&eigenvectors[1]);
113        eigenvectors
114    } else {
115        let mut eigenvectors = eigenvalues
116            .iter()
117            .map(|&eigenvalue| eigenvector_symmetric(eigenvalue, tensor))
118            .collect::<TensorRank2<3, I, I, Dimensionless>>();
119        eigenvectors[0].normalize();
120        let proj1 = &eigenvectors[1] * &eigenvectors[0];
121        for i in 0..3 {
122            let projected = proj1 * eigenvectors[0][i];
123            eigenvectors[1][i] -= projected;
124        }
125        eigenvectors[1].normalize();
126        eigenvectors[2] = eigenvectors[0].cross(&eigenvectors[1]);
127        eigenvectors
128    }
129}
130
131fn orthogonal_unit_vector<I>(
132    vector: &TensorRank1<3, I, Dimensionless>,
133) -> TensorRank1<3, I, Dimensionless> {
134    let axis = vector
135        .iter()
136        .enumerate()
137        .min_by(|(_, a), (_, b)| a.abs().partial_cmp(&b.abs()).unwrap())
138        .map(|(i, _)| i)
139        .unwrap();
140    let mut other = TensorRank1::<3, I, Dimensionless>::zero();
141    other[axis] = Quantity::new(1.0);
142    vector.cross(&other).normalized()
143}
144
145fn eigenvector_symmetric<I>(
146    eigenvalue: TensorRank0,
147    tensor: &TensorRank2<3, I, I, Dimensionless>,
148) -> TensorRank1<3, I, Dimensionless> {
149    let m = tensor - TensorRank2::identity() * eigenvalue;
150    [m[1].cross(&m[2]), m[0].cross(&m[2]), m[0].cross(&m[1])]
151        .into_iter()
152        .max_by(|a, b| a.norm().partial_cmp(&b.norm()).unwrap())
153        .unwrap()
154        .normalized()
155}
156
157/// Shared by [`super::logarithm`] and [`super::power`]; see [`solve_cubic_symmetric`].
158pub(super) fn reconstruct_symmetric<I>(
159    eigenvalues: TensorRank0List<3>,
160    eigenvectors: TensorRank2<3, I, I, Dimensionless>,
161) -> TensorRank2<3, I, I, Dimensionless> {
162    let mut tensor = TensorRank2::zero();
163    eigenvalues
164        .iter()
165        .zip(eigenvectors.iter())
166        .for_each(|(eigenvalue, eigenvector)| {
167            tensor
168                .iter_mut()
169                .zip(eigenvector.iter())
170                .for_each(|(tensor_i, eigenvector_i)| {
171                    tensor_i.iter_mut().zip(eigenvector.iter()).for_each(
172                        |(tensor_ij, eigenvector_j)| {
173                            *tensor_ij += eigenvalue * eigenvector_i * eigenvector_j
174                        },
175                    )
176                })
177        });
178    tensor
179}
conspire/math/tensor/rank_2/eigen/mod.rs

conspire/math/tensor/rank_2/eigen/
mod.rs
