SVD::solve quietly produces invalid solution for certain matrices

[RReverser]

Example (playground):

let matrix = nalgebra::matrix![
   1., 0., 2.;
   3., 0., 4.;
   5., 0., 6.
];

let svd = matrix.svd(true, true);
let solution = svd.solve(&matrix, 0.0).unwrap();

println!("solution for self: {:.3}", solution);
/*
Output:

solution for self: 
  ┌                                                                   ┐
  │                 1.000                 0.000                -0.000 │
  │ 18014398509481988.000                 0.000 18014398509481988.000 │
  │                 4.000                 0.000                 5.000 │
  └                                                                   ┘
*/

println!("apply solution: {:.3}", matrix * solution);
/*
Output:

apply solution: 
  ┌                      ┐
  │  9.000  0.000 10.000 │
  │ 19.000  0.000 20.000 │
  │ 29.000  0.000 30.000 │
  └                      ┘
*/

I'm not deeply familiar with matrix decompositions, so I don't know if this is a fundamental limitation of SVD and whether it's possible to fix it to produce a correct result in this case, but at least I would expect the .solve() to return an error in this case instead of quietly returning invalid solution.

[Andlon]

The example matrix you use is clearly singular (it has a zero column). From a numerical perspective, there is no meaningful way of distinguishing between a singular matrix and a badly-conditioned but invertible matrix when working with floating-point numbers. Tiny round-off errors will make numbers that should be zero in exact arithmetic non-zero in floating-point arithmetic, which means that there's no way to determine if a matrix truly is singular or not.

However, SVD does have the means to still solve possibly-singular systems by basically solving a least-squares version of the problem. The way it does so is by truncating sufficiently small singular values (this is a numerical variant of the pseudo-inverse). However, you explicitly disable truncating small singular values by passing 0.0 as the second argument.

You'll see that you get the solution that you expect if you replace your solve call with the following:

let solution = svd.solve(&matrix, 1e-14).unwrap();

How would you know that 1e-14 is an appropriate value to pass here though? Well, you couldn't possibly know that unless you're already intimately familiar with the matrix in question. This suggests to me that the API of SVD is rather flawed.

In my opinion, solve should assume that your linear system is in fact square and invertible (and not catastrophically ill-conditioned). Then there should be a separate method like solve_least_squares which would:

1. work for non-square systems
2. give the approximate minimum-norm solution for inconsistent systems

We could do the second by truncating singular values smaller than $\epsilon \sigma\max$, where $\sigma\max$ is the largest singular value. This does however break down for the zero matrix, so we might need to think carefully about that.

[Andlon]

In fact, maybe we should deprecate solve and provide solve_invertible and solve_least_squares or something like that, to make the intention clearer.

[RReverser]

Hmm. In my original use-case I do need solving overdetermined systems for least squares, and due to #672 SVD is my only option.

But yeah, when testing, I'm finding great difficulty in finding an epsilon that would work for arbitrary inputs. I also wasn't sure what epsilon stands for here, I thought it's the precision of the result (aka "how hard solve should try to find a solution"), and from docs I got an impression that 0.0 would be equivane to f64::EPSILON which seemed like what I want.

In my original issue, found by running proptest on my code, the actual values were non-zero, but they were sufficiently small to cause a problem, and I thought replacing them with zeroes would provide for a better minimal repro.

I still kind of feel that, if this can't be fixed due to nature of SVD, solve() should return an error in this case, since this is very obviously not even a least-squares solution, but I'm not sure if it's possible to detect the issue without multiplying matrices back inside solve() and checking the result.

[RReverser]

We could do the second by truncating singular values smaller than ϵσmax, where σmax is the largest singular value. This does however break down for the zero matrix, so we might need to think carefully about that.

Hmm this does seem to work beautifully for this particular matrix (although it still doesn't produce the expected identity matrix as a solution). Just to be sure, this is what you meant?

let solution = svd.solve(&matrix, std::f64::EPSILON * svd.singular_values.max()).unwrap();

Maybe that should be the implied default when eps == 0.0?

[RReverser]

Unfortunately, even for near-zero values this method is too problematic, e.g.: https://play.rust-lang.org/?version=stable&mode=debug&edition=2021&gist=d4e2d88ede0b5ce43e19ecf931c834b8

I guess my only hope is that someone implements #672 someday.

[Andlon]

SVD least squares produces the minimum norm solution. So obtaining the identity here is not the expected solution, the solution that is returned (with a zero middle column) is what you would expect, because its norm is smaller than the identity matrix, and both matrices are valid solutions to the initial problem.

For your example with the near-zero value, the value you've chosen is so small (~9e-20) that it's effectively noise, and gets treated as such by the truncation inside SVD's solve function.

I would also consider it a "coincidence" that QR actually returns the identity matrix here, because the identity matrix is just one of an infinite number of solutions.

[RReverser]

Sadly, this is the kind of values I do expect in realistic scenarios of systems I'm solving, and they're well off from f64::EPSILON, so I expected I'd be fine...

because the identity matrix is just one of an infinite number of solutions.

Hm, even in case of small but non-zero value?

[Andlon]

Sadly, this is the kind of values I do expect in realistic scenarios of systems I'm solving, and they're well off from f64::EPSILON, so I expected I'd be fine...

Well, if you have quantities than differ in magnitude by more orders of magnitude than digits of accuracy we can't really expect too much.

because the identity matrix is just one of an infinite number of solutions.

Hm, even in case of small but non-zero value?

Sorry, with this statement I was thinking about the exactly zero case. But it is still likely more or less a coincidence that QR returns the identity in this case too, because the value is so small as to be indistinguishable from zero if it's added to or subtracted to any numbers.

In any case, if you really need this very badly scaled linear system to give you a "correct" solution, you might consider introducing a diagonal scaling matrix $D$ and modify your system to

$$ A D Y = A $$

so you'd solve the system

$$ B Y = A $$

where $B = AD$. Then solve $X = D^{-1} Y$. Since the matrix is diagonal you can compute its inverse almost exactly, even with extreme differences in magnitude.

The entries of $D$ might for example be

$$ D_{ii} := \| a_i \|^{-1} $$

where $a_i$ is the $i$-th column of $A$. This would basically balance the scale of the columns of the matrix.

I didn't actually try this out, might be some gotchas I didn't consider.

[RReverser]

FWIW those matrices in my case consist of points on a unit sphere, so yeah, it's very much expected that some coordinates will be close to +-1 while others equal or nearly equal to 0.

Thanks for the ideas, I'll look into that after checking out some other possible solutions.

Should I keep this issue open for tracking improvements you suggested in https://github.com/dimforge/nalgebra/issues/1117#issuecomment-1150792084?

[Andlon]

We could perhaps make a new issue with a more actionable set of suggestions, but yeah, let's keep this issue open for the time being, as you've raised some valid concerns.