Module peroxide::structure::matrix

Matrix for Scientific computation

Declare matrix

• You can declare matrix by various ways.
  • R’s way - Default
  • MATLAB’s way
  • Python’s way
  • Other macro

R’s way

• Description: Same as R - matrix(Vec<f64>, Row, Col, Shape)

• Type: matrix(Vec<T>, usize, usize, Shape) where T: std::convert::Into<f64> + Copy

  • Shape: Enum for matrix shape - Row & Col

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix(c!(1,2,3,4), 2, 2, Row);
      a.print();
      //       c[0] c[1]
      // r[0]     1    2
      // r[1]     3    4
  
      let b = matrix(c!(1,2,3,4), 2, 2, Col);
      b.print();
      //       c[0] c[1]
      // r[0]     1    3
      // r[1]     2    4
  }

MATLAB’s way

• Description: Similar to MATLAB (But should use &str)

• Type: ml_matrix(&str)

  use peroxide::fuga::*;
  
  fn main() {
      let a = ml_matrix("1 2; 3 4");
      a.print();
      //       c[0] c[1]
      // r[0]     1    2
      // r[1]     3    4
  }

Python’s way

• Description: Declare matrix as vector of vectors.

• Type: py_matrix(Vec<Vec<T>>) where T: std::convert::Into<f64> + Copy

  use peroxide::fuga::*;
  
  fn main() {
      let a = py_matrix(vec![vec![1, 2], vec![3, 4]]);
      a.print();
      //       c[0] c[1]
      // r[0]     1    2
      // r[1]     3    4
  }

Other macro

• Description: R-like macro to declare matrix

• For R,

  # R
  a = matrix(1:4, nrow = 2, ncol = 2, byrow = T)
  print(a)
  #      [,1] [,2]
  # [1,]    1    2
  # [2,]    3    4
  

• For Peroxide,

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix!(1;4;1, 2, 2, Row);
      a.print();
      //       c[0] c[1]
      // r[0]     1    2
      // r[1]     3    4
  }

Basic Method for Matrix

There are some useful methods for Matrix

• row(&self, index: usize) -> Vec<f64> : Extract specific row as Vec<f64>

• col(&self, index: usize) -> Vec<f64> : Extract specific column as Vec<f64>

• diag(&self) -> Vec<f64>: Extract diagonal components as Vec<f64>

• swap(&self, usize, usize, Shape): Swap two rows or columns (unsafe function)

• subs_col(&mut self, usize, Vec<f64>): Substitute column with Vec<f64>

• subs_row(&mut self, usize, Vec<f64>): Substitute row with Vec<f64>

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let mut a = ml_matrix("1 2; 3 4");
  
      a.row(0).print(); // [1, 2]
      a.col(0).print(); // [1, 3]
      a.diag().print(); // [1, 4]
      unsafe {
          a.swap(0, 1, Row);
      }
      a.print();
      //      c[0] c[1]
      // r[0]    3    4
      // r[1]    1    2
  
      let mut b = ml_matrix("1 2;3 4");
      b.subs_col(0, &c!(5, 6));
      b.subs_row(1, &c!(7, 8));
      b.print();
      //       c[0] c[1]
      // r[0]    5    2
      // r[1]    7    8
  }

Read & Write

In peroxide, we can write matrix to csv

CSV (Not recommended)

• csv feature should be required

• write(&self, file_path: &str): Write matrix to csv

• write_with_header(&self, file_path, header: Vec<&str>): Write with header

  use peroxide::fuga::*;
  
  fn main() {
      let a = ml_matrix("1 2;3 4");
      a.write("example_data/matrix.csv").expect("Can't write file");
  
      let b = ml_matrix("1 2; 3 4; 5 6");
      b.write_with_header("example_data/header.csv", vec!["odd", "even"])
          .expect("Can't write header file");
      println!("Complete!")
  }

Also, you can read matrix from csv.

• Type: read(&str, bool, char) -> Result<Matrix, Box<Error>>

• Description: read(file_path, is_header, delimiter)

  use peroxide::fuga::*;
  
  fn main() {
      // `csv` feature should be required
      let a = Matrix::read("example_data/matrix.csv", false, ',')
          .expect("Can't read matrix.csv file");
      a.print();
      //       c[0] c[1]
      // r[0]     1    2
      // r[1]     3    4
  }

Convert to DataFrame (Recommended)

To write columns or rows, DataFrame and nc feature could be the best choice.

• nc feature should be required - netcdf or libnetcdf are prerequisites.

#[macro_use]
extern crate peroxide;
use peroxide::fuga::*;

fn main() {
    let a = matrix(c!(1,2,3,4,5,6), 3, 2, Col);

    // Construct DataFrame
    let mut df = DataFrame::new(vec![]);
    df.push("x", Series::new(a.col(0)));
    df.push("y", Series::new(a.col(1)));

    // Write nc file (`nc` feature should be required)
    df.write_nc("data.nc").expect("Can't write data.nc");

    // Read nc file (`nc` feature should be required)
    let dg = DataFrame::read_nc("data.nc").expect("Can't read data.nc");
    let x: Vec<f64> = dg["x"].to_vec();
    let y: Vec<f64> = dg["y"].to_vec();

    assert_eq!(a.col(0), x);
    assert_eq!(a.col(1), y);
}

Concatenation

There are two options to concatenate matrices.

• cbind: Concatenate two matrices by column direction.
• rbind: Concatenate two matrices by row direction.

use peroxide::fuga::*;

fn main() -> Result<(), Box<dyn Error>> {
    let a = ml_matrix("1 2;3 4");
    let b = ml_matrix("5 6;7 8");

    cbind(a.clone(), b.clone())?.print();
    //      c[0] c[1] c[2] c[3]
    // r[0]    1    2    5    7
    // r[1]    3    4    6    8

    rbind(a, b)?.print();
    //      c[0] c[1]
    // r[0]    1    2
    // r[1]    3    4
    // r[2]    5    6
    // r[3]    7    8

    Ok(())
}

Matrix operations

• In peroxide, can use basic operations between matrices. I’ll show you by examples.

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix!(1;4;1, 2, 2, Row);
      (a.clone() + 1).print(); // -, *, / are also available
      //      c[0] c[1]
      // r[0]    2    3
      // r[1]    4    5
  
      let b = matrix!(5;8;1, 2, 2, Row);
      (a.clone() + b.clone()).print(); // - is also available
      //      c[0] c[1]
      // r[0]    6    8
      // r[1]   10   12
  
      (a.clone() * b.clone()).print(); // Matrix multiplication
      //      c[0] c[1]
      // r[0]   19   22
      // r[1]   43   50
  }

• clone is too annoying - We can use reference operations!

  use peroxide::fuga::*;
  
  fn main() {
      let a = ml_matrix("1 2;3 4");
      let b = ml_matrix("5 6;7 8");
  
      (&a + 1).print();
      (&a + &b).print();
      (&a - &b).print();
      (&a * &b).print();
  }

Extract & modify components

• In peroxide, matrix data is saved as linear structure.

• But you can use two-dimensional index to extract or modify components.

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let mut a = matrix!(1;4;1, 2, 2, Row);
      a[(0,0)].print();   // 1
      a[(0,0)] = 2f64;    // Modify component
      a.print();
      //       c[0] c[1]
      //  r[0]    2    2
      //  r[1]    3    4
  }

Conversion to Vec<f64>

• Just use row or col method (I already showed at Basic method section).

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix!(1;4;1, 2, 2, Row);
      a.row(0).print(); // [1, 2]
      assert_eq!(a.row(0), vec![1f64, 2f64]);
  }

Useful constructor

• zeros(usize, usize): Construct matrix which elements are all zero

• eye(usize): Identity matrix

• rand(usize, usize): Construct random uniform matrix (from 0 to 1)

• rand_with_rng(usize, usize, &mut Rng): Construct random matrix with user-defined RNG

  use peroxide::fuga::*;
  
  fn main() {
      let a = zeros(2, 2);
      assert_eq!(a, ml_matrix("0 0;0 0"));
  
      let b = eye(2);
      assert_eq!(b, ml_matrix("1 0;0 1"));
  
      let c = rand(2, 2);
      c.print(); // Random 2x2 matrix
  }

Linear Algebra

Transpose

• Caution: Transpose does not consume the original value.

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix!(1;4;1, 2, 2, Row);
      a.transpose().print();
      // Or you can use shorter one
      a.t().print();
      //      c[0] c[1]
      // r[0]    1    3
      // r[1]    2    4
  }

LU Decomposition

• Peroxide uses complete pivoting for LU decomposition - Very stable

• Since there are lots of causes to generate error, you should use Option

• lu returns PQLU

  • PQLU has four field - p, q, l , u
  • p means row permutations
  • q means column permutations
  • l means lower triangular matrix
  • u menas upper triangular matrix

• The structure of PQLU is as follows:

  use peroxide::fuga::*;
  
  #[derive(Debug, Clone)]
  pub struct PQLU {
      pub p: Vec<usize>,
      pub q: Vec<usize>,
      pub l: Matrix,
      pub u: Matrix,
  }

• Example of LU decomposition:

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix(c!(1,2,3,4), 2, 2, Row);
      let pqlu = a.lu();
      let (p,q,l,u) = (pqlu.p, pqlu.q, pqlu.l, pqlu.u);
      assert_eq!(p, vec![1]); // swap 0 & 1 (Row)
      assert_eq!(q, vec![1]); // swap 0 & 1 (Col)
      assert_eq!(l, matrix(c!(1,0,0.5,1),2,2,Row));
      //      c[0] c[1]
      // r[0]    1    0
      // r[1]  0.5    1
      assert_eq!(u, matrix(c!(4,3,0,-0.5),2,2,Row));
      //      c[0] c[1]
      // r[0]    4    3
      // r[1]    0 -0.5
  }

Determinant

• Peroxide uses LU decomposition to obtain determinant ($ \mathcal{O}(n^3) $)

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix!(1;4;1, 2, 2, Row);
      assert_eq!(a.det(), -2f64);
  }

Inverse matrix

• Peroxide uses LU decomposition (via GECP) to obtain inverse matrix.
• It needs two sub functions - inv_l, inv_u
  • For inverse of L, U, I use block partitioning. For example, for lower triangular matrix : $$ \begin{aligned} L &= \begin{pmatrix} L_1 & \mathbf{0} \\ L_2 & L_3 \end{pmatrix} \\ L^{-1} &= \begin{pmatrix} L_1^{-1} & \mathbf{0} \\ -L_3^{-1}L_2 L_1^{-1} & L_3^{-1} \end{pmatrix} \end{aligned} $$

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix!(1;4;1, 2, 2, Row);
      a.inv().print();
      //      c[0] c[1]
      // r[0]   -2    1
      // r[1]  1.5 -0.5
  }

Tips for LU, Det, Inverse

• If you save self.lu() rather than the direct use of self.det() or self.lu() then you can get better performance (via memoization)

use peroxide::fuga::*;

fn main() {
    let a = ml_matrix("1 2;3 4");
    let pqlu = a.lu();  // Memoization of LU
    pqlu.det().print(); // Same as a.det() but do not need an additional LU
    pqlu.inv().print(); // Same as a.inv() but do not need an additional LU
}

QR Decomposition (O3 feature only)

• Use dgeqrf of LAPACK

• Return QR structure.

  pub struct QR {
      pub q: Matrix,
      pub r: Matrix,
  }

• Example

  use peroxide::fuga::*;
  
  let a = ml_matrix("
      1 -1 4;
      1 4 -2;
      1 4 2;
      1 -1 0
  ");
  
  // QR decomposition
  let qr = a.qr();
  
  qr.q().print();
  //         c[0]    c[1]    c[2]
  // r[0]    -0.5     0.5    -0.5
  // r[1]    -0.5 -0.5000  0.5000
  // r[2]    -0.5 -0.5000    -0.5
  // r[3]    -0.5     0.5  0.5000
  
  qr.r().print();
  //      c[0] c[1] c[2]
  // r[0]   -2   -3   -2
  // r[1]    0   -5    2
  // r[2]    0    0   -4

Singular Value Decomposition (O3 feature only)

• Use dgesvd of LAPACK

• Return SVD structure

  use peroxide::fuga::*;
  
  pub struct SVD {
      pub s: Vec<f64>,
      pub u: Matrix,
      pub vt: Matrix,
  }

• Example

  use peroxide::fuga::*;
  
  fn main() {
      let a = ml_matrix("3 2 2;2 3 -2");
      #[cfg(feature="O3")]
      {
          // Full SVD
          let svd = a.svd();
          assert!(eq_vec(&vec![5f64, 3f64], &svd.s, 1e-7));
  
          // Or Truncated SVD
          let svd2 = svd.truncated();
          assert!(eq_vec(&vec![5f64, 3f64], &svd2.s, 1e-7));
      }
      a.print();
  }

Cholesky Decomposition

• Use dpotrf of LAPACK

• Return Matrix (But there can be panic! - Not symmetric or Not positive definite)

• Example

  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = ml_matrix("1 2;2 5");
      #[cfg(feature = "O3")]
      {
          let u = a.cholesky(Upper);
          assert_eq!(u, ml_matrix("1 2;0 1"));
  
          let l = a.cholesky(Lower);
          assert_eq!(l, ml_matrix("1 0;2 1"));
      }
      a.print();
  }

Moore-Penrose Pseudo Inverse

• $ X^\dagger = \left(X^T X\right)^{-1} X^T $

• For O3 feature, peroxide use SVD to obtain pseudo inverse

  #[macro_use]
  extern crate peroxide;
  use peroxide::fuga::*;
  
  fn main() {
      let a = matrix!(1;4;1, 2, 2, Row);
      let pinv_a = a.pseudo_inv();
      let inv_a = a.inv();
  
      assert_eq!(inv_a, pinv_a); // Nearly equal (not actually equal)
  }

Re-exports

• pub use self::Shape::Col;
• pub use self::Shape::Row;

Structs

• DGEQRF
  Temporary data structure from dgeqrf
• DGESVD
• DGETRF
  Temporary data structure from dgetrf
• DPOTRF
• Matrix
  R-like matrix structure

Enums

• LAPACK_STATUS
• POSITIVE_STATUS
• SVD_STATUS
• Shape
  To select matrices’ binding.

Traits

• Error
  Error is a trait representing the basic expectations for error values, i.e., values of type E in Result<T, E>.

Functions

• blas_mul
  Matrix multiplication with BLAS
• combine
  Combine separated matrix to one matrix
• diag
• gemm
  GEMM wrapper for Matrixmultiply
• gemv
  General Matrix-Vector multiplication
• gen_householder
• gevm
  General Vector-Matrix multiplication
• inv_l
  Inverse of Lower matrix
• inv_u
  Inverse of upper triangular matrix
• lapack_dgeqrf
  Peroxide version of dgeqrf
• lapack_dgesvd
• lapack_dgetrf
  Peroxide version of dgetrf
• lapack_dgetri
  Peroxide version of dgetri
• lapack_dgetrs
  Peroxide version of dgetrs
• lapack_dpotrf
• matrix
  R-like matrix constructor
• ml_matrix
  Matlab-like matrix constructor
• py_matrix
  Python-like matrix constructor
• r_matrix
  R-like matrix constructor (Explicit ver.)

Type Aliases

• Perms