Implement SafeSVD algorithm

[VictorVanthilt]

When performing a singular value decomposition it can be beneficial to try to use LAPACK_DivideAndConquer for its speed but if it fails, use the more stable but slower LAPACK_QRIteration algorithm.

Thereto I implemented a very minimal SafeSVD algorithm that makes an extra copy of the input data, tries LAPACK_DivideAndConquer on the copy, and if it fails performs LAPACK_QRIteration on the original data.

Feel free to comment on what other magic functions I need to add to make things like Mooncake happy.
I haven't added any tests yet, does anyone know a pathological matrix that will fail with YALAPACK.gesdd! but succeed with YALAPACK.gesvd!?

[kshyatt]

You shouldn't need to add anything for mooncake btw, the algorithm is marked as a non differentiable type

[kshyatt]

BTW can this also support the AMD algorithms ROCSOLVER_DivideAndConquer and ROCSOLVER_QRIteration?

[VictorVanthilt]

any reason why DivideAndConquer and Bisection are left out of:
const ROCSOLVER_SVDAlgorithm = Union{ROCSOLVER_QRIteration, ROCSOLVER_Jacobi}?

[kshyatt]

Oh I got it mixed up, neither gpu library has gesdd yet. So no worries!

It might nevertheless be nice to do this in a way that's easily extensible to non LAPACK backends for the future, in case either does add gesdd

[codecov]

Codecov Report

❌ Patch coverage is 46.66667% with 8 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/implementations/svd.jl	46.66%	8 Missing ⚠️

Files with missing lines	Coverage Δ
src/MatrixAlgebraKit.jl	100.00% <ø> (ø)
src/interface/decompositions.jl	73.68% <ø> (ø)
src/implementations/svd.jl	74.39% <46.66%> (-18.28%)	⬇️

... and 3 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[lkdvos]

Given the remark about the GPU implementations, it might really make sense to just spend the time and make it a bit more generic:

struct FallbackAlgorithm{A <: AbstractAlgorithm, B <: AbstractAlgorithm} <: AbstractAlgorithm
    alg1::A
    alg2::B
end

for f in (:svd_compact, :svd_full)
    f! = Symbol(f, :!)

    # hijack non-mutating function to avoid double-copy
    @eval function $f(A, alg::FallbackAlgorithm)
        Ac = copy_input($f, A)
        out = initialize_output($f!, A, alg.alg1)
        return try
            $f!(Ac, out, alg.alg1)
        catch
            copy!(Ac, A)
            $f!(Ac, out, alg.alg2)
        end
    end

    # implement mutating function to make copy
    @eval function $f!(A, out, alg::FallbackAlgorithm)
        Ac = copy_input($f, A)
        return try
            $f!(Ac, out, alg.alg1)
        catch
            $f!(A, out, alg.alg2)
        end
    end
end

I think our implementations of the other functions might even be generic enough to make this "just work".

This also needs some additional tests, and possibly some select_algorithm overloads or additional constructors to make this a bit more convenient

[lkdvos]

Possibly might find some failure cases here: Reference-LAPACK/lapack#672 (comment)

[Jutho]

While I am typically among the first to consider possible generalisations, in this case I am not sure if this is really the right approach. In fact, I would go further and say that I think this should probably be called something like LAPACK_SafeSVD (but could become the default for BLAS matrices), and would like to write a specialised YALAPACK call for this, that doesn't really throw the error but simply does the LAPACK style info checking, and possibly reuse some of the work space between the gesvd! and gesdd! calls.

[kshyatt]

LAPACK_SafeSVD also works from my PoV. I'm just against giving something a generic name that's really LAPACK only.