Make the lapack SVD struct public?

[Makogan]

Hey, could you please make the svd trait public so that one can use the lapack crate to make an actual SVD decomposition? Given that currently the methods called SVD in regular nalgebra are not actually computing an SVD decomposition, it would be nice to at least be able to construct the SVD in generic code.

[Finomnis]

"Given that currently the methods called SVD in regular nalgebra are not actually computing an SVD decomposition" - don't you think this is the point that should be addressed, instead of asking for opening possibilities for a workaround?

Also, don't just throw random statements around without backing them up with examples. Of course the svd functions are supposed to calculate an SVD. If they never do that, someone would have already noticed. And I tried, for a random 3x3 matrix it does calculate the proper SVD no problem. So don't just simply throw statements like this into the room.

I suppose you found a corner case where they don't, so please provide said cornercase so that people can actually try to fix it. Don't tell people that their library is bad and needs to be partially replaced with a different library. This is just really disrespectful.

[Makogan]

Yes it works for square matrices but for efficiency purposes it is not computing the full SVD for rectangular matrices, for example:

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);
}

Mathematically the SVD decomposition is defined as computing square matrices for U and V_T which is not what the svd method returns, it only computes the first columns which are necessary for solving the linear system. Applications which need the full SVD with U and V_T being orthogonal matrices cannot use nalgebras implementation.

So as I said, nalgebra is not computing the SVD decomposition, it;s computing submatrices of the SVD decomposition.

[Finomnis]

What exactly are you referring to when you want to make something public? You say struct in the title, but trait in the question.

This library has no SVD trait. It only has an SVD struct, and that one is already public.

[Makogan]

the struct

use nalgebra_lapack::SVD;

SVD::new(&matrix);

Cannot be used in a generic context, because it will demand that the generic scalar type implement the following trait:

T: nalgebra_lapack::svd::SVDScalar<R, R>,

This trait however has been made private and thus no generic code outside of nalgebra can explicitely put that trait restriction on the type. So one cannot write a function that is generic over the matrix type and that also uses the struct to compute the full SVD.

[Finomnis]

I'm unsure what you mean; it should work just like this:

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 │
  └       ┘

Would you mind being a bit more specific why this doesn't work for your usecase and what exactly you would like to change and why?

[Finomnis]

@Makogan In case you just didn't know how to use the library, that's no problem :) I would never blame anyone for that.

I apologize about my harsh tone earlier, I'm afraid I just came across too many XY problems this week and it seems this one was one too many for me. No hard feelings though.

[Makogan]

I'm unsure what you mean; it should work just like this:

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 │
  └       ┘

Would you mind being a bit more specific why this doesn't work for your usecase and what exactly you would like to change and why?

You are right that it works in non generic contexts. I am saying that in generic contexts the compiler will require a private trait. I have tried to make an MVE but rust playground doesn't seem to have nalgebra lapack exposed so it might not work exactly:

use nalgebra as na;
use nalgebra::*;
use na::allocator::Allocator;
use nalgebra_lapack::SVD;

pub fn null_basis<T: RealField, D: Dim, E: Dim>(
    matrix: &OMatrix<T, D, E>,
)
where
    DefaultAllocator: Allocator<T, D> + Allocator<T, E>,
    DefaultAllocator: Allocator<T, D, E> + Allocator<T, D, E>,
    DefaultAllocator: Allocator<T, na::base::dimension::Const<1>, D> + Allocator<T, D, E>,
    DefaultAllocator: Allocator<T, na::base::dimension::Const<1>, E> + Allocator<T, D, E>,
    DefaultAllocator: Allocator<T, <D as DimMin<E>>::Output, E>,
    DefaultAllocator: Allocator<T, <D as DimMin<E>>::Output>,
    DefaultAllocator: Allocator<T, D, <D as DimMin<E>>::Output>,
    DefaultAllocator:
        Allocator<T, <<D as DimMin<E>>::Output as DimSub<Const<1>>>::Output>,
    DefaultAllocator: nalgebra::allocator::Allocator<T, E, E>,
    DefaultAllocator: Allocator<T, <E as DimMin<E>>::Output>,
    T: Scalar + ClosedMul + Copy,
    D: DimMin<E>,
    E: DimMin<E>,
    <D as DimMin<E>>::Output: DimSub<Const<1>>,
{
    let svd = nalgebra_lapack::SVD::new(matrix.clone());
}

fn main() {

}

[Makogan]

I have something similar in my code which returns the following error:

error[E0277]: the trait bound `T: nalgebra_lapack::svd::SVDScalar<E, E>` is not satisfied
   --> crates/algebra/src/lib.rs:254:21
    |
254 | {let svd = SVD::new(matrix.clone());
    |            -------- ^^^^^^^^^^^^^^ the trait `nalgebra_lapack::svd::SVDScalar<E, E>` is not implemented for `T`
    |            |
    |            required by a bound introduced by this call
    |

[Finomnis]

"rust playground doesn't seem to have nalgebra lapack" - Make your own local project then with cargo new. That's what I did to run the previously shown code.

I agree with you, though, that it might be benefitial to make this trait public. I'm not sure if a dev is listening in, though.

Until then, you can always make your own trait:

use na::allocator::Allocator;
use nalgebra as na;
use nalgebra::*;

pub trait OurSVDScalar<R: DimMin<C>, C: Dim>: Scalar
where
    DefaultAllocator: Allocator<Self, R, R>
        + Allocator<Self, R, C>
        + Allocator<Self, DimMinimum<R, C>>
        + Allocator<Self, C, C>,
{
    /// Computes the SVD decomposition of `m`.
    fn compute(m: OMatrix<Self, R, C>) -> Option<nalgebra_lapack::SVD<Self, R, C>>;
}

impl<R: DimMin<C>, C: Dim> OurSVDScalar<R, C> for f32
where
    DefaultAllocator: Allocator<Self, R, R>
        + Allocator<Self, R, C>
        + Allocator<Self, DimMinimum<R, C>>
        + Allocator<Self, C, C>,
{
    fn compute(m: OMatrix<Self, R, C>) -> Option<nalgebra_lapack::SVD<Self, R, C>> {
        nalgebra_lapack::SVD::new(m)
    }
}

impl<R: DimMin<C>, C: Dim> OurSVDScalar<R, C> for f64
where
    DefaultAllocator: Allocator<Self, R, R>
        + Allocator<Self, R, C>
        + Allocator<Self, DimMinimum<R, C>>
        + Allocator<Self, C, C>,
{
    fn compute(m: OMatrix<Self, R, C>) -> Option<nalgebra_lapack::SVD<Self, R, C>> {
        nalgebra_lapack::SVD::new(m)
    }
}

pub fn null_basis<T: OurSVDScalar<R, C>, R: DimMin<C>, C: Dim>(matrix: &OMatrix<T, R, C>)
where
    DefaultAllocator: Allocator<T, R, R>
        + Allocator<T, R, C>
        + Allocator<T, DimMinimum<R, C>>
        + Allocator<T, C, C>,
{
    let svd = <T as OurSVDScalar<R, C>>::compute(matrix.clone());
}

[Makogan]

I will try it out, thank you for the suggestion.