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<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Coefficients<Self::Elem>
96    where
97        S: DataMut<Elem = Self::Elem>;
98
99    /// Calculate the coefficient to the current basis basis
100    ///
101    /// - This will be faster than `decompose` because the construction of the residual vector may
102    ///   requires more Calculation
103    ///
104    fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Coefficients<Self::Elem>
105    where
106        S: Data<Elem = Self::Elem>;
107
108    /// Add new vector if the residual is larger than relative tolerance
109    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
110    where
111        S: Data<Elem = Self::Elem>;
112
113    /// Add new vector if the residual is larger than relative tolerance,
114    /// and return the residual vector
115    fn div_append<S>(&mut self, a: &mut ArrayBase<S, Ix1>) -> AppendResult<Self::Elem>
116    where
117        S: DataMut<Elem = Self::Elem>;
118
119    /// Get Q-matrix of generated basis
120    fn get_q(&self) -> Q<Self::Elem>;
121}
122
123pub enum AppendResult<A> {
124    Added(Coefficients<A>),
125    Dependent(Coefficients<A>),
126}
127
128impl<A: Scalar> AppendResult<A> {
129    pub fn into_coeff(self) -> Coefficients<A> {
130        match self {
131            AppendResult::Added(c) => c,
132            AppendResult::Dependent(c) => c,
133        }
134    }
135
136    pub fn is_dependent(&self) -> bool {
137        match self {
138            AppendResult::Added(_) => false,
139            AppendResult::Dependent(_) => true,
140        }
141    }
142
143    pub fn coeff(&self) -> &Coefficients<A> {
144        match self {
145            AppendResult::Added(c) => c,
146            AppendResult::Dependent(c) => c,
147        }
148    }
149
150    pub fn residual_norm(&self) -> A::Real {
151        let c = self.coeff();
152        c[c.len() - 1].abs()
153    }
154}
155
156/// Strategy for linearly dependent vectors appearing in iterative QR decomposition
157#[derive(Clone, Copy, Debug, Eq, PartialEq)]
158pub enum Strategy {
159    /// Terminate iteration if dependent vector comes
160    Terminate,
161
162    /// Skip dependent vector
163    Skip,
164
165    /// Orthogonalize dependent vector without adding to Q,
166    /// i.e. R must be non-square like following:
167    ///
168    /// ```text
169    /// x x x x x
170    /// 0 x x x x
171    /// 0 0 0 x x
172    /// 0 0 0 0 x
173    /// ```
174    Full,
175}
176
177/// Online QR decomposition using arbitrary orthogonalizer
178pub fn qr<A, S>(
179    iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
180    mut ortho: impl Orthogonalizer<Elem = A>,
181    strategy: Strategy,
182) -> (Q<A>, R<A>)
183where
184    A: Scalar + Lapack,
185    S: Data<Elem = A>,
186{
187    assert_eq!(ortho.len(), 0);
188
189    let mut coefs = Vec::new();
190    for a in iter {
191        match ortho.append(a.into_owned()) {
192            AppendResult::Added(coef) => coefs.push(coef),
193            AppendResult::Dependent(coef) => match strategy {
194                Strategy::Terminate => break,
195                Strategy::Skip => continue,
196                Strategy::Full => coefs.push(coef),
197            },
198        }
199    }
200    let n = ortho.len();
201    let m = coefs.len();
202    let mut r = Array2::zeros((n, m).f());
203    for j in 0..m {
204        for i in 0..n {
205            if i < coefs[j].len() {
206                r[(i, j)] = coefs[j][i];
207            }
208        }
209    }
210    (ortho.get_q(), r)
211}
ndarray_linalg/krylov/mod.rs

ndarray_linalg/krylov/
mod.rs
