Module peroxide::structure::matrix

source · [−]
Expand description

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)

    extern crate peroxide;
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

    extern crate peroxide;
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,

    a = matrix(1:4:1, 2, 2, Row)
print(a)

  • 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 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

    extern crate peroxide;
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)

    extern crate peroxide;
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
}

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.

    extern crate peroxide;
use peroxide::fuga::*;

fn main() {
    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
}

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!

    extern crate peroxide;
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

  • 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)

    extern crate peroxide;
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:

    extern crate peroxide;
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)
extern crate peroxide;
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

    ⓘ
    extern crate peroxide;
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

    extern crate peroxide;
use peroxide::fuga::*;

pub struct SVD {
    pub s: Vec<f64>,
    pub u: Matrix,
    pub vt: Matrix,
}

  • Example

    extern crate peroxide;
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
PQLU
Data structure for Complete Pivoting LU decomposition
QR
SVD
WAZD

Enums

Form
LAPACK_STATUS
POSITIVE_STATUS
SVD_STATUS
Shape
To select matrices’ binding.
SolveKind
UPLO

Traits

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

Functions

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
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.)
solve

Type Definitions

Perms