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

1use crate::math::Quantity;
2use crate::units::Dimensionless;
3
4use super::{
5    super::{
6        Rank2, Tensor, TensorArray, TensorError,
7        rank_0::{TensorRank0, list::TensorRank0List},
8        rank_4::TensorRank4,
9    },
10    TensorRank2,
11    eigen::reconstruct_symmetric,
12};
13use crate::math::assert::Assert;
14
15impl<I> TensorRank2<3, I, I, Dimensionless> {
16    /// Returns the matrix power of the 3x3 symmetric tensor.
17    pub fn powm(&self, exponent: TensorRank0) -> Result<Self, TensorError> {
18        if self.is_diagonal() {
19            let mut powm = TensorRank2::zero();
20            powm.iter_mut()
21                .enumerate()
22                .zip(self.iter())
23                .for_each(|((i, powm_i), self_i)| powm_i[i] = self_i[i].powf(exponent));
24            Ok(powm)
25        } else {
26            let tensor = self - &TensorRank2::identity();
27            let norm = tensor.norm();
28            if norm < 1e-2 {
29                let num_terms = if norm < 1e-4 {
30                    2
31                } else if norm < 1e-3 {
32                    3
33                } else {
34                    5
35                };
36                let mut powm = TensorRank2::identity();
37                let mut term = tensor.clone();
38                let mut coefficient = exponent;
39                (1..=num_terms).for_each(|k| {
40                    powm += &term * coefficient;
41                    term *= &tensor;
42                    coefficient *= (exponent - k as TensorRank0) / (k as TensorRank0 + 1.0);
43                });
44                Ok(powm)
45            } else if self.is_symmetric() {
46                let (eigenvalues, eigenvectors) = self.eigen()?;
47                Self::powm_from_eigen(&eigenvalues, &eigenvectors, exponent)
48            } else {
49                panic!("Matrix power only implemented for symmetric cases")
50            }
51        }
52    }
53    /// Returns the matrix power from an eigendecomposition obtained from [`Self::eigen`].
54    pub fn powm_from_eigen(
55        eigenvalues: &TensorRank0List<3>,
56        eigenvectors: &Self,
57        exponent: TensorRank0,
58    ) -> Result<Self, TensorError> {
59        let powered: TensorRank0List<3> = eigenvalues
60            .iter()
61            .map(|eigenvalue| eigenvalue.powf(exponent))
62            .collect();
63        if powered.iter().any(|value| !value.is_finite()) {
64            return Err(TensorError::NotPositiveDefinite);
65        }
66        Ok(reconstruct_symmetric(powered, eigenvectors.clone()))
67    }
68    /// Returns the derivative of the matrix power of the 3x3 symmetric tensor.
69    pub fn dpowm(
70        &self,
71        exponent: TensorRank0,
72    ) -> Result<TensorRank4<3, I, I, I, I, Dimensionless>, TensorError> {
73        if self.is_diagonal() {
74            let mut dpowm = TensorRank4::zero();
75            dpowm.iter_mut().enumerate().for_each(|(i, dpowm_i)| {
76                dpowm_i.iter_mut().enumerate().for_each(|(j, dpowm_ij)| {
77                    dpowm_ij.iter_mut().enumerate().for_each(|(k, dpowm_ijk)| {
78                        dpowm_ijk
79                            .iter_mut()
80                            .enumerate()
81                            .filter(|(l, _)| i == k && &j == l)
82                            .for_each(|(_, dpowm_ijkl)| {
83                                *dpowm_ijkl = if Assert::default()
84                                    .eq_within_tols(self[i][i], &self[j][j])
85                                    .is_ok()
86                                {
87                                    exponent * self[j][j].powf(exponent - 1.0)
88                                } else {
89                                    (self[i][i].powf(exponent) - self[j][j].powf(exponent))
90                                        / (self[i][i] - self[j][j])
91                                }
92                            })
93                    })
94                })
95            });
96            Ok(dpowm)
97        } else if self.is_symmetric() {
98            let (eigenvalues, eigenvectors) = self.eigen()?;
99            Self::dpowm_from_eigen(&eigenvalues, &eigenvectors, exponent)
100        } else {
101            panic!("Matrix power only implemented for symmetric cases")
102        }
103    }
104    /// Returns the derivative of the matrix power from an eigendecomposition obtained from
105    /// [`Self::eigen`].
106    pub fn dpowm_from_eigen(
107        eigenvalues: &TensorRank0List<3>,
108        eigenvectors: &Self,
109        exponent: TensorRank0,
110    ) -> Result<TensorRank4<3, I, I, I, I, Dimensionless>, TensorError> {
111        let divided_difference: Self = eigenvalues
112            .iter()
113            .map(|eigenvalue_i| {
114                eigenvalues
115                    .iter()
116                    .map(|eigenvalue_j| {
117                        if Assert::default()
118                            .eq_within_tols(eigenvalue_i, eigenvalue_j)
119                            .is_ok()
120                        {
121                            exponent * eigenvalue_j.powf(exponent - 1.0)
122                        } else {
123                            (eigenvalue_i.powf(exponent) - eigenvalue_j.powf(exponent))
124                                / (eigenvalue_i - eigenvalue_j)
125                        }
126                    })
127                    .collect()
128            })
129            .collect();
130        if divided_difference
131            .iter()
132            .flat_map(|row| row.iter())
133            .any(|value| !value.value().is_finite())
134        {
135            return Err(TensorError::NotPositiveDefinite);
136        }
137        let eigenvectors_transposed = eigenvectors.transpose();
138        Ok(eigenvectors_transposed.iter().map(|eigenvector_i|
139            eigenvectors_transposed.iter().map(|eigenvector_j|
140                eigenvectors_transposed.iter().map(|eigenvector_k|
141                    eigenvectors_transposed.iter().map(|eigenvector_l|
142                        eigenvector_i.iter().zip(eigenvector_k.iter().zip(divided_difference.iter())).map(|(eigenvector_ip, (eigenvector_kp, divided_difference_p))|
143                            eigenvector_j.iter().zip(eigenvector_l.iter().zip(divided_difference_p.iter())).map(|(eigenvector_jq, (eigenvector_lq, divided_difference_pq))|
144                                eigenvector_ip * eigenvector_kp * divided_difference_pq * eigenvector_jq * eigenvector_lq
145                            ).sum::<Quantity>()
146                        ).sum::<Quantity>()
147                    ).collect()
148                ).collect()
149            ).collect()
150        ).collect())
151    }
152}
153
154/// A cached eigendecomposition of a symmetric tensor, letting [`Self::powm`]/[`Self::dpowm`] be
155/// evaluated at several exponents while paying for only one cubic eigensolve.
156pub enum Spectrum<I> {
157    Eigen(TensorRank0List<3>, TensorRank2<3, I, I, Dimensionless>),
158    Fallback(TensorRank2<3, I, I, Dimensionless>),
159}
160
161impl<I> Spectrum<I> {
162    /// Caches the eigendecomposition of the tensor, when one is needed.
163    pub fn new(tensor: &TensorRank2<3, I, I, Dimensionless>) -> Result<Self, TensorError> {
164        if tensor.is_diagonal() || (tensor - &TensorRank2::identity()).norm() < 1e-2 {
165            Ok(Self::Fallback(tensor.clone()))
166        } else {
167            let (eigenvalues, eigenvectors) = tensor.eigen()?;
168            Ok(Self::Eigen(eigenvalues, eigenvectors))
169        }
170    }
171    /// Returns the matrix power at the given exponent, reusing the cached eigendecomposition.
172    pub fn powm(
173        &self,
174        exponent: TensorRank0,
175    ) -> Result<TensorRank2<3, I, I, Dimensionless>, TensorError> {
176        match self {
177            Self::Eigen(eigenvalues, eigenvectors) => {
178                TensorRank2::powm_from_eigen(eigenvalues, eigenvectors, exponent)
179            }
180            Self::Fallback(tensor) => tensor.powm(exponent),
181        }
182    }
183    /// Returns the derivative of the matrix power at the given exponent, reusing the cached
184    /// eigendecomposition.
185    pub fn dpowm(
186        &self,
187        exponent: TensorRank0,
188    ) -> Result<TensorRank4<3, I, I, I, I, Dimensionless>, TensorError> {
189        match self {
190            Self::Eigen(eigenvalues, eigenvectors) => {
191                TensorRank2::dpowm_from_eigen(eigenvalues, eigenvectors, exponent)
192            }
193            Self::Fallback(tensor) => tensor.dpowm(exponent),
194        }
195    }
196}
conspire/math/tensor/rank_2/power/mod.rs

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