mdarray-linalg 0.2.0

Linear algebra operations for mdarray, with multiple exchangeable backends
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100

use approx::assert_relative_eq;
use mdarray::DArray;
use num_complex::Complex;
use rand::prelude::*;

use super::common::{assert_complex_matrix_eq, assert_matrix_eq, naive_matmul};
use crate::{qr::QR, utils::pretty_print};

pub fn test_qr_random_matrix(bd: &impl QR<f64, usize, usize>) {
    let (m, n) = (5, 5);
    let mut rng = rand::rng();

    let a = DArray::<f64, 2>::from_fn([m, n], |_| rng.random::<f64>());
    test_qr_reconstruction(bd, &a);
}

pub fn test_qr_structured_matrix(bd: &impl QR<f64, usize, usize>) {
    let (m, n) = (3, 3);

    let a = DArray::<f64, 2>::from_fn([m, n], |i| (i[0] * i[1] + 1) as f64);
    test_qr_reconstruction(bd, &a);
}

pub fn test_qr_complex_matrix(bd: &impl QR<Complex<f64>, usize, usize>) {
    let (m, n) = (3, 3);

    let mut a = DArray::<Complex<f64>, 2>::from_fn([m, n], |i| {
        Complex::new((i[0] + 1) as f64, (i[1] + 1) as f64)
    });

    a[[1, 2]] = Complex::new(1., 5.); // destroy symmetry

    let mut q = DArray::<Complex<f64>, 2>::zeros([m, m]);
    let mut r = DArray::<Complex<f64>, 2>::zeros([m, n]);

    bd.qr_write(&mut a.clone(), &mut q, &mut r);
    let reconstructed = naive_matmul(&q, &r);
    assert_complex_matrix_eq!(a, reconstructed);

    let (q, r) = bd.qr(&mut a.clone());
    let reconstructed = naive_matmul(&q, &r);
    assert_complex_matrix_eq!(a, reconstructed);

    pretty_print(&a);
    pretty_print(&reconstructed);
}

pub fn test_qr_rectangular_matrix(bd: &impl QR<f64, usize, usize>) {
    // M > N
    {
        let (m, n) = (5, 3);
        let mut rng = rand::rng();
        let a = DArray::<f64, 2>::from_fn([m, n], |_| rng.random::<f64>());
        test_qr_reconstruction(bd, &a);
    }

    // M < N
    {
        let (m, n) = (3, 5);
        let mut rng = rand::rng();
        let a = DArray::<f64, 2>::from_fn([m, n], |_| rng.random::<f64>());
        test_qr_reconstruction(bd, &a);
    }
}

pub fn test_qr_reconstruction<T>(bd: &impl QR<T, usize, usize>, a: &DArray<T, 2>)
where
    T: num_traits::float::FloatConst
        + Default
        + Copy
        + std::fmt::Debug
        + approx::AbsDiffEq<Epsilon = f64>
        + std::fmt::Display
        + approx::RelativeEq
        + num_traits::Float
        + std::convert::From<i8>,
{
    let (m, n) = *a.shape();
    let mut q = DArray::<T, 2>::zeros([m, m]);
    let mut r = DArray::<T, 2>::zeros([m, n]);

    bd.qr_write(&mut a.clone(), &mut q, &mut r);
    let reconstructed = naive_matmul(&q, &r);

    pretty_print(&q);
    pretty_print(&r);

    pretty_print(a);
    pretty_print(&reconstructed);

    assert_matrix_eq!(a, reconstructed);

    let (q, r) = bd.qr(&mut a.clone());
    let reconstructed = naive_matmul(&q, &r);

    pretty_print(a);
    pretty_print(&reconstructed);

    assert_matrix_eq!(a, reconstructed);
}