CODEWITHSUNDEEP
mathrustlinear-algebracomputational-geometrynalgebra
Makogan
Makogan

Reputation: 9632

SVD decomposition of nalgebra behaves weirdly

I am trying to compute the SVD of matrices, as a toy example I used a vector.

I ran my code: https://play.rust-lang.org/?version=stable&mode=debug&edition=2021&gist=00110137fbcbda27c95e9d8791bb97cf

use nalgebra::*;

fn main() {
    let x = dmatrix![
        1., 0., 0.
    ];

    println!("{}", x);

    let svd = x.transpose().svd(true, true);

    // dbg!(&svd);

    let rank = svd.rank(1e-9);

    dbg!(rank);

    let x = svd.v_t.as_ref().unwrap().transpose();
    println!("v\n{}", svd.v_t.as_ref().unwrap().transpose());
    println!("u\n{}", svd.u.as_ref().unwrap().transpose());
    let null: OMatrix<_, _, _> = x.columns(rank, 3 - rank).clone_owned();

    println!("{}", null);
}

And I get this output before it crashes:

  ┌       ┐
  │ 1 0 0 │
  └       ┘


v

  ┌   ┐
  │ 1 │
  └   ┘


u

  ┌       ┐
  │ 1 0 0 │
  └       ┘

Which is nonsense. u should be a 3x3 matrix by definition, and it should be the identity matrix in this case, where are the missing vectors???

Upvotes: 0

Views: 425

Answers (2)

Mario S. Mommer
Mario S. Mommer

Reputation: 51

There are a number of ways of defining the SVD, in one of them you can have U just span the range of the linear function (and V the orthogonal complement to its nullspace). So this code is returning a correct SVD with respect to such a definition.

The advantage of such an abridged formulation is that you don't have to invent a basis for the parts of the spaces that do not play any role for your matrix X. In this case, you would have to complement U with a basis for the orthogonal complement of [1,0,0]^T to get the full U. Also, this basis would be any arbitrary basis of that two-dimensional space.

Upvotes: 3

Finomnis
Finomnis

Reputation: 22808

I think nalgebra::linalg::SVD seems to be a performance-optimized version that deviates from the 'classic' definition of an SVD for MxN matrices.

If you want this behaviour, use nalgebra_lapack::SVD instead:

use nalgebra::dmatrix;
use nalgebra_lapack::SVD;

fn main() {
    let x = dmatrix![1., 0., 0.];

    let svd = SVD::new(x).unwrap();

    println!("U: {}", svd.u);
    println!("Σ: {}", svd.singular_values);
    println!("V_t: {}", svd.vt);
}

U: 
  ┌   ┐
  │ 1 │
  └   ┘

Σ: 
  ┌   ┐
  │ 1 │
  └   ┘

V_t: 
  ┌       ┐
  │ 1 0 0 │
  │ 0 1 0 │
  │ 0 0 1 │
  └       ┘

Upvotes: 3

Related Questions