0

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

Makogan
  • 8,208
  • 7
  • 44
  • 112
  • Seems maybe related to this [issue](https://github.com/dimforge/nalgebra/issues/1117) as your `x` too is singular. Though it's been too long since I've done some serious maths for me to judge either way. – cafce25 Aug 19 '23 at 07:14
  • The matrix beign singular should not matter, all matrices accept SVD decompositions. – Makogan Aug 19 '23 at 07:43
  • *"it crashes"* - please provide more information. What happens? A panic? A segfault? Your PC starts burning? If a message is attached, what does the message say? – Finomnis Aug 19 '23 at 08:31
  • It crashes because out of bounds access, the problem is that the SVD decomposition is not returning what it should. The dimensions of the matrices are wrong – Makogan Aug 19 '23 at 09:01
  • 1
    @Makogan Are you sure there's one single definition of an MxN SVD? Because online calculators also disagree with each other; https://keisan.casio.com/exec/system/15076953160460# says that a `[1;0;0]` svd is `[-1;0;0], [1], [-1]`, which is exactly what nalgebra returns. While https://www.emathhelp.net/en/calculators/linear-algebra/svd-calculator/ actually returns what you claim it should be; that said, many of the numbers are meaningless because they fall into a `0` eigenvalue and could be removed entirely. I think that's what nalgebra does - it strips meaningless values. – Finomnis Aug 19 '23 at 09:20

2 Answers2

2

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 │
  └       ┘
Finomnis
  • 18,094
  • 1
  • 20
  • 27
2

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.

Mario S. Mommer
  • 41
  • 3
  • And in many cases you need that basis. SVD is one of the numerical stable ways of obtaining the orthogonal complement for a span of vectors. – Makogan Aug 21 '23 at 21:03
  • 2
    If you only need the basis for the orthogonal complement of a span of vectors you can use cheaper methods than an SVD that are about as stable (eg using Householder transforms). – Mario S. Mommer Aug 22 '23 at 08:47