Implementation of new linear algebra traits

examples/rlst_interpolative_decomposition.rs

@@ -0,0 +1,125 @@
+//! Demo interpolative decomposition of a matrix.
+
+use rlst::dense::linalg::interpolative_decomposition::{Accuracy, MatrixIdDecomposition};
+pub use rlst::prelude::*;
+
+//Function that creates a low rank matrix by calculating a kernel given a random point distribution on an unit sphere.
+fn low_rank_matrix(n: usize, arr: &mut Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2>) {
+    //Obtain n equally distributed angles 0<phi<pi and n equally distributed angles 0<theta<2pi
+    let pi: f64 = std::f64::consts::PI;
+
+    let mut angles1 = rlst_dynamic_array2!(f64, [n, 1]);
+    angles1.fill_from_seed_equally_distributed(0);
+    angles1.scale_inplace(pi);
+
+    let mut angles2 = rlst_dynamic_array2!(f64, [n, 1]);
+    angles2.fill_from_seed_equally_distributed(1);
+    angles2.scale_inplace(2.0 * pi);
+
+    //Calculate n points on a sphere given phi and theta
+    let mut points = rlst_dynamic_array2!(f64, [n, 3]);
+
+    for i in 0..n {
+        let phi: f64 = angles1.get_value([i, 0]).unwrap();
+        let theta: f64 = angles2.get_value([i, 0]).unwrap();
+        *points.get_mut([i, 0]).unwrap() = phi.sin() * theta.cos();
+        *points.get_mut([i, 1]).unwrap() = phi.sin() * theta.sin();
+        *points.get_mut([i, 2]).unwrap() = phi.cos();
+    }
+
+    for i in 0..arr.shape()[0] {
+        let point1 = points.r().into_subview([i, 0], [1, 3]);
+        for j in 0..n {
+            if i != j {
+                let point2 = points.r().into_subview([j, 0], [1, 3]);
+                //Calculate distance between all combinations of different points.
+                let res = point1.r().sub(point2.r());
+                //Calculate kernel e^{-|points1-points2|^2}
+                *arr.get_mut([i, j]).unwrap() = 1.0 / ((res.view_flat().norm_2().pow(2.0)).exp());
+            } else {
+                //If points are equal, set the value to 1
+                *arr.get_mut([i, j]).unwrap() = 1.0;
+            }
+        }
+    }
+}
+
+fn test_fat_matrix(slice: usize, n: usize, tol: f64) {
+    let mut arr: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice, n]);
+    low_rank_matrix(n, &mut arr);
+
+    let res: IdDecomposition<f64> = arr
+        .r_mut()
+        .into_id_alloc(Accuracy::Tol(tol), TransMode::NoTrans)
+        .unwrap();
+
+    println!("The skeleton of the matrix is given by");
+    res.skel.pretty_print();
+
+    println!("The interpolative decomposition matrix is:");
+    res.id_mat.pretty_print();
+
+    println!("The rank of this matrix is {}\n", res.rank);
+
+    //We extract the residuals of the matrix
+    let mut perm_mat: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice, slice]);
+    res.get_p(perm_mat.r_mut());
+    let perm_arr: DynamicArray<f64, 2> =
+        empty_array::<f64, 2>().simple_mult_into_resize(perm_mat.r_mut(), arr.r());
+    let mut a_rs: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice - res.rank, n]);
+    a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice - res.rank, n]));
+
+    //We compute an approximation of the residual columns of the matrix
+    let a_rs_app: DynamicArray<f64, 2> =
+        empty_array().simple_mult_into_resize(res.id_mat.r(), res.skel);
+
+    let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2();
+    println!("Interpolative Decomposition L2 absolute error: {}", error);
+}
+
+fn test_skinny_matrix(slice: usize, n: usize, tol: f64) {
+    let mut arr: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice, n]);
+    low_rank_matrix(n, &mut arr);
+
+    let mut arr_trans = empty_array();
+    arr_trans.fill_from_resize(arr.r().transpose());
+
+    let res: IdDecomposition<f64> = arr_trans
+        .r_mut()
+        .into_id_alloc(Accuracy::Tol(tol), TransMode::Trans)
+        .unwrap();
+
+    println!("The skeleton of the matrix is given by");
+    res.skel.pretty_print();
+
+    println!("The interpolative decomposition matrix is:");
+    res.id_mat.pretty_print();
+
+    println!("The rank of this matrix is {}\n", res.rank);
+
+    //We extract the residuals of the matrix
+    let mut perm_mat: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice, slice]);
+    res.get_p(perm_mat.r_mut());
+    let perm_arr: DynamicArray<f64, 2> =
+        empty_array::<f64, 2>().simple_mult_into_resize(perm_mat.r_mut(), arr.r());
+    let mut a_rs: DynamicArray<f64, 2> = rlst_dynamic_array2!(f64, [slice - res.rank, n]);
+    a_rs.fill_from(perm_arr.into_subview([res.rank, 0], [slice - res.rank, n]));
+
+    //We compute an approximation of the residual columns of the matrix
+    let a_rs_app: DynamicArray<f64, 2> =
+        empty_array().simple_mult_into_resize(res.id_mat.r(), res.skel);
+
+    let error: f64 = a_rs.r().sub(a_rs_app.r()).view_flat().norm_2();
+    println!("Interpolative Decomposition L2 absolute error: {}", error);
+}
+
+pub fn main() {
+    let n: usize = 100;
+    let slice: usize = 50;
+    //Tolerance given for the
+    let tol: f64 = 1e-4;
+
+    //Create a low rank matrix
+    test_fat_matrix(slice, n, tol);
+    test_skinny_matrix(slice, n, tol);
+}

examples/rlst_null_space.rs

@@ -0,0 +1,32 @@
+//! Demo the null space of a matrix.
+use rlst::dense::linalg::null_space::Method;
+pub use rlst::prelude::*;
+
+//Here we compute the null space (B) of the rowspace of a short matrix (A).
+pub fn main() {
+    let mut arr = rlst_dynamic_array2!(f64, [3, 4]);
+    arr.fill_from_seed_equally_distributed(0);
+    let tol = 1e-15;
+    let null_res = arr.r_mut().into_null_alloc(tol, Method::Svd).unwrap();
+    let res: Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2> =
+        empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r());
+    //Result
+    println!(
+        "Nullification via SVD. Value of |A*B|_2, where B=null(A): {}",
+        res.view_flat().norm_2()
+    );
+
+    let mut arr = rlst_dynamic_array2!(f64, [4, 3]);
+    arr.fill_from_seed_equally_distributed(0);
+    let shape = arr.shape();
+    let mut arr_qr = rlst_dynamic_array2!(f64, [shape[1], shape[0]]);
+    arr_qr.fill_from(arr.r().transpose());
+    let null_res = arr_qr.r_mut().into_null_alloc(tol, Method::Qr).unwrap();
+    let res: Array<f64, BaseArray<f64, VectorContainer<f64>, 2>, 2> =
+        empty_array().simple_mult_into_resize(arr.r_mut(), null_res.null_space_arr.r());
+    //Result
+    println!(
+        "Nullification via QR. Value of |A*B|_2, where B=null(A): {}",
+        res.view_flat().norm_2()
+    );
+}

proc-macro/Cargo.toml

@@ -19,7 +19,7 @@ categories = ["mathematics", "science"]
 proc-macro = true
 
 [dependencies]
-syn = { version = "1.0", features = ["full"]}
+syn = { version = "2.0.100", features = ["full"]}
 quote = "1.0"
 darling = "0.20"
 

src/dense/linalg.rs

@@ -1,14 +1,19 @@
 //! Linear algebra routines
 
-use crate::{MatrixInverse, MatrixPseudoInverse, MatrixQr, MatrixSvd, RlstScalar};
+use crate::{
+    MatrixId, MatrixInverse, MatrixNull, MatrixPseudoInverse, MatrixQr, MatrixSvd, RlstScalar,
+};
 
 use self::lu::MatrixLu;
 
+pub mod interpolative_decomposition;
 pub mod inverse;
 pub mod lu;
+pub mod null_space;
 pub mod pseudo_inverse;
 pub mod qr;
 pub mod svd;
+pub mod triangular_arrays;
 
 /// Return true if stride is column major as required by Lapack.
 pub fn assert_lapack_stride(stride: [usize; 2]) {
@@ -20,10 +25,21 @@ pub fn assert_lapack_stride(stride: [usize; 2]) {
 }
 
 /// Marker trait for objects that support Matrix decompositions.
-pub trait LinAlg: MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse {}
+pub trait LinAlg:
+    MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse + MatrixId + MatrixNull
+{
+}
 
-impl<T: RlstScalar + MatrixInverse + MatrixQr + MatrixSvd + MatrixLu + MatrixPseudoInverse> LinAlg
-    for T
+impl<
+        T: RlstScalar
+            + MatrixInverse
+            + MatrixQr
+            + MatrixSvd
+            + MatrixLu
+            + MatrixPseudoInverse
+            + MatrixId
+            + MatrixNull,
+    > LinAlg for T
 {
 }