ndarray_linalg/krylov/
mod.rs

1//! Krylov subspace methods
2
3use crate::types::*;
4use ndarray::*;
5
6pub mod arnoldi;
7pub mod householder;
8pub mod mgs;
9
10pub use arnoldi::{arnoldi_householder, arnoldi_mgs, Arnoldi};
11pub use householder::{householder, Householder};
12pub use mgs::{mgs, MGS};
13
14/// Q-matrix
15///
16/// - Maybe **NOT** square
17/// - Unitary for existing columns
18///
19pub type Q<A> = Array2<A>;
20
21/// R-matrix
22///
23/// - Maybe **NOT** square
24/// - Upper triangle
25///
26pub type R<A> = Array2<A>;
27
28/// H-matrix
29///
30/// - Maybe **NOT** square
31/// - Hessenberg matrix
32///
33pub type H<A> = Array2<A>;
34
35/// Array type for coefficients to the current basis
36///
37/// - The length must be `self.len() + 1`
38/// - Last component is the residual norm
39///
40pub type Coefficients<A> = Array1<A>;
41
42/// Trait for creating orthogonal basis from iterator of arrays
43///
44/// Panic
45/// -------
46/// - if the size of the input array mismatches to the dimension
47///
48/// Example
49/// -------
50///
51/// ```rust
52/// # use ndarray::*;
53/// # use ndarray_linalg::{krylov::*, *};
54/// let mut mgs = MGS::new(3, 1e-9);
55/// let coef = mgs.append(array![0.0, 1.0, 0.0]).into_coeff();
56/// close_l2(&coef, &array![1.0], 1e-9);
57///
58/// let coef = mgs.append(array![1.0, 1.0, 0.0]).into_coeff();
59/// close_l2(&coef, &array![1.0, 1.0], 1e-9);
60///
61/// // Fail if the vector is linearly dependent
62/// assert!(mgs.append(array![1.0, 2.0, 0.0]).is_dependent());
63///
64/// // You can get coefficients of dependent vector
65/// if let AppendResult::Dependent(coef) = mgs.append(array![1.0, 2.0, 0.0]) {
66///     close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9);
67/// }
68/// ```
69pub trait Orthogonalizer {
70    type Elem: Scalar;
71
72    /// Dimension of input array
73    fn dim(&self) -> usize;
74
75    /// Number of cached basis
76    fn len(&self) -> usize;
77
78    /// check if the basis spans entire space
79    fn is_full(&self) -> bool {
80        self.len() == self.dim()
81    }
82
83    fn is_empty(&self) -> bool {
84        self.len() == 0
85    }
86
87    fn tolerance(&self) -> <Self::Elem as Scalar>::Real;
88
89    /// Decompose given vector into the span of current basis and
90    /// its tangent space
91    ///
92    /// - `a` becomes the tangent vector
93    /// - The Coefficients to the current basis is returned.
94    ///
95    fn decompose(&self, a: &mut ArrayRef<Self::Elem, Ix1>) -> Coefficients<Self::Elem>;
96
97    /// Calculate the coefficient to the current basis basis
98    ///
99    /// - This will be faster than `decompose` because the construction of the residual vector may
100    ///   requires more Calculation
101    ///
102    fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Coefficients<Self::Elem>
103    where
104        S: Data<Elem = Self::Elem>;
105
106    /// Add new vector if the residual is larger than relative tolerance
107    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
108    where
109        S: Data<Elem = Self::Elem>;
110
111    /// Add new vector if the residual is larger than relative tolerance,
112    /// and return the residual vector
113    fn div_append(&mut self, a: &mut ArrayRef<Self::Elem, Ix1>) -> AppendResult<Self::Elem>;
114
115    /// Get Q-matrix of generated basis
116    fn get_q(&self) -> Q<Self::Elem>;
117}
118
119pub enum AppendResult<A> {
120    Added(Coefficients<A>),
121    Dependent(Coefficients<A>),
122}
123
124impl<A: Scalar> AppendResult<A> {
125    pub fn into_coeff(self) -> Coefficients<A> {
126        match self {
127            AppendResult::Added(c) => c,
128            AppendResult::Dependent(c) => c,
129        }
130    }
131
132    pub fn is_dependent(&self) -> bool {
133        match self {
134            AppendResult::Added(_) => false,
135            AppendResult::Dependent(_) => true,
136        }
137    }
138
139    pub fn coeff(&self) -> &Coefficients<A> {
140        match self {
141            AppendResult::Added(c) => c,
142            AppendResult::Dependent(c) => c,
143        }
144    }
145
146    pub fn residual_norm(&self) -> A::Real {
147        let c = self.coeff();
148        c[c.len() - 1].abs()
149    }
150}
151
152/// Strategy for linearly dependent vectors appearing in iterative QR decomposition
153#[derive(Clone, Copy, Debug, Eq, PartialEq)]
154pub enum Strategy {
155    /// Terminate iteration if dependent vector comes
156    Terminate,
157
158    /// Skip dependent vector
159    Skip,
160
161    /// Orthogonalize dependent vector without adding to Q,
162    /// i.e. R must be non-square like following:
163    ///
164    /// ```text
165    /// x x x x x
166    /// 0 x x x x
167    /// 0 0 0 x x
168    /// 0 0 0 0 x
169    /// ```
170    Full,
171}
172
173/// Online QR decomposition using arbitrary orthogonalizer
174pub fn qr<A, S>(
175    iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
176    mut ortho: impl Orthogonalizer<Elem = A>,
177    strategy: Strategy,
178) -> (Q<A>, R<A>)
179where
180    A: Scalar + Lapack,
181    S: Data<Elem = A>,
182{
183    assert_eq!(ortho.len(), 0);
184
185    let mut coefs = Vec::new();
186    for a in iter {
187        match ortho.append(a.into_owned()) {
188            AppendResult::Added(coef) => coefs.push(coef),
189            AppendResult::Dependent(coef) => match strategy {
190                Strategy::Terminate => break,
191                Strategy::Skip => continue,
192                Strategy::Full => coefs.push(coef),
193            },
194        }
195    }
196    let n = ortho.len();
197    let m = coefs.len();
198    let mut r = Array2::zeros((n, m).f());
199    for j in 0..m {
200        for i in 0..n {
201            if i < coefs[j].len() {
202                r[(i, j)] = coefs[j][i];
203            }
204        }
205    }
206    (ortho.get_q(), r)
207}
ndarray_linalg/krylov/mod.rs

ndarray_linalg/krylov/
mod.rs
