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.

[cqc-alec]

Issue resurfaces with this matrix:

    let m = M4C::new(
        Complex::new(-0.4999999999999969, 0.00000000000000009847209127792333),
        Complex::new(-0.5000000000000002, -0.0000000000000005145817176957074),
        Complex::new(-0.5000000000000031, -0.00000000000000011326570795854838),
        Complex::new(0.4999999999999997, 0.00000000000000028317808990753396),
        Complex::new(-0.49999999999999994, 0.0000000000000002324728614724245),
        Complex::new(-0.5000000000000029, 0.00000000000000009174755495786714),
        Complex::new(-0.49999999999999994, 0.00000000000000037035002208515673),
        Complex::new(0.4999999999999969, -0.0000000000000000779181294373365),
        Complex::new(-0.5000000000000031, -0.00000000000000010725752877840523),
        Complex::new(-0.4999999999999997, -0.0000000000000004179869031467993),
        Complex::new(-0.499999999999997, -0.000000000000000027676985470314004),
        Complex::new(0.5000000000000002, 0.0000000000000004129804445544165),
        Complex::new(0.5000000000000001, -0.000000000000000490409530759043),
        Complex::new(0.4999999999999969, -0.00000000000000008015506216445231),
        Complex::new(0.5000000000000001, -0.00000000000000032732792387061645),
        Complex::new(-0.5000000000000029, 0.00000000000000010807420148690568),
    );

This time the issue is with division by very small numbers in compute_2x2_uptrig_svd(). I'll try to fix this.

[cqc-alec]

OK, I think the last sequence of commits fixes the SVD issues -- at least, it makes all tests pass in my downstream project, which was apparently particularly sensitive to these edge cases!

I have replaced the function compute_2x2_uptrig_svd() with a new port from LAPACK.

@sebcrozet Would you be prepared to review this PR?