algebraeon_rings/matrix/
gram_schmidt.rs

1use super::*;
2
3// #[derive(Debug, Clone)]
4// struct FiniteDimensionalInnerProductSpaceStructure<Ring: ComplexSubsetStructure> {
5//     ring: Ring,
6//     rows: usize,
7//     cols: usize,
8//     // Let n = self.rows * self.cols.
9//     // Then self.values is an n by n matrix containing the inner products
10//     // of the elementary basis vectors when expressed as metamatricies.
11//     values: Matrix<Ring::Set>,
12// }
13
14// impl<Ring: ComplexSubsetStructure> PartialEq for FiniteDimensionalInnerProductSpaceStructure<Ring> {
15//     fn eq(&self, other: &Self) -> bool {
16//         self.ring == other.ring
17//             && self.rows == other.rows
18//             && self.cols == other.cols
19//             && MatrixStructure::new(self.ring.clone()).equal(&self.values, &other.values)
20//     }
21// }
22
23// impl<Ring: ComplexSubsetStructure> Eq for FiniteDimensionalInnerProductSpaceStructure<Ring> {}
24
25impl<FS: ComplexConjugateSignature, FSB: BorrowedStructure<FS>> MatrixStructure<FS, FSB> {
26    pub fn conjugate(&self, mat: &Matrix<FS::Set>) -> Matrix<FS::Set> {
27        mat.apply_map(|x| self.ring().conjugate(x))
28    }
29
30    pub fn conjugate_transpose(&self, mat: &Matrix<FS::Set>) -> Matrix<FS::Set> {
31        self.conjugate(mat).transpose()
32    }
33}
34
35impl<FS: ComplexConjugateSignature + RingSignature, FSB: BorrowedStructure<FS>>
36    MatrixStructure<FS, FSB>
37{
38    // dot product of a and conj(b)
39    pub fn inner_product(&self, a: &Matrix<FS::Set>, b: &Matrix<FS::Set>) -> FS::Set {
40        self.dot(a, &self.conjugate(b))
41    }
42}
43
44impl<FS: ComplexConjugateSignature + FieldSignature, FSB: BorrowedStructure<FS>>
45    MatrixStructure<FS, FSB>
46{
47    //return mat=LQ where L is lower triangular and Q is row-orthogonal (not orthonormal)
48    pub fn gram_schmidt_row_orthogonalization_algorithm(
49        &self,
50        mut mat: Matrix<FS::Set>,
51    ) -> (Matrix<FS::Set>, Matrix<FS::Set>) {
52        #[cfg(debug_assertions)]
53        let original_mat = mat.clone();
54
55        let mut lt = self.ident(mat.rows());
56        for i in 0..mat.rows() {
57            for j in 0..i {
58                //col_i = col_i - (projection of col_j onto col_i)
59                let lambda = self
60                    .ring()
61                    .try_divide(
62                        &self.inner_product(&mat.get_row_submatrix(i), &mat.get_row_submatrix(j)),
63                        &self.inner_product(&mat.get_row_submatrix(j), &mat.get_row_submatrix(j)),
64                    )
65                    .unwrap();
66                //col_i -= lambda col_j
67                let row_opp = ElementaryOpp::new_row_opp(
68                    self.ring().clone(),
69                    ElementaryOppType::AddRowMul {
70                        i,
71                        j,
72                        x: self.ring().neg(&lambda),
73                    },
74                );
75                row_opp.apply(&mut lt);
76                row_opp.apply(&mut mat);
77            }
78        }
79
80        for i in 0..mat.rows() {
81            for j in (i + 1)..mat.rows() {
82                debug_assert!(self.ring().is_zero(
83                    &self.inner_product(&mat.get_row_submatrix(i), &mat.get_row_submatrix(j)),
84                ));
85            }
86        }
87
88        #[cfg(debug_assertions)]
89        assert!(self.equal(&mat, &self.mul(&lt, &original_mat).unwrap()));
90
91        (lt, mat)
92    }
93
94    //return mat=QR where Q is col-orthogonal (not orthonormal) and R is upper triangular and
95    pub fn gram_schmidt_col_orthogonalization_algorithm(
96        &self,
97        mat: Matrix<FS::Set>,
98    ) -> (Matrix<FS::Set>, Matrix<FS::Set>) {
99        let (l, q) = self.gram_schmidt_row_orthogonalization_algorithm(mat.transpose());
100        (q.transpose(), l.transpose())
101    }
102
103    pub fn gram_schmidt_row_orthogonalization(&self, mat: Matrix<FS::Set>) -> Matrix<FS::Set> {
104        self.gram_schmidt_row_orthogonalization_algorithm(mat).1
105    }
106
107    pub fn gram_schmidt_col_orthogonalization(&self, mat: Matrix<FS::Set>) -> Matrix<FS::Set> {
108        self.gram_schmidt_col_orthogonalization_algorithm(mat).0
109    }
110}
111
112impl<
113    FS: ComplexConjugateSignature + PositiveRealNthRootSignature + FieldSignature,
114    FSB: BorrowedStructure<FS>,
115> MatrixStructure<FS, FSB>
116{
117    //return L*mat=Q where L is lower triangular and Q is orthonormal
118    pub fn lq_decomposition_algorithm(
119        &self,
120        mat: Matrix<FS::Set>,
121    ) -> (Matrix<FS::Set>, Matrix<FS::Set>) {
122        let (mut lt, mut mat) = self.gram_schmidt_row_orthogonalization_algorithm(mat);
123
124        for r in 0..mat.rows() {
125            let row = mat.get_row_submatrix(r);
126            let lensq = self.inner_product(&row, &row);
127            let row_opp = ElementaryOpp::new_row_opp(
128                self.ring().clone(),
129                ElementaryOppType::UnitMul {
130                    row: r,
131                    unit: self
132                        .ring()
133                        .try_reciprocal(&self.ring().nth_root(&lensq, 2).unwrap())
134                        .unwrap(),
135                },
136            );
137            row_opp.apply(&mut lt);
138            row_opp.apply(&mut mat);
139        }
140
141        debug_assert!(
142            self.equal(
143                &self.ident(mat.rows()),
144                &self
145                    .mul(&mat, &self.conjugate(&mat.transpose_ref()))
146                    .unwrap()
147            )
148        );
149
150        (lt, mat)
151    }
152
153    //return mat*R=Q where Q is col-orthogonal (not orthonormal) and R is upper triangular
154    pub fn qr_decomposition_algorithm(
155        &self,
156        mat: Matrix<FS::Set>,
157    ) -> (Matrix<FS::Set>, Matrix<FS::Set>) {
158        let (l, q) = self.lq_decomposition_algorithm(mat.transpose());
159        (q.transpose(), l.transpose())
160    }
161
162    pub fn gram_schmidt_row_orthonormalization(&self, mat: Matrix<FS::Set>) -> Matrix<FS::Set> {
163        self.lq_decomposition_algorithm(mat).1
164    }
165
166    pub fn gram_schmidt_col_orthonormalization(&self, mat: Matrix<FS::Set>) -> Matrix<FS::Set> {
167        self.qr_decomposition_algorithm(mat).0
168    }
169}
170
171impl<F: MetaType> Matrix<F>
172where
173    F::Signature: ComplexConjugateSignature + FieldSignature,
174{
175    pub fn gram_schmidt_row_orthogonalization_algorithm(self) -> (Matrix<F>, Matrix<F>) {
176        Self::structure().gram_schmidt_row_orthogonalization_algorithm(self)
177    }
178
179    pub fn gram_schmidt_col_orthogonalization_algorithm(self) -> (Matrix<F>, Matrix<F>) {
180        Self::structure().gram_schmidt_col_orthogonalization_algorithm(self)
181    }
182
183    pub fn gram_schmidt_row_orthogonalization(self) -> Matrix<F> {
184        Self::structure().gram_schmidt_row_orthogonalization(self)
185    }
186
187    pub fn gram_schmidt_col_orthogonalization(self) -> Matrix<F> {
188        Self::structure().gram_schmidt_col_orthogonalization(self)
189    }
190}
191
192impl<F: MetaType> Matrix<F>
193where
194    F::Signature: ComplexConjugateSignature + PositiveRealNthRootSignature + FieldSignature,
195{
196    pub fn lq_decomposition_algorithm(self) -> (Matrix<F>, Matrix<F>) {
197        Self::structure().lq_decomposition_algorithm(self)
198    }
199
200    pub fn qr_decomposition_algorithm(self) -> (Matrix<F>, Matrix<F>) {
201        Self::structure().qr_decomposition_algorithm(self)
202    }
203
204    pub fn gram_schmidt_row_orthonormalization(self) -> Matrix<F> {
205        Self::structure().gram_schmidt_row_orthonormalization(self)
206    }
207
208    pub fn gram_schmidt_col_orthonormalization(self) -> Matrix<F> {
209        Self::structure().gram_schmidt_col_orthonormalization(self)
210    }
211}
212
213#[cfg(test)]
214mod tests {
215    use super::*;
216    use crate::isolated_algebraic::{ComplexAlgebraic, RealAlgebraic};
217    use std::str::FromStr;
218
219    #[test]
220    fn rational_gram_schmidt() {
221        let mat = Matrix::<Rational>::from_rows(vec![
222            vec![Rational::from(1), Rational::from(-1), Rational::from(3)],
223            vec![Rational::from(1), Rational::from(0), Rational::from(5)],
224            vec![Rational::from(1), Rational::from(2), Rational::from(6)],
225        ]);
226
227        let mat_expected_gs = Matrix::from_rows(vec![
228            vec![
229                Rational::from_str("1").unwrap(),
230                Rational::from_str("-4/3").unwrap(),
231                Rational::from_str("-3/7").unwrap(),
232            ],
233            vec![
234                Rational::from_str("1").unwrap(),
235                Rational::from_str("-1/3").unwrap(),
236                Rational::from_str("9/14").unwrap(),
237            ],
238            vec![
239                Rational::from_str("1").unwrap(),
240                Rational::from_str("5/3").unwrap(),
241                Rational::from_str("-3/14").unwrap(),
242            ],
243        ]);
244
245        mat.pprint();
246        mat_expected_gs.pprint();
247
248        let (mat_actual_gs, mat_actual_ut) =
249            mat.clone().gram_schmidt_col_orthogonalization_algorithm();
250
251        mat_actual_gs.pprint();
252        mat_actual_ut.pprint();
253        Matrix::mul(&mat, &mat_actual_ut).unwrap().pprint();
254
255        assert_eq!(mat.gram_schmidt_col_orthogonalization(), mat_expected_gs);
256    }
257
258    #[allow(clippy::erasing_op)]
259    #[test]
260    fn complex_gram_schmidt() {
261        let i = &ComplexAlgebraic::i().into_ergonomic();
262
263        let mat = Matrix::<ComplexAlgebraic>::from_rows(vec![
264            vec![(1 + 0 * i).into_verbose(), (1 * i).into_verbose()],
265            vec![(1 + 0 * i).into_verbose(), (1 + 0 * i).into_verbose()],
266        ]);
267        mat.pprint();
268        mat.gram_schmidt_col_orthogonalization().pprint();
269
270        let mat = Matrix::<ComplexAlgebraic>::from_rows(vec![
271            vec![
272                (-2 + 2 * i).into_verbose(),
273                (7 + 3 * i).into_verbose(),
274                (7 + 3 * i).into_verbose(),
275            ],
276            vec![
277                (3 + 3 * i).into_verbose(),
278                (-2 + 4 * i).into_verbose(),
279                (6 + 2 * i).into_verbose(),
280            ],
281            vec![
282                (2 + 2 * i).into_verbose(),
283                (8 + 0 * i).into_verbose(),
284                (-9 + 1 * i).into_verbose(),
285            ],
286        ]);
287        mat.pprint();
288        mat.clone().gram_schmidt_col_orthogonalization().pprint();
289    }
290
291    #[test]
292    fn complex_gram_schmidt_normalized() {
293        let one = &RealAlgebraic::one().into_ergonomic();
294
295        let mat = Matrix::<RealAlgebraic>::from_rows(vec![
296            vec![(1 * one).into_verbose(), (1 * one).into_verbose()],
297            vec![(1 * one).into_verbose(), (2 * one).into_verbose()],
298        ]);
299        mat.pprint();
300        mat.gram_schmidt_col_orthonormalization().pprint();
301
302        let i = &ComplexAlgebraic::i().into_ergonomic();
303        let mat = Matrix::<ComplexAlgebraic>::from_rows(vec![
304            vec![(-2 + 2 * i).into_verbose(), (-9 + 1 * i).into_verbose()],
305            vec![(3 + 3 * i).into_verbose(), (-2 + 4 * i).into_verbose()],
306        ]);
307        mat.pprint();
308        mat.clone().gram_schmidt_col_orthonormalization().pprint();
309    }
310}
algebraeon_rings/matrix/gram_schmidt.rs

algebraeon_rings/matrix/
gram_schmidt.rs
