Function russell_lab::sv_decomp

pub fn sv_decomp(
    s: &mut Vector, 
    u: &mut Matrix, 
    vt: &mut Matrix, 
    a: &mut Matrix
) -> Result<(), &'static str>

Computes the singular value decomposition (SVD) of a matrix

Finds u, s, and v, such that:

  a  :=  u   ⋅   s   ⋅   vᵀ
(m,n)  (m,m)   (m,n)   (n,n)

Output

• s – min(m,n) vector with the diagonal elements
• u – (m,m) orthogonal matrix
• vt – (n,n) orthogonal matrix with the transpose of v

Input

• a – (m,n) matrix, symmetric or not [will be modified]

Note

1. The matrix a will be modified

Examples

First - 2 x 3 rectangular matrix

// import
use russell_lab::*;

// set matrix
let mut a = Matrix::from(&[
    [3.0, 2.0,  2.0],
    [2.0, 3.0, -2.0],
]);

// allocate output structures
let (m, n) = a.dims();
let min_mn = if m < n { m } else { n };
let mut s = Vector::new(min_mn);
let mut u = Matrix::new(m, m);
let mut vt = Matrix::new(n, n);

// perform SVD
sv_decomp(&mut s, &mut u, &mut vt, &mut a)?;

// define correct data
let s_correct = "┌       ┐\n\
                 │ 5.000 │\n\
                 │ 3.000 │\n\
                 └       ┘";
let u_correct = "┌               ┐\n\
                 │ -0.707 -0.707 │\n\
                 │ -0.707  0.707 │\n\
                 └               ┘";
let vt_correct = "┌                      ┐\n\
                  │ -0.707 -0.707 -0.000 │\n\
                  │ -0.236  0.236 -0.943 │\n\
                  │ -0.667  0.667  0.333 │\n\
                  └                      ┘";

// check solution
assert_eq!(format!("{:.3}", s), s_correct);
assert_eq!(format!("{:.3}", u), u_correct);
assert_eq!(format!("{:.3}", vt), vt_correct);

// check SVD: a == u * s * vt
let mut usv = Matrix::new(m, n);
for i in 0..m {
    for j in 0..n {
        for k in 0..min_mn {
            usv[i][j] += u[i][k] * s[k] * vt[k][j];
        }
    }
}
let usv_correct = "┌                               ┐\n\
                   │  3.000000  2.000000  2.000000 │\n\
                   │  2.000000  3.000000 -2.000000 │\n\
                   └                               ┘";
assert_eq!(format!("{:.6}", usv), usv_correct);

Second - 4 x 2 rectangular matrix

// import
use russell_lab::*;

// set matrix
let mut a = Matrix::from(&[
    [2.0, 4.0],
    [1.0, 3.0],
    [0.0, 0.0],
    [0.0, 0.0],
]);

// allocate output structures
let (m, n) = a.dims();
let min_mn = if m < n { m } else { n };
let mut s = Vector::new(min_mn);
let mut u = Matrix::new(m, m);
let mut vt = Matrix::new(n, n);

// perform SVD
sv_decomp(&mut s, &mut u, &mut vt, &mut a)?;

// define correct data
let s_correct = "┌      ┐\n\
                 │ 5.46 │\n\
                 │ 0.37 │\n\
                 └      ┘";
let u_correct = "┌                         ┐\n\
                 │ -0.82 -0.58  0.00  0.00 │\n\
                 │ -0.58  0.82  0.00  0.00 │\n\
                 │  0.00  0.00  1.00  0.00 │\n\
                 │  0.00  0.00  0.00  1.00 │\n\
                 └                         ┘";
let vt_correct = "┌             ┐\n\
                  │ -0.40 -0.91 │\n\
                  │ -0.91  0.40 │\n\
                  └             ┘";

// check solution
assert_eq!(format!("{:.2}", s), s_correct);
assert_eq!(format!("{:.2}", u), u_correct);
assert_eq!(format!("{:.2}", vt), vt_correct);

// check SVD: a == u * s * vt
let mut usv = Matrix::new(m, n);
for i in 0..m {
    for j in 0..n {
        for k in 0..min_mn {
            usv[i][j] += u[i][k] * s[k] * vt[k][j];
        }
    }
}
let usv_correct = "┌     ┐\n\
                   │ 2 4 │\n\
                   │ 1 3 │\n\
                   │ 0 0 │\n\
                   │ 0 0 │\n\
                   └     ┘";
assert_eq!(format!("{}", usv), usv_correct);