sciforge_lib/maths/linalg/
factorization.rs1pub fn qr_gram_schmidt(a: &[Vec<f64>]) -> (Vec<Vec<f64>>, Vec<Vec<f64>>) {
2 let m = a.len();
3 let n = if m > 0 { a[0].len() } else { 0 };
4 let mut q = vec![vec![0.0; m]; n];
5 let mut r = vec![vec![0.0; n]; n];
6 for j in 0..n {
7 let mut v: Vec<f64> = (0..m).map(|i| a[i][j]).collect();
8 for i in 0..j {
9 let dot: f64 = (0..m).map(|k| q[i][k] * v[k]).sum();
10 r[i][j] = dot;
11 for (vk, &qik) in v.iter_mut().zip(q[i].iter()) {
12 *vk -= dot * qik;
13 }
14 }
15 let norm = v.iter().map(|x| x * x).sum::<f64>().sqrt();
16 r[j][j] = norm;
17 if norm > 1e-30 {
18 for (qjk, &vk) in q[j].iter_mut().zip(v.iter()) {
19 *qjk = vk / norm;
20 }
21 }
22 }
23 let q_matrix = (0..m).map(|i| (0..n).map(|j| q[j][i]).collect()).collect();
24 (q_matrix, r)
25}
26
27pub fn qr_solve(a: &[Vec<f64>], b: &[f64]) -> Vec<f64> {
28 let (q, r) = qr_gram_schmidt(a);
29 let m = q.len();
30 let n = r.len();
31 let mut qtb = vec![0.0; n];
32 for (j, qtb_j) in qtb.iter_mut().enumerate() {
33 for i in 0..m {
34 *qtb_j += q[i][j] * b[i];
35 }
36 }
37 let mut x = vec![0.0; n];
38 for i in (0..n).rev() {
39 let mut sum = 0.0;
40 for (&rij, &xj) in r[i][(i + 1)..].iter().zip(&x[(i + 1)..]) {
41 sum += rij * xj;
42 }
43 x[i] = (qtb[i] - sum) / r[i][i];
44 }
45 x
46}
47
48pub fn householder_reflection(v: &[f64]) -> Vec<Vec<f64>> {
49 let n = v.len();
50 let dot: f64 = v.iter().map(|x| x * x).sum();
51 let mut h = vec![vec![0.0; n]; n];
52 for i in 0..n {
53 for j in 0..n {
54 h[i][j] = if i == j { 1.0 } else { 0.0 } - 2.0 * v[i] * v[j] / dot.max(1e-30);
55 }
56 }
57 h
58}
59
60pub fn rank(a: &[Vec<f64>], tol: f64) -> usize {
61 let m = a.len();
62 let n = if m > 0 {
63 a[0].len()
64 } else {
65 return 0;
66 };
67 let mut work: Vec<Vec<f64>> = a.to_vec();
68 let mut r = 0;
69 for col in 0..n {
70 let mut pivot = r;
71 for row in (r + 1)..m {
72 if work[row][col].abs() > work[pivot][col].abs() {
73 pivot = row;
74 }
75 }
76 if work[pivot][col].abs() < tol {
77 continue;
78 }
79 work.swap(r, pivot);
80 let div = work[r][col];
81 work[r][col..n].iter_mut().for_each(|v| *v /= div);
82 for row in 0..m {
83 if row == r {
84 continue;
85 }
86 let factor = work[row][col];
87 let src = work[r][col..n].to_vec();
88 work[row][col..n]
89 .iter_mut()
90 .zip(&src)
91 .for_each(|(d, &s)| *d -= factor * s);
92 }
93 r += 1;
94 }
95 r
96}
97
98pub fn null_space_dimension(a: &[Vec<f64>], tol: f64) -> usize {
99 let n = if a.is_empty() { 0 } else { a[0].len() };
100 n - rank(a, tol)
101}
102
103pub fn matrix_norm_frobenius(a: &[Vec<f64>]) -> f64 {
104 a.iter()
105 .flat_map(|row| row.iter())
106 .map(|x| x * x)
107 .sum::<f64>()
108 .sqrt()
109}
110
111pub fn matrix_norm_inf(a: &[Vec<f64>]) -> f64 {
112 a.iter()
113 .map(|row| row.iter().map(|x| x.abs()).sum::<f64>())
114 .fold(0.0_f64, f64::max)
115}
116
117pub fn is_symmetric(a: &[Vec<f64>], tol: f64) -> bool {
118 for (i, ai) in a.iter().enumerate() {
119 for (j, aj) in a.iter().enumerate().skip(i + 1) {
120 if (ai[j] - aj[i]).abs() > tol {
121 return false;
122 }
123 }
124 }
125 true
126}
127
128pub fn is_positive_definite(a: &[Vec<f64>]) -> bool {
129 let n = a.len();
130 let mut l = vec![vec![0.0; n]; n];
131 for i in 0..n {
132 for j in 0..=i {
133 let mut sum = 0.0;
134 for (&a, &b) in l[i][..j].iter().zip(&l[j][..j]) {
135 sum += a * b;
136 }
137 if i == j {
138 let val = a[i][i] - sum;
139 if val <= 0.0 {
140 return false;
141 }
142 l[i][j] = val.sqrt();
143 } else {
144 l[i][j] = (a[i][j] - sum) / l[j][j];
145 }
146 }
147 }
148 true
149}
150
151pub fn is_orthogonal(a: &[Vec<f64>], tol: f64) -> bool {
152 let n = a.len();
153 for i in 0..n {
154 for j in 0..n {
155 let mut dot = 0.0;
156 for ak in a.iter() {
157 dot += ak[i] * ak[j];
158 }
159 let expected = if i == j { 1.0 } else { 0.0 };
160 if (dot - expected).abs() > tol {
161 return false;
162 }
163 }
164 }
165 true
166}
167
168pub fn is_diagonal(a: &[Vec<f64>], tol: f64) -> bool {
169 for (i, ai) in a.iter().enumerate() {
170 for (j, &aij) in ai.iter().enumerate() {
171 if i != j && aij.abs() > tol {
172 return false;
173 }
174 }
175 }
176 true
177}
178
179pub fn is_upper_triangular(a: &[Vec<f64>], tol: f64) -> bool {
180 for (i, ai) in a.iter().enumerate() {
181 for &aij in &ai[..i] {
182 if aij.abs() > tol {
183 return false;
184 }
185 }
186 }
187 true
188}
189
190pub fn is_lower_triangular(a: &[Vec<f64>], tol: f64) -> bool {
191 for (i, ai) in a.iter().enumerate() {
192 for &aij in &ai[(i + 1)..] {
193 if aij.abs() > tol {
194 return false;
195 }
196 }
197 }
198 true
199}
200
201pub fn matrix_norm_1(a: &[Vec<f64>]) -> f64 {
202 let n = a.len();
203 let m = if n > 0 { a[0].len() } else { 0 };
204 let mut max_col = 0.0f64;
205 for j in 0..m {
206 let col_sum: f64 = a.iter().map(|row| row[j].abs()).sum();
207 if col_sum > max_col {
208 max_col = col_sum;
209 }
210 }
211 max_col
212}
213
214pub fn matrix_condition_frobenius(a: &[Vec<f64>]) -> f64 {
215 let n = a.len();
216 let mut inv = vec![vec![0.0; n]; n];
217 let mut aug: Vec<Vec<f64>> = (0..n)
218 .map(|i| {
219 let mut row = a[i].clone();
220 row.push(0.0);
221 row
222 })
223 .collect();
224 for col_idx in 0..n {
225 for (i, aug_i) in aug.iter_mut().enumerate() {
226 aug_i[n] = if i == col_idx { 1.0 } else { 0.0 };
227 }
228 let mut work = aug.clone();
229 for c in 0..n {
230 let mut pivot = c;
231 for r in (c + 1)..n {
232 if work[r][c].abs() > work[pivot][c].abs() {
233 pivot = r;
234 }
235 }
236 work.swap(c, pivot);
237 if work[c][c].abs() < 1e-30 {
238 continue;
239 }
240 for r in (c + 1)..n {
241 let f = work[r][c] / work[c][c];
242 let (top, bot) = work.split_at_mut(r);
243 for (d, &s) in bot[0][c..=n].iter_mut().zip(&top[c][c..=n]) {
244 *d -= f * s;
245 }
246 }
247 }
248 for i in (0..n).rev() {
249 let mut s = 0.0;
250 for (&wij, inv_j) in work[i][(i + 1)..].iter().zip(&inv[(i + 1)..]) {
251 s += wij * inv_j[col_idx];
252 }
253 inv[i][col_idx] = if work[i][i].abs() > 1e-30 {
254 (work[i][n] - s) / work[i][i]
255 } else {
256 0.0
257 };
258 }
259 }
260 let norm_a = matrix_norm_frobenius(a);
261 let norm_inv = matrix_norm_frobenius(&inv);
262 norm_a * norm_inv
263}
264
265pub fn gram_matrix(a: &[Vec<f64>]) -> Vec<Vec<f64>> {
266 let m = a.len();
267 let mut g = vec![vec![0.0; m]; m];
268 for i in 0..m {
269 for j in i..m {
270 let mut dot = 0.0;
271 for (&ai, &aj) in a[i].iter().zip(&a[j]) {
272 dot += ai * aj;
273 }
274 g[i][j] = dot;
275 g[j][i] = dot;
276 }
277 }
278 g
279}
280
281pub fn column_space_basis(a: &[Vec<f64>], tol: f64) -> Vec<Vec<f64>> {
282 let m = a.len();
283 let n = if m > 0 {
284 a[0].len()
285 } else {
286 return Vec::new();
287 };
288 let mut work = a.to_vec();
289 let mut pivot_cols = Vec::new();
290 let mut r = 0;
291 for col in 0..n {
292 let mut pivot = r;
293 for row in (r + 1)..m {
294 if work[row][col].abs() > work[pivot][col].abs() {
295 pivot = row;
296 }
297 }
298 if work[pivot][col].abs() < tol {
299 continue;
300 }
301 work.swap(r, pivot);
302 let div = work[r][col];
303 work[r][col..n].iter_mut().for_each(|v| *v /= div);
304 for row in 0..m {
305 if row == r {
306 continue;
307 }
308 let f = work[row][col];
309 let src = work[r][col..n].to_vec();
310 work[row][col..n]
311 .iter_mut()
312 .zip(&src)
313 .for_each(|(d, &s)| *d -= f * s);
314 }
315 pivot_cols.push(col);
316 r += 1;
317 }
318 pivot_cols
319 .iter()
320 .map(|&col| (0..m).map(|i| a[i][col]).collect())
321 .collect()
322}