sciforge_lib/maths/tensor/
linalg.rs

1use super::ops::transpose;
2use super::storage::Tensor;
3
4pub fn trace(t: &Tensor) -> f64 {
5    assert!(t.rank() == 2 && t.shape()[0] == t.shape()[1]);
6    (0..t.shape()[0]).map(|i| t.get(&[i, i])).sum()
7}
8
9pub fn matmul(a: &Tensor, b: &Tensor) -> Tensor {
10    assert!(a.rank() == 2 && b.rank() == 2);
11    assert_eq!(a.shape()[1], b.shape()[0]);
12    let (m, n, p) = (a.shape()[0], a.shape()[1], b.shape()[1]);
13    Tensor::from_fn(&[m, p], |idx| {
14        (0..n)
15            .map(|k| a.get(&[idx[0], k]) * b.get(&[k, idx[1]]))
16            .sum()
17    })
18}
19
20pub fn minor(t: &Tensor, row: usize, col: usize) -> Tensor {
21    let n = t.shape()[0];
22    let mut data = Vec::with_capacity((n - 1) * (n - 1));
23    for i in 0..n {
24        if i == row {
25            continue;
26        }
27        for j in 0..n {
28            if j == col {
29                continue;
30            }
31            data.push(t.get(&[i, j]));
32        }
33    }
34    Tensor::from_vec(&[n - 1, n - 1], data)
35}
36
37pub fn determinant(t: &Tensor) -> f64 {
38    assert!(t.rank() == 2 && t.shape()[0] == t.shape()[1]);
39    let n = t.shape()[0];
40    if n == 1 {
41        return t.get(&[0, 0]);
42    }
43    if n == 2 {
44        return t.get(&[0, 0]) * t.get(&[1, 1]) - t.get(&[0, 1]) * t.get(&[1, 0]);
45    }
46    let mut det = 0.0;
47    for j in 0..n {
48        let sign = if j % 2 == 0 { 1.0 } else { -1.0 };
49        det += sign * t.get(&[0, j]) * determinant(&minor(t, 0, j));
50    }
51    det
52}
53
54pub fn inverse(t: &Tensor) -> Option<Tensor> {
55    assert!(t.rank() == 2 && t.shape()[0] == t.shape()[1]);
56    let n = t.shape()[0];
57    let det = determinant(t);
58    if det.abs() < 1e-15 {
59        return None;
60    }
61    let cofactors = Tensor::from_fn(&[n, n], |idx| {
62        let sign = if (idx[0] + idx[1]) % 2 == 0 {
63            1.0
64        } else {
65            -1.0
66        };
67        sign * determinant(&minor(t, idx[0], idx[1]))
68    });
69    let adjugate = transpose(&cofactors, &[1, 0]);
70    Some(adjugate.scale(1.0 / det))
71}
72
73pub fn eigenvalues_2x2(t: &Tensor) -> (f64, f64) {
74    assert!(t.rank() == 2 && t.shape() == [2, 2]);
75    let a = t.get(&[0, 0]);
76    let b = t.get(&[0, 1]);
77    let c = t.get(&[1, 0]);
78    let d = t.get(&[1, 1]);
79    let tr = a + d;
80    let det = a * d - b * c;
81    let disc = (tr * tr - 4.0 * det).max(0.0).sqrt();
82    ((tr + disc) / 2.0, (tr - disc) / 2.0)
83}
84
85pub fn lu_decompose(t: &Tensor) -> (Tensor, Tensor) {
86    assert!(t.rank() == 2 && t.shape()[0] == t.shape()[1]);
87    let n = t.shape()[0];
88    let mut l = Tensor::identity(n);
89    let mut u = t.clone();
90    for j in 0..n {
91        for i in (j + 1)..n {
92            let pivot = u.get(&[j, j]);
93            if pivot.abs() < 1e-30 {
94                continue;
95            }
96            let factor = u.get(&[i, j]) / pivot;
97            l.set(&[i, j], factor);
98            for k in j..n {
99                let val = u.get(&[i, k]) - factor * u.get(&[j, k]);
100                u.set(&[i, k], val);
101            }
102        }
103    }
104    (l, u)
105}
106
107pub fn solve(a: &Tensor, b: &[f64]) -> Vec<f64> {
108    let (l, u) = lu_decompose(a);
109    let n = a.shape()[0];
110    let mut y = vec![0.0; n];
111    for i in 0..n {
112        y[i] = b[i];
113        for j in 0..i {
114            y[i] -= l.get(&[i, j]) * y[j];
115        }
116    }
117    let mut x = vec![0.0; n];
118    for i in (0..n).rev() {
119        x[i] = y[i];
120        for j in (i + 1)..n {
121            x[i] -= u.get(&[i, j]) * x[j];
122        }
123        let pivot = u.get(&[i, i]);
124        if pivot.abs() > 1e-30 {
125            x[i] /= pivot;
126        }
127    }
128    x
129}
sciforge_lib/maths/tensor/linalg.rs

sciforge_lib/maths/tensor/
linalg.rs
