mathhook_core/matrices/decomposition/
qr.rs

1//! QR decomposition algorithms
2//!
3//! This module provides QR decomposition using the Gram-Schmidt process
4//! for orthogonalization and solving least squares problems.
5
6use crate::core::Expression;
7use crate::matrices::types::*;
8use crate::matrices::unified::Matrix;
9use crate::simplify::Simplify;
10
11/// QR decomposition implementation
12impl Matrix {
13    /// Perform QR decomposition using Gram-Schmidt process
14    ///
15    /// Decomposes matrix A into A = QR where:
16    /// - Q is orthogonal (Q^T * Q = I)
17    /// - R is upper triangular
18    ///
19    /// # Examples
20    ///
21    /// ```
22    /// use mathhook_core::matrices::Matrix;
23    ///
24    /// let matrix = Matrix::from_arrays([
25    ///     [1, 1],
26    ///     [0, 1]
27    /// ]);
28    ///
29    /// let qr = matrix.qr_decomposition().unwrap();
30    /// let (q_rows, q_cols) = qr.q.dimensions();
31    /// assert_eq!(q_rows, 2);
32    /// assert_eq!(q_cols, 2);
33    /// ```
34    pub fn qr_decomposition(&self) -> Option<QRDecomposition> {
35        match self {
36            Matrix::Identity(data) => Some(QRDecomposition {
37                q: Matrix::identity(data.size),
38                r: Matrix::identity(data.size),
39            }),
40            Matrix::Zero(data) => Some(QRDecomposition {
41                q: Matrix::identity(data.rows),
42                r: Matrix::zero(data.rows, data.cols),
43            }),
44            _ => {
45                // General QR decomposition using Gram-Schmidt process
46                self.gram_schmidt_qr()
47            }
48        }
49    }
50
51    /// Gram-Schmidt QR decomposition implementation
52    fn gram_schmidt_qr(&self) -> Option<QRDecomposition> {
53        let (rows, cols) = self.dimensions();
54        let mut q_columns: Vec<Vec<Expression>> = Vec::new();
55        let mut r_elements = vec![vec![Expression::integer(0); cols]; cols];
56
57        // Convert columns to vectors for processing
58        for j in 0..cols {
59            let mut column: Vec<Expression> = (0..rows).map(|i| self.get_element(i, j)).collect();
60
61            // Orthogonalize against previous columns
62            for k in 0..j {
63                let dot_product = self.vector_dot(&column, &q_columns[k]);
64                r_elements[k][j] = dot_product.clone();
65
66                // column = column - dot_product * q_k
67                for elem in column.iter_mut().zip(&q_columns[k]) {
68                    let (col_elem, q_elem) = elem;
69                    let old_val = col_elem.clone();
70                    let subtract_val = Expression::mul(vec![dot_product.clone(), q_elem.clone()]);
71                    *col_elem = Expression::add(vec![
72                        old_val,
73                        Expression::mul(vec![Expression::integer(-1), subtract_val]),
74                    ])
75                    .simplify();
76                }
77            }
78
79            // Normalize the column
80            let norm = self.vector_norm(&column);
81            if norm.is_zero() {
82                return None; // Linearly dependent columns
83            }
84
85            r_elements[j][j] = norm.clone();
86
87            // Normalize column to get q_j
88            let q_column: Vec<Expression> = column
89                .iter()
90                .map(|elem| {
91                    // Use canonical form for division: a / b = a * b^(-1)
92                    Expression::mul(vec![
93                        elem.clone(),
94                        Expression::pow(norm.clone(), Expression::integer(-1)),
95                    ])
96                    .simplify()
97                })
98                .collect();
99            q_columns.push(q_column);
100        }
101
102        // Build Q and R matrices
103        let q_rows: Vec<Vec<Expression>> = (0..rows)
104            .map(|i| {
105                (0..cols)
106                    .map(|j| {
107                        if j < q_columns.len() {
108                            q_columns[j][i].clone()
109                        } else {
110                            Expression::integer(0)
111                        }
112                    })
113                    .collect()
114            })
115            .collect();
116
117        Some(QRDecomposition {
118            q: Matrix::dense(q_rows),
119            r: Matrix::dense(r_elements),
120        })
121    }
122
123    /// Compute dot product of two vectors (internal helper for QR decomposition)
124    ///
125    /// This is a private helper method. Result: v1 · v2 = sum(v1[i] * v2[i])
126    fn vector_dot(&self, v1: &[Expression], v2: &[Expression]) -> Expression {
127        if v1.len() != v2.len() {
128            return Expression::integer(0);
129        }
130
131        let products: Vec<Expression> = v1
132            .iter()
133            .zip(v2.iter())
134            .map(|(a, b)| Expression::mul(vec![a.clone(), b.clone()]))
135            .collect();
136
137        Expression::add(products).simplify()
138    }
139
140    /// Compute norm of a vector (internal helper for QR decomposition)
141    ///
142    /// This is a private helper method. Result: ||v|| = sqrt(sum(v[i]²))
143    fn vector_norm(&self, v: &[Expression]) -> Expression {
144        let sum_of_squares: Vec<Expression> = v
145            .iter()
146            .map(|x| Expression::pow(x.clone(), Expression::integer(2)))
147            .collect();
148
149        let sum = Expression::add(sum_of_squares).simplify();
150        Expression::pow(sum, Expression::rational(1, 2))
151    }
152}
mathhook_core/matrices/decomposition/qr.rs

mathhook_core/matrices/decomposition/
qr.rs
