rusty_compression/
random_matrix.rs

1//! Generation of random matrices for various types
2
3use ndarray::Array2;
4use ndarray_linalg::{Lapack, SVDDCInto, Scalar, UVTFlag};
5use num::complex::Complex;
6use num::traits::cast::cast;
7use num::Float;
8use rand::Rng;
9use rand_distr::{Distribution, Normal};
10
11pub trait RandomMatrix
12where
13    Self: Scalar + Lapack,
14{
15    /// Generate a random Gaussian matrix.
16    ///
17    /// # Arguments
18    ///
19    /// * `dimension`: Tuple (rows, cols) specifying the number of rows and columns.
20    /// * `rng`: The random number generator to use.
21    fn random_gaussian<R: Rng>(dimension: (usize, usize), rng: &mut R) -> Array2<Self>;
22
23    /// Generate a random matrix with orthogonal rows or columns.
24    ///
25    /// This function creates a normally distributed (m, n) random matrix,
26    /// orthogonalizes it and returns the resulting orthogonal matrix.
27    ///
28    /// If m > n then the returned matrix has orthogonal columns. If n > m
29    /// the returned matrix has orthogonalized rows.
30    ///
31    /// # Arguments
32    ///
33    /// * `dimension`: Tuple (rows, cols) specifying the number of rows and columns.
34    /// * `rng`: The random number generator to use.
35    fn random_orthogonal_matrix<R: Rng>(dimension: (usize, usize), rng: &mut R) -> Array2<Self> {
36        let mut m = dimension.0;
37        let mut n = dimension.1;
38
39        // Always ensure that we form the QR decomp for a long and skinny matrix
40        if dimension.1 > dimension.0 {
41            std::mem::swap(&mut m, &mut n);
42        }
43
44        let mat = Self::random_gaussian((m, n), rng);
45
46        let (u, _, _) = mat
47            .svddc_into(UVTFlag::Some)
48            .expect("`compress_svd_rank_based`: SVD computation failed.");
49
50        // If we originally had more columns than rows, conjugate transpose again.
51        if dimension.1 > dimension.0 {
52            u.unwrap().t().map(|item| item.conj())
53        } else {
54            u.unwrap()
55        }
56    }
57
58    /// Generate a random approximate low-rank matrix.
59    ///
60    /// This function generates a random approximate low-rank matrix
61    /// with singular values logarithmically distributed between
62    /// 'sigma_max` and `sigma_min`.
63    ///
64    /// # Arguments
65    ///
66    /// * `dimension`: Tuple (rows, cols) specifying the number of rows and columns.
67    /// * `sigma_max`: Maximum singular value.
68    /// * `sigma_min`: Minimum singular value.
69    /// * `rng`: The random number generator to use.
70    fn random_approximate_low_rank_matrix<R: Rng>(
71        dimension: (usize, usize),
72        sigma_max: f64,
73        sigma_min: f64,
74        rng: &mut R,
75    ) -> Array2<Self> {
76        use ndarray::Array;
77
78        assert!(
79            sigma_min < sigma_max,
80            "`sigma_min` must be smaller than `sigma_max`"
81        );
82        assert!(sigma_min > 0.0, "`sigma_min` must be positive.");
83
84        let min_dim = std::cmp::min(dimension.0, dimension.1);
85
86        let u = Self::random_orthogonal_matrix((dimension.0, min_dim), rng);
87        let vt = Self::random_orthogonal_matrix((min_dim, dimension.1), rng);
88        let singvals = Array::geomspace(sigma_min, sigma_max, min_dim)
89            .unwrap()
90            .map(|&item| cast::<f64, Self>(item).unwrap());
91        let sigma = Array2::from_diag(&singvals);
92        u.dot(&sigma.dot(&vt))
93    }
94}
95
96impl RandomMatrix for f64 {
97    fn random_gaussian<R: Rng>(dimension: (usize, usize), rng: &mut R) -> Array2<f64> {
98        random_gaussian_real::<f64, R>(dimension, rng)
99    }
100}
101
102impl RandomMatrix for f32 {
103    fn random_gaussian<R: Rng>(dimension: (usize, usize), rng: &mut R) -> Array2<f32> {
104        random_gaussian_real::<f32, R>(dimension, rng)
105    }
106}
107
108impl RandomMatrix for Complex<f64> {
109    fn random_gaussian<R: Rng>(dimension: (usize, usize), rng: &mut R) -> Array2<Complex<f64>> {
110        random_gaussian_complex::<f64, R>(dimension, rng)
111    }
112}
113
114impl RandomMatrix for Complex<f32> {
115    fn random_gaussian<R: Rng>(dimension: (usize, usize), rng: &mut R) -> Array2<Complex<f32>> {
116        random_gaussian_complex::<f32, R>(dimension, rng)
117    }
118}
119
120fn random_gaussian_real<T: Float, R: Rng>(dimension: (usize, usize), rng: &mut R) -> Array2<T> {
121    let mut mat = Array2::<T>::zeros(dimension);
122    let normal = Normal::new(0.0, 1.0).unwrap();
123    mat.map_inplace(|item| *item = cast::<f64, T>(normal.sample(rng)).unwrap());
124    mat
125}
126
127/// Generate a random Gaussian matrix.
128///
129/// # Arguments
130///
131/// * `dimension`: Tuple (rows, cols) specifying the number of rows and columns.
132/// * `rng`: The random number generator to use.
133fn random_gaussian_complex<T: Float, R: Rng>(
134    dimension: (usize, usize),
135    rng: &mut R,
136) -> Array2<Complex<T>> {
137    let mut mat = Array2::<Complex<T>>::zeros(dimension);
138    let normal = Normal::new(0.0, 1.0).unwrap();
139    mat.map_inplace(|item| {
140        let re = cast::<f64, T>(normal.sample(rng)).unwrap();
141        let im = cast::<f64, T>(normal.sample(rng)).unwrap();
142        *item = Complex::new(re, im);
143    });
144    mat
145}
rusty_compression/random_matrix.rs

rusty_compression/
random_matrix.rs
