Implement SafeSVD algorithm

src/MatrixAlgebraKit.jl

@@ -33,7 +33,7 @@ export left_orth!, right_orth!, left_null!, right_null!
 
 export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert,
     LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations,
-    LAPACK_DivideAndConquer, LAPACK_Jacobi
+    LAPACK_DivideAndConquer, LAPACK_Jacobi, SafeSVD
 export GLA_HouseholderQR, GLA_QRIteration, GS_QRIteration
 export LQViaTransposedQR
 export PolarViaSVD, PolarNewton

src/implementations/svd.jl

@@ -131,6 +131,16 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm)
         isempty(alg_kwargs) ||
             throw(ArgumentError("invalid keyword arguments for LAPACK_DivideAndConquer"))
         YALAPACK.gesdd!(A, view(S, 1:minmn, 1), U, Vᴴ)
+    elseif alg isa SafeSVD
+        isempty(alg_kwargs) ||
+            throw(ArgumentError("invalid keyword arguments for SafeSVD"))
+        # extra copy to avoid modifying input if sdd goes wrong
+        A′ = copy(A)
+        try
+            YALAPACK.gesdd!(A′, view(S, 1:minmn, 1), U, Vᴴ)
+        catch
+            YALAPACK.gesvd!(A, view(S, 1:minmn, 1), U, Vᴴ)
+        end
     elseif alg isa LAPACK_Bisection
         throw(ArgumentError("LAPACK_Bisection is not supported for full SVD"))
     elseif alg isa LAPACK_Jacobi
@@ -164,6 +174,16 @@ function svd_compact!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm)
         isempty(alg_kwargs) ||
             throw(ArgumentError("invalid keyword arguments for LAPACK_DivideAndConquer"))
         YALAPACK.gesdd!(A, diagview(S), U, Vᴴ)
+    elseif alg isa SafeSVD
+        isempty(alg_kwargs) ||
+            throw(ArgumentError("invalid keyword arguments for SafeSVD"))
+        # extra copy to avoid modifying input if sdd goes wrong
+        A′ = copy(A)
+        try
+            YALAPACK.gesdd!(A′, diagview(S), U, Vᴴ)
+        catch
+            YALAPACK.gesvd!(A, diagview(S), U, Vᴴ)
+        end
     elseif alg isa LAPACK_Bisection
         YALAPACK.gesvdx!(A, diagview(S), U, Vᴴ; alg_kwargs...)
     elseif alg isa LAPACK_Jacobi

src/interface/decompositions.jl

@@ -153,11 +153,24 @@ see also [`gaugefix!`](@ref).
 """
 @algdef LAPACK_Jacobi
 
+"""
+    SafeSVD(; fixgauge::Bool = true)
+
+Algorithm type to denote a combination of LAPACK_DivideAndConquer and LAPACK_QRIteration
+for computing the singular value decomposition of a general matrix.
+This algorithm first tries to perform the SVD using the faster LAPACK_DivideAndConquer method,
+and falls back to the more stable LAPACK_QRIteration method if numerical issues are detected.
+The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the singular vectors,
+see also [`gaugefix!`](@ref), [`LAPACK_DivideAndConquer`](@ref) and [`LAPACK_QRIteration`](@ref).
+"""
+@algdef SafeSVD
+
 const LAPACK_SVDAlgorithm = Union{
     LAPACK_QRIteration,
     LAPACK_Bisection,
     LAPACK_DivideAndConquer,
     LAPACK_Jacobi,
+    SafeSVD,
 }
 
 # =========================