fix: SVD loss of accuracy with small non-zero singular values

[ABorgna]

Addresses a numerical stability issue in the SVD implementation that caused loss of accuracy with small non-zero singular values. We now use a relative epsilon rather than a fixed one.

Fixes #1172. Adds a now-fixed regression test from the issue.

I ran the benchmarks to see if computing the Frobenius norm has an impact on performance, but it doesn't seem to be significant.

[cqc-alec]

This indeed fixes the problem in many cases! Unfortunately, I've just hit an example for which the problem persists. If we use the following matrix in the test:

    let m = M4C::new(
        Complex::new(0.5163130224597328, 0.2640110414676673),
        Complex::new(-0.10845827476820835, 0.34220148642893244),
        Complex::new(0.14618038920991627, 0.3278663576677311),
        Complex::new(0.4834191243671928, -0.32029192524071315),
        Complex::new(0.25011023587834597, -0.5693169162970136),
        Complex::new(0.3731433798035375, 0.09455746917888716),
        Complex::new(0.34176716325705053, -0.17712228258452342),
        Complex::new(-0.3732958190779979, -0.49731992995651203),
        Complex::new(0.3732958190779976, -0.49731992995651314),
        Complex::new(0.34176716325704987, 0.17712228258452062),
        Complex::new(0.37314337980353635, -0.09455746917888491),
        Complex::new(-0.25011023587834524, -0.5693169162970155),
        Complex::new(0.4834191243671925, 0.32029192524071587),
        Complex::new(-0.14618038920991572, 0.32786635766773),
        Complex::new(0.10845827476820777, 0.3422014864289315),
        Complex::new(0.5163130224597317, -0.26401104146766996),
    );

then we find (m - u * sigma * v_t).norm() is about 0.12.

[ABorgna]

Uhm, I tried that matrix in the test and I see the deviation is just 4.03e-15.

Are you using the latest code in this branch?

[cqc-alec]

Uhm, I tried that matrix in the test and I see the deviation is just 4.03e-15.

Are you using the latest code in this branch?

I was, but on further investigation it looks like there's a difference in behaviour between platforms: the issue manifests on Linux x86 but not on MacOS arm64.

[ABorgna]

It seems that matrix gets the same discrepant results on main. I'll see if I can find the issue.

[ABorgna]

Added some relative epsilon scaling based on LAPACK's DBDSQR, and included your matrix as a test.

Now it works on both aarch64-apple-darwin and x86_64-unknown-linux-gnu.

[cqc-alec]

Great progress, but not quite there yet, as this example shows:

    let m = M4C::new(
        Complex::new(0.5507428877419177, -0.22418277708063508),
        Complex::new(0.3304118701150042, -0.19300867894787538),
        Complex::new(0.38115865977169494, -0.033799854188208585),
        Complex::new(-0.5788888093103887, 0.13588006620948034),
        Complex::new(-0.5788888093104021, -0.13588006620946957),
        Complex::new(-0.3811586597717028, -0.03379985418820041),
        Complex::new(-0.33041187011500606, -0.19300867894788695),
        Complex::new(0.5507428877419225, 0.22418277708065218),
        Complex::new(-0.5507428877419267, 0.22418277708064818),
        Complex::new(-0.3304118701150135, 0.1930086789478659),
        Complex::new(-0.38115865977169644, 0.03379985418822192),
        Complex::new(0.5788888093104035, -0.13588006620947532),
        Complex::new(-0.5788888093103961, -0.1358800662094592),
        Complex::new(-0.38115865977169044, -0.0337998541882125),
        Complex::new(-0.33041187011500334, -0.19300867894787002),
        Complex::new(0.5507428877419116, 0.22418277708065648),
    );

At least on Linux x86 this gives a discrepancy of about 0.015.