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, r_column) in r_elements.iter_mut().enumerate() {
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, q_col) in q_columns.iter().enumerate() {
63                let dot_product = self.vector_dot(&column, q_col);
64                r_column[k] = dot_product.clone();
65
66                // column = column - dot_product * q_k
67                for (col_elem, q_elem) in column.iter_mut().zip(q_col) {
68                    let old_val = col_elem.clone();
69                    let subtract_val = Expression::mul(vec![dot_product.clone(), q_elem.clone()]);
70                    *col_elem = Expression::add(vec![
71                        old_val,
72                        Expression::mul(vec![Expression::integer(-1), subtract_val]),
73                    ])
74                    .simplify();
75                }
76            }
77
78            // Normalize the column
79            let norm = self.vector_norm(&column);
80            if norm.is_zero() {
81                return None; // Linearly dependent columns
82            }
83
84            r_column[j] = norm.clone();
85
86            // Normalize column to get q_j
87            let q_column: Vec<Expression> = column
88                .iter()
89                .map(|elem| {
90                    // Use canonical form for division: a / b = a * b^(-1)
91                    Expression::mul(vec![
92                        elem.clone(),
93                        Expression::pow(norm.clone(), Expression::integer(-1)),
94                    ])
95                    .simplify()
96                })
97                .collect();
98            q_columns.push(q_column);
99        }
100
101        // Build Q and R matrices
102        let q_rows: Vec<Vec<Expression>> = (0..rows)
103            .map(|i| {
104                (0..cols)
105                    .map(|j| {
106                        if j < q_columns.len() {
107                            q_columns[j][i].clone()
108                        } else {
109                            Expression::integer(0)
110                        }
111                    })
112                    .collect()
113            })
114            .collect();
115
116        Some(QRDecomposition {
117            q: Matrix::dense(q_rows),
118            r: Matrix::dense(r_elements),
119        })
120    }
121
122    /// Compute dot product of two vectors (internal helper for QR decomposition)
123    ///
124    /// This is a private helper method. Result: v1 · v2 = sum(v1[i] * v2[i])
125    fn vector_dot(&self, v1: &[Expression], v2: &[Expression]) -> Expression {
126        if v1.len() != v2.len() {
127            return Expression::integer(0);
128        }
129
130        let products: Vec<Expression> = v1
131            .iter()
132            .zip(v2.iter())
133            .map(|(a, b)| Expression::mul(vec![a.clone(), b.clone()]))
134            .collect();
135
136        Expression::add(products).simplify()
137    }
138
139    /// Compute norm of a vector (internal helper for QR decomposition)
140    ///
141    /// This is a private helper method. Result: ||v|| = sqrt(sum(v[i]²))
142    fn vector_norm(&self, v: &[Expression]) -> Expression {
143        let sum_of_squares: Vec<Expression> = v
144            .iter()
145            .map(|x| Expression::pow(x.clone(), Expression::integer(2)))
146            .collect();
147
148        let sum = Expression::add(sum_of_squares).simplify();
149        Expression::pow(sum, Expression::rational(1, 2))
150    }
151}
mathhook_core/matrices/decomposition/qr.rs

mathhook_core/matrices/decomposition/
qr.rs
