diff --git a/docs/src/user_interface/decompositions.md b/docs/src/user_interface/decompositions.md index 1ae1fefe..60959cc4 100644 --- a/docs/src/user_interface/decompositions.md +++ b/docs/src/user_interface/decompositions.md @@ -72,6 +72,11 @@ eigh_trunc eigh_vals ``` +!!! note "Gauge Degrees of Freedom" + The eigenvectors returned by these functions have residual phase degrees of freedom. + By default, MatrixAlgebraKit applies a gauge fixing convention to ensure reproducible results. + See [Gauge choices](@ref sec_gaugefix) for more details. + Alongside these functions, we provide a LAPACK-based implementation for dense arrays, as provided by the following algorithms: ```@autodocs; canonical=false @@ -89,6 +94,11 @@ eig_trunc eig_vals ``` +!!! note "Gauge Degrees of Freedom" + The eigenvectors returned by these functions have residual phase degrees of freedom. + By default, MatrixAlgebraKit applies a gauge fixing convention to ensure reproducible results. + See [Gauge choices](@ref sec_gaugefix) for more details. + Alongside these functions, we provide a LAPACK-based implementation for dense arrays, as provided by the following algorithms: ```@autodocs; canonical=false @@ -137,6 +147,11 @@ svd_vals svd_trunc ``` +!!! note "Gauge Degrees of Freedom" + The singular vectors returned by these functions have residual phase degrees of freedom. + By default, MatrixAlgebraKit applies a gauge fixing convention to ensure reproducible results. + See [Gauge choices](@ref sec_gaugefix) for more details. + MatrixAlgebraKit again ships with LAPACK-based implementations for dense arrays: ```@autodocs; canonical=false @@ -371,3 +386,48 @@ norm(A * N1') < 1e-14 && norm(A * N2') < 1e-14 && # output true ``` + +## [Gauge choices](@id sec_gaugefix) + +Both eigenvalue and singular value decompositions have residual gauge degrees of freedom even when the eigenvalues or singular values are unique. +These arise from the fact that even after normalization, the eigenvectors and singular vectors are only determined up to a phase factor. + +### Phase Ambiguity in Decompositions + +For the eigenvalue decomposition `A * V = V * D`, if `v` is an eigenvector with eigenvalue `λ` and `|v| = 1`, then so is `e^(iθ) * v` for any real phase `θ`. +When `λ` is non-degenerate (i.e., has multiplicity 1), the eigenvector is unique up to this phase. + +Similarly, for the singular value decomposition `A = U * Σ * Vᴴ`, the singular vectors `u` and `v` corresponding to a non-degenerate singular value `σ` are unique only up to a common phase. +We can replace `u → e^(iθ) * u` and `vᴴ → e^(-iθ) * vᴴ` simultaneously. + +### Gauge Fixing Convention + +To remove this phase ambiguity and ensure reproducible results, MatrixAlgebraKit implements a gauge fixing convention by default. +The convention ensures that **the entry with the largest magnitude in each eigenvector or left singular vector is real and positive**. + +For eigenvectors, this means that for each column `v` of `V`, we multiply by `conj(sign(v[i]))` where `i` is the index of the entry with largest absolute value. + +For singular vectors, we apply the phase factor to both `u` and `v` based on the entry with largest magnitude in `u`. +This preserves the decomposition `A = U * Σ * Vᴴ` while fixing the gauge. + +### Disabling Gauge Fixing + +Gauge fixing is enabled by default for all eigenvalue and singular value decompositions. +If you prefer to obtain the raw results from the underlying computational routines without gauge fixing, you can disable it using the `gaugefix` keyword argument: + +```julia +# With gauge fixing (default) +D, V = eigh_full(A) + +# Without gauge fixing +D, V = eigh_full(A; gaugefix = false) +``` + +The same keyword is available for `eig_full`, `eig_trunc`, `svd_full`, `svd_compact`, and `svd_trunc` functions. +Additionally, the default value can also be controlled with a global toggle using [`MatrixAlgebraKit.default_gaugefix`](@ref). + +```@docs; canonical=false +MatrixAlgebraKit.gaugefix! +MatrixAlgebraKit.default_gaugefix +``` + diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index 7b078381..e09fccfa 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -1,7 +1,7 @@ module MatrixAlgebraKitGenericLinearAlgebraExt using MatrixAlgebraKit -using MatrixAlgebraKit: sign_safe, check_input, diagview +using MatrixAlgebraKit: sign_safe, check_input, diagview, gaugefix!, default_gaugefix using GenericLinearAlgebra: svd!, svdvals!, eigen!, eigvals!, Hermitian, qr! using LinearAlgebra: I, Diagonal, lmul! @@ -13,18 +13,26 @@ for f! in (:svd_compact!, :svd_full!, :svd_vals!) @eval MatrixAlgebraKit.initialize_output(::typeof($f!), A::AbstractMatrix, ::GLA_QRIteration) = nothing end -function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, ::GLA_QRIteration) +function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIteration) F = svd!(A) U, S, Vᴴ = F.U, Diagonal(F.S), F.Vt - return MatrixAlgebraKit.gaugefix!(svd_compact!, U, S, Vᴴ, size(A)...) + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + do_gauge_fix && gaugefix!(svd_compact!, U, Vᴴ) + + return U, S, Vᴴ end -function MatrixAlgebraKit.svd_full!(A::AbstractMatrix, USVᴴ, ::GLA_QRIteration) +function MatrixAlgebraKit.svd_full!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIteration) F = svd!(A; full = true) U, Vᴴ = F.U, F.Vt S = MatrixAlgebraKit.zero!(similar(F.S, (size(U, 2), size(Vᴴ, 1)))) diagview(S) .= F.S - return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, size(A)...) + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + do_gauge_fix && gaugefix!(svd_full!, U, Vᴴ) + + return U, S, Vᴴ end function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix, S, ::GLA_QRIteration) diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl index 4a846f85..d3a04741 100644 --- a/src/MatrixAlgebraKit.jl +++ b/src/MatrixAlgebraKit.jl @@ -76,7 +76,6 @@ include("common/safemethods.jl") include("common/view.jl") include("common/regularinv.jl") include("common/matrixproperties.jl") -include("common/gauge.jl") include("yalapack.jl") include("algorithms.jl") @@ -93,6 +92,8 @@ include("interface/schur.jl") include("interface/polar.jl") include("interface/orthnull.jl") +include("common/gauge.jl") # needs to be defined after the functions are + include("implementations/projections.jl") include("implementations/truncation.jl") include("implementations/qr.jl") diff --git a/src/common/defaults.jl b/src/common/defaults.jl index 332807da..9f48c941 100644 --- a/src/common/defaults.jl +++ b/src/common/defaults.jl @@ -41,3 +41,20 @@ default_pullback_rank_atol(A) = eps(norm(A, Inf))^(3 / 4) Default tolerance for deciding to warn if the provided `A` is not hermitian. """ default_hermitian_tol(A) = eps(norm(A, Inf))^(3 / 4) + + +const DEFAULT_GAUGEFIX = Ref(true) + +@doc """ + default_gaugefix() -> current_value + default_gaugefix(new_value::Bool) -> previous_value + +Global toggle for enabling or disabling the default behavior of gauge fixing the output of the eigen- and singular value decompositions. +""" default_gaugefix + +default_gaugefix() = DEFAULT_GAUGEFIX[] +function default_gaugefix(new_value::Bool) + previous_value = DEFAULT_GAUGEFIX[] + DEFAULT_GAUGEFIX[] = new_value + return previous_value +end diff --git a/src/common/gauge.jl b/src/common/gauge.jl index a9bf985b..e855548f 100644 --- a/src/common/gauge.jl +++ b/src/common/gauge.jl @@ -1,8 +1,53 @@ -function gaugefix!(V::AbstractMatrix) +""" + gaugefix!(f_eig, V) + gaugefix!(f_svd, U, Vᴴ) + +Fix the residual gauge degrees of freedom in the eigen or singular vectors, that are +obtained from the decomposition performed by `f_eig` or `f_svd`. +This is achieved by ensuring that the entry with the largest magnitude in `V` or `U` +is real and positive. +""" gaugefix! + + +function gaugefix!(::Union{typeof(eig_full!), typeof(eigh_full!), typeof(gen_eig_full!)}, V::AbstractMatrix) for j in axes(V, 2) v = view(V, :, j) - s = conj(sign(_argmaxabs(v))) - @inbounds v .*= s + s = sign(_argmaxabs(v)) + @inbounds v .*= conj(s) end return V end + +function gaugefix!(::typeof(svd_full!), U, Vᴴ) + m, n = size(U, 2), size(Vᴴ, 1) + for j in 1:max(m, n) + if j <= min(m, n) + u = view(U, :, j) + v = view(Vᴴ, j, :) + s = sign(_argmaxabs(u)) + u .*= conj(s) + v .*= s + elseif j <= m + u = view(U, :, j) + s = sign(_argmaxabs(u)) + u .*= conj(s) + else + v = view(Vᴴ, j, :) + s = sign(_argmaxabs(v)) + v .*= conj(s) + end + end + return (U, Vᴴ) +end + +function gaugefix!(::Union{typeof(svd_compact!), typeof(svd_trunc!)}, U, Vᴴ) + @assert axes(U, 2) == axes(Vᴴ, 1) + for j in axes(U, 2) + u = view(U, :, j) + v = view(Vᴴ, j, :) + s = sign(_argmaxabs(u)) + u .*= conj(s) + v .*= s + end + return (U, Vᴴ) +end diff --git a/src/implementations/eig.jl b/src/implementations/eig.jl index 8532c29a..bdb6981a 100644 --- a/src/implementations/eig.jl +++ b/src/implementations/eig.jl @@ -81,28 +81,37 @@ end function eig_full!(A::AbstractMatrix, DV, alg::LAPACK_EigAlgorithm) check_input(eig_full!, A, DV, alg) D, V = DV + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_Simple - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_Simple (geev) does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_Simple")) YALAPACK.geev!(A, D.diag, V) else # alg isa LAPACK_Expert - YALAPACK.geevx!(A, D.diag, V; alg.kwargs...) + YALAPACK.geevx!(A, D.diag, V; alg_kwargs...) end - # TODO: make this controllable using a `gaugefix` keyword argument - V = gaugefix!(V) + + do_gauge_fix && (V = gaugefix!(eig_full!, V)) + return D, V end function eig_vals!(A::AbstractMatrix, D, alg::LAPACK_EigAlgorithm) check_input(eig_vals!, A, D, alg) V = similar(A, complex(eltype(A)), (size(A, 1), 0)) + + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_Simple - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_Simple (geev) does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_Simple")) YALAPACK.geev!(A, D, V) else # alg isa LAPACK_Expert - YALAPACK.geevx!(A, D, V; alg.kwargs...) + YALAPACK.geevx!(A, D, V; alg_kwargs...) end + return D end @@ -135,23 +144,30 @@ _gpu_geev!(A, D, V) = throw(MethodError(_gpu_geev!, (A, D, V))) function eig_full!(A::AbstractMatrix, DV, alg::GPU_EigAlgorithm) check_input(eig_full!, A, DV, alg) D, V = DV + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa GPU_Simple - isempty(alg.kwargs) || - @warn "GPU_Simple (geev) does not accept any keyword arguments" + isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_Simple" _gpu_geev!(A, D.diag, V) end - # TODO: make this controllable using a `gaugefix` keyword argument - V = gaugefix!(V) + + do_gauge_fix && (V = gaugefix!(eig_full!, V)) + return D, V end function eig_vals!(A::AbstractMatrix, D, alg::GPU_EigAlgorithm) check_input(eig_vals!, A, D, alg) V = similar(A, complex(eltype(A)), (size(A, 1), 0)) + + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa GPU_Simple - isempty(alg.kwargs) || - @warn "GPU_Simple (geev) does not accept any keyword arguments" + isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_Simple" _gpu_geev!(A, D, V) end + return D end diff --git a/src/implementations/eigh.jl b/src/implementations/eigh.jl index 13d0b9d3..d4a67dd6 100644 --- a/src/implementations/eigh.jl +++ b/src/implementations/eigh.jl @@ -91,32 +91,41 @@ function eigh_full!(A::AbstractMatrix, DV, alg::LAPACK_EighAlgorithm) check_input(eigh_full!, A, DV, alg) D, V = DV Dd = D.diag + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_MultipleRelativelyRobustRepresentations - YALAPACK.heevr!(A, Dd, V; alg.kwargs...) + YALAPACK.heevr!(A, Dd, V; alg_kwargs...) elseif alg isa LAPACK_DivideAndConquer - YALAPACK.heevd!(A, Dd, V; alg.kwargs...) + YALAPACK.heevd!(A, Dd, V; alg_kwargs...) elseif alg isa LAPACK_Simple - YALAPACK.heev!(A, Dd, V; alg.kwargs...) + YALAPACK.heev!(A, Dd, V; alg_kwargs...) else # alg isa LAPACK_Expert - YALAPACK.heevx!(A, Dd, V; alg.kwargs...) + YALAPACK.heevx!(A, Dd, V; alg_kwargs...) end - # TODO: make this controllable using a `gaugefix` keyword argument - V = gaugefix!(V) + + do_gauge_fix && (V = gaugefix!(eigh_full!, V)) + return D, V end function eigh_vals!(A::AbstractMatrix, D, alg::LAPACK_EighAlgorithm) check_input(eigh_vals!, A, D, alg) V = similar(A, (size(A, 1), 0)) + + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_MultipleRelativelyRobustRepresentations - YALAPACK.heevr!(A, D, V; alg.kwargs...) + YALAPACK.heevr!(A, D, V; alg_kwargs...) elseif alg isa LAPACK_DivideAndConquer - YALAPACK.heevd!(A, D, V; alg.kwargs...) + YALAPACK.heevd!(A, D, V; alg_kwargs...) elseif alg isa LAPACK_QRIteration # == LAPACK_Simple - YALAPACK.heev!(A, D, V; alg.kwargs...) + YALAPACK.heev!(A, D, V; alg_kwargs...) else # alg isa LAPACK_Bisection == LAPACK_Expert - YALAPACK.heevx!(A, D, V; alg.kwargs...) + YALAPACK.heevx!(A, D, V; alg_kwargs...) end + return D end @@ -158,35 +167,44 @@ function eigh_full!(A::AbstractMatrix, DV, alg::GPU_EighAlgorithm) check_input(eigh_full!, A, DV, alg) D, V = DV Dd = D.diag + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa GPU_Jacobi - _gpu_heevj!(A, Dd, V; alg.kwargs...) + _gpu_heevj!(A, Dd, V; alg_kwargs...) elseif alg isa GPU_DivideAndConquer - _gpu_heevd!(A, Dd, V; alg.kwargs...) + _gpu_heevd!(A, Dd, V; alg_kwargs...) elseif alg isa GPU_QRIteration # alg isa GPU_QRIteration == GPU_Simple - _gpu_heev!(A, Dd, V; alg.kwargs...) + _gpu_heev!(A, Dd, V; alg_kwargs...) elseif alg isa GPU_Bisection # alg isa GPU_Bisection == GPU_Expert - _gpu_heevx!(A, Dd, V; alg.kwargs...) + _gpu_heevx!(A, Dd, V; alg_kwargs...) else throw(ArgumentError("Unsupported eigh algorithm")) end - # TODO: make this controllable using a `gaugefix` keyword argument - V = gaugefix!(V) + + do_gauge_fix && (V = gaugefix!(eigh_full!, V)) + return D, V end function eigh_vals!(A::AbstractMatrix, D, alg::GPU_EighAlgorithm) check_input(eigh_vals!, A, D, alg) V = similar(A, (size(A, 1), 0)) + + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa GPU_Jacobi - _gpu_heevj!(A, D, V; alg.kwargs...) + _gpu_heevj!(A, D, V; alg_kwargs...) elseif alg isa GPU_DivideAndConquer - _gpu_heevd!(A, D, V; alg.kwargs...) + _gpu_heevd!(A, D, V; alg_kwargs...) elseif alg isa GPU_QRIteration - _gpu_heev!(A, D, V; alg.kwargs...) + _gpu_heev!(A, D, V; alg_kwargs...) elseif alg isa GPU_Bisection - _gpu_heevx!(A, D, V; alg.kwargs...) + _gpu_heevx!(A, D, V; alg_kwargs...) else throw(ArgumentError("Unsupported eigh algorithm")) end + return D end diff --git a/src/implementations/gen_eig.jl b/src/implementations/gen_eig.jl index 07bf6e0c..94dfe47e 100644 --- a/src/implementations/gen_eig.jl +++ b/src/implementations/gen_eig.jl @@ -57,27 +57,36 @@ end function gen_eig_full!(A::AbstractMatrix, B::AbstractMatrix, WV, alg::LAPACK_EigAlgorithm) check_input(gen_eig_full!, A, B, WV, alg) W, V = WV + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_Simple - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_Simple (ggev) does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_Simple")) YALAPACK.ggev!(A, B, W.diag, V, similar(W.diag, eltype(A))) else # alg isa LAPACK_Expert throw(ArgumentError("LAPACK_Expert is not supported for ggev")) end - # TODO: make this controllable using a `gaugefix` keyword argument - V = gaugefix!(V) + + do_gauge_fix && (V = gaugefix!(gen_eig_full!, V)) + return W, V end function gen_eig_vals!(A::AbstractMatrix, B::AbstractMatrix, W, alg::LAPACK_EigAlgorithm) check_input(gen_eig_vals!, A, B, W, alg) V = similar(A, complex(eltype(A)), (size(A, 1), 0)) + + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_Simple - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_Simple (ggev) does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_Simple")) YALAPACK.ggev!(A, B, W, V, similar(W, eltype(A))) else # alg isa LAPACK_Expert throw(ArgumentError("LAPACK_Expert is not supported for ggev")) end + return W end diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index ba54f233..2066b2cd 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -1,4 +1,4 @@ -# Inputs +# Input # ------ copy_input(::typeof(svd_full), A::AbstractMatrix) = copy!(similar(A, float(eltype(A))), A) copy_input(::typeof(svd_compact), A) = copy_input(svd_full, A) @@ -105,51 +105,6 @@ function initialize_output(::typeof(svd_vals!), A::Diagonal, ::DiagonalAlgorithm return eltype(A) <: Real ? diagview(A) : similar(A, real(eltype(A)), size(A, 1)) end -function gaugefix!(::typeof(svd_full!), U, S, Vᴴ, m::Int, n::Int) - for j in 1:max(m, n) - if j <= min(m, n) - u = view(U, :, j) - v = view(Vᴴ, j, :) - s = conj(sign(_argmaxabs(u))) - u .*= s - v .*= conj(s) - elseif j <= m - u = view(U, :, j) - s = conj(sign(_argmaxabs(u))) - u .*= s - else - v = view(Vᴴ, j, :) - s = conj(sign(_argmaxabs(v))) - v .*= s - end - end - return (U, S, Vᴴ) -end - -# Gauge fixation -# -------------- -function gaugefix!(::typeof(svd_compact!), U, S, Vᴴ, m::Int, n::Int) - for j in 1:size(U, 2) - u = view(U, :, j) - v = view(Vᴴ, j, :) - s = conj(sign(_argmaxabs(u))) - u .*= s - v .*= conj(s) - end - return (U, S, Vᴴ) -end - -function gaugefix!(::typeof(svd_trunc!), U, S, Vᴴ, m::Int, n::Int) - for j in 1:min(m, n) - u = view(U, :, j) - v = view(Vᴴ, j, :) - s = conj(sign(_argmaxabs(u))) - u .*= s - v .*= conj(s) - end - return (U, S, Vᴴ) -end - # Implementation # -------------- function svd_full!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm) @@ -164,13 +119,17 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm) one!(Vᴴ) return USVᴴ end + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_QRIteration - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_QRIteration does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_QRIteration")) YALAPACK.gesvd!(A, view(S, 1:minmn, 1), U, Vᴴ) elseif alg isa LAPACK_DivideAndConquer - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_DivideAndConquer does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_DivideAndConquer")) YALAPACK.gesdd!(A, view(S, 1:minmn, 1), U, Vᴴ) elseif alg isa LAPACK_Bisection throw(ArgumentError("LAPACK_Bisection is not supported for full SVD")) @@ -179,60 +138,71 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm) else throw(ArgumentError("Unsupported SVD algorithm")) end + for i in 2:minmn S[i, i] = S[i, 1] S[i, 1] = zero(eltype(S)) end - # TODO: make this controllable using a `gaugefix` keyword argument - gaugefix!(svd_full!, U, S, Vᴴ, m, n) + + do_gauge_fix && gaugefix!(svd_full!, U, Vᴴ) + return USVᴴ end function svd_compact!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm) check_input(svd_compact!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_QRIteration - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_QRIteration does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_QRIteration")) YALAPACK.gesvd!(A, S.diag, U, Vᴴ) elseif alg isa LAPACK_DivideAndConquer - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_DivideAndConquer does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_DivideAndConquer")) YALAPACK.gesdd!(A, S.diag, U, Vᴴ) elseif alg isa LAPACK_Bisection - YALAPACK.gesvdx!(A, S.diag, U, Vᴴ; alg.kwargs...) + YALAPACK.gesvdx!(A, S.diag, U, Vᴴ; alg_kwargs...) elseif alg isa LAPACK_Jacobi - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_Jacobi does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_Jacobi")) YALAPACK.gesvj!(A, S.diag, U, Vᴴ) else throw(ArgumentError("Unsupported SVD algorithm")) end - # TODO: make this controllable using a `gaugefix` keyword argument - gaugefix!(svd_compact!, U, S, Vᴴ, size(A)...) + + do_gauge_fix && gaugefix!(svd_compact!, U, Vᴴ) + return USVᴴ end function svd_vals!(A::AbstractMatrix, S, alg::LAPACK_SVDAlgorithm) check_input(svd_vals!, A, S, alg) U, Vᴴ = similar(A, (0, 0)), similar(A, (0, 0)) + + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa LAPACK_QRIteration - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_QRIteration does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_QRIteration")) YALAPACK.gesvd!(A, S, U, Vᴴ) elseif alg isa LAPACK_DivideAndConquer - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_DivideAndConquer does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_DivideAndConquer")) YALAPACK.gesdd!(A, S, U, Vᴴ) elseif alg isa LAPACK_Bisection - YALAPACK.gesvdx!(A, S, U, Vᴴ; alg.kwargs...) + YALAPACK.gesvdx!(A, S, U, Vᴴ; alg_kwargs...) elseif alg isa LAPACK_Jacobi - isempty(alg.kwargs) || - throw(ArgumentError("LAPACK_Jacobi does not accept any keyword arguments")) + isempty(alg_kwargs) || + throw(ArgumentError("invalid keyword arguments for LAPACK_Jacobi")) YALAPACK.gesvj!(A, S, U, Vᴴ) else throw(ArgumentError("Unsupported SVD algorithm")) end + return S end @@ -365,21 +335,25 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) one!(Vᴴ) return USVᴴ end + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa GPU_QRIteration - isempty(alg.kwargs) || - @warn "GPU_QRIteration does not accept any keyword arguments" + isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_QRIteration" _gpu_gesvd_maybe_transpose!(A, view(S, 1:minmn, 1), U, Vᴴ) elseif alg isa GPU_SVDPolar - _gpu_Xgesvdp!(A, view(S, 1:minmn, 1), U, Vᴴ; alg.kwargs...) + _gpu_Xgesvdp!(A, view(S, 1:minmn, 1), U, Vᴴ; alg_kwargs...) elseif alg isa GPU_Jacobi - _gpu_gesvdj!(A, view(S, 1:minmn, 1), U, Vᴴ; alg.kwargs...) + _gpu_gesvdj!(A, view(S, 1:minmn, 1), U, Vᴴ; alg_kwargs...) else throw(ArgumentError("Unsupported SVD algorithm")) end diagview(S) .= view(S, 1:minmn, 1) view(S, 2:minmn, 1) .= zero(eltype(S)) - # TODO: make this controllable using a `gaugefix` keyword argument - gaugefix!(svd_full!, U, S, Vᴴ, m, n) + + do_gauge_fix && gaugefix!(svd_full!, U, Vᴴ) + return USVᴴ end @@ -387,32 +361,39 @@ function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Ran check_input(svd_trunc!, A, USVᴴ, alg.alg) U, S, Vᴴ = USVᴴ _gpu_Xgesvdr!(A, S.diag, U, Vᴴ; alg.alg.kwargs...) - # TODO: make this controllable using a `gaugefix` keyword argument - gaugefix!(svd_trunc!, U, S, Vᴴ, size(A)...) + # TODO: make sure that truncation is based on maxrank, otherwise this might be wrong - USVᴴtrunc, ind = truncate(svd_trunc!, (U, S, Vᴴ), alg.trunc) - Strunc = diagview(USVᴴtrunc[2]) + (Utr, Str, Vᴴtr), _ = truncate(svd_trunc!, (U, S, Vᴴ), alg.trunc) + # normal `truncation_error!` does not work here since `S` is not the full singular value spectrum - ϵ = sqrt(norm(A)^2 - norm(Strunc)^2) # is there a more accurate way to do this? - return USVᴴtrunc..., ϵ + ϵ = sqrt(norm(A)^2 - norm(diagview(Str))^2) # is there a more accurate way to do this? + + do_gauge_fix = get(alg.alg.kwargs, :gaugefix, default_gaugefix())::Bool + do_gauge_fix && gaugefix!(svd_trunc!, Utr, Vᴴtr) + + return Utr, Str, Vᴴtr, ϵ end function svd_compact!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) check_input(svd_compact!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ + + do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa GPU_QRIteration - isempty(alg.kwargs) || - @warn "GPU_QRIteration does not accept any keyword arguments" + isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_QRIteration" _gpu_gesvd_maybe_transpose!(A, S.diag, U, Vᴴ) elseif alg isa GPU_SVDPolar - _gpu_Xgesvdp!(A, S.diag, U, Vᴴ; alg.kwargs...) + _gpu_Xgesvdp!(A, S.diag, U, Vᴴ; alg_kwargs...) elseif alg isa GPU_Jacobi - _gpu_gesvdj!(A, S.diag, U, Vᴴ; alg.kwargs...) + _gpu_gesvdj!(A, S.diag, U, Vᴴ; alg_kwargs...) else throw(ArgumentError("Unsupported SVD algorithm")) end - # TODO: make this controllable using a `gaugefix` keyword argument - gaugefix!(svd_compact!, U, S, Vᴴ, size(A)...) + + do_gauge_fix && gaugefix!(svd_compact!, U, Vᴴ) + return USVᴴ end _argmaxabs(x) = reduce(_largest, x; init = zero(eltype(x))) @@ -421,16 +402,19 @@ _largest(x, y) = abs(x) < abs(y) ? y : x function svd_vals!(A::AbstractMatrix, S, alg::GPU_SVDAlgorithm) check_input(svd_vals!, A, S, alg) U, Vᴴ = similar(A, (0, 0)), similar(A, (0, 0)) + + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + if alg isa GPU_QRIteration - isempty(alg.kwargs) || - @warn "GPU_QRIteration does not accept any keyword arguments" + isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_QRIteration" _gpu_gesvd_maybe_transpose!(A, S, U, Vᴴ) elseif alg isa GPU_SVDPolar - _gpu_Xgesvdp!(A, S, U, Vᴴ; alg.kwargs...) + _gpu_Xgesvdp!(A, S, U, Vᴴ; alg_kwargs...) elseif alg isa GPU_Jacobi - _gpu_gesvdj!(A, S, U, Vᴴ; alg.kwargs...) + _gpu_gesvdj!(A, S, U, Vᴴ; alg_kwargs...) else throw(ArgumentError("Unsupported SVD algorithm")) end + return S end diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index 1bdf1534..df5d55e0 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -49,18 +49,22 @@ of `R` are non-negative. # General Eigenvalue Decomposition # ------------------------------- """ - LAPACK_Simple() + LAPACK_Simple(; gaugefix::Bool = true) Algorithm type to denote the simple LAPACK driver for computing the Schur or non-Hermitian eigenvalue decomposition of a matrix. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, +see also [`gaugefix!`](@ref). """ @algdef LAPACK_Simple """ - LAPACK_Expert() + LAPACK_Expert(; gaugefix::Bool = true) Algorithm type to denote the expert LAPACK driver for computing the Schur or non-Hermitian eigenvalue decomposition of a matrix. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, +see also [`gaugefix!`](@ref). """ @algdef LAPACK_Expert @@ -77,37 +81,45 @@ eigenvalue decomposition of a non-Hermitian matrix. # Hermitian Eigenvalue Decomposition # ---------------------------------- """ - LAPACK_QRIteration() + LAPACK_QRIteration(; gaugefix::Bool = true) Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the QR Iteration algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_QRIteration """ - LAPACK_Bisection() + LAPACK_Bisection(; gaugefix::Bool = true) Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Bisection algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_Bisection """ - LAPACK_DivideAndConquer() + LAPACK_DivideAndConquer(; gaugefix::Bool = true) Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Divide and Conquer algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_DivideAndConquer """ - LAPACK_MultipleRelativelyRobustRepresentations() + LAPACK_MultipleRelativelyRobustRepresentations(; gaugefix::Bool = true) Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix using the Multiple Relatively Robust Representations algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, +see also [`gaugefix!`](@ref). """ @algdef LAPACK_MultipleRelativelyRobustRepresentations @@ -119,21 +131,25 @@ const LAPACK_EighAlgorithm = Union{ } """ - GLA_QRIteration() + GLA_QRIteration(; gaugefix::Bool = true) Algorithm type to denote the GenericLinearAlgebra.jl implementation for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef GLA_QRIteration # Singular Value Decomposition # ---------------------------- """ - LAPACK_Jacobi() + LAPACK_Jacobi(; gaugefix::Bool = true) Algorithm type to denote the LAPACK driver for computing the singular value decomposition of a general matrix using the Jacobi algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the singular vectors, +see also [`gaugefix!`](@ref). """ @algdef LAPACK_Jacobi @@ -199,34 +215,40 @@ end CUSOLVER_HouseholderQR(; positive = false) Algorithm type to denote the standard CUSOLVER algorithm for computing the QR decomposition of -a matrix using Householder reflectors. The keyword `positive=true` can be used to ensure that +a matrix using Householder reflectors. The keyword `positive = true` can be used to ensure that the diagonal elements of `R` are non-negative. """ @algdef CUSOLVER_HouseholderQR """ - CUSOLVER_QRIteration() + CUSOLVER_QRIteration(; gaugefix::Bool = true) Algorithm type to denote the CUSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the QR Iteration algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_QRIteration """ - CUSOLVER_SVDPolar() + CUSOLVER_SVDPolar(; gaugefix::Bool = true) Algorithm type to denote the CUSOLVER driver for computing the singular value decomposition of a general matrix by using Halley's iterative algorithm to compute the polar decompositon, followed by the hermitian eigenvalue decomposition of the positive definite factor. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the singular +vectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_SVDPolar """ - CUSOLVER_Jacobi() + CUSOLVER_Jacobi(; gaugefix::Bool = true) Algorithm type to denote the CUSOLVER driver for computing the singular value decomposition of a general matrix using the Jacobi algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the singular +vectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_Jacobi @@ -248,21 +270,25 @@ for more information. does_truncate(::TruncatedAlgorithm{<:CUSOLVER_Randomized}) = true """ - CUSOLVER_Simple() + CUSOLVER_Simple(; gaugefix::Bool = true) Algorithm type to denote the simple CUSOLVER driver for computing the non-Hermitian eigenvalue decomposition of a matrix. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, +see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_Simple const CUSOLVER_EigAlgorithm = Union{CUSOLVER_Simple} """ - CUSOLVER_DivideAndConquer() + CUSOLVER_DivideAndConquer(; gaugefix::Bool = true) Algorithm type to denote the CUSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Divide and Conquer algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_DivideAndConquer @@ -283,37 +309,45 @@ the diagonal elements of `R` are non-negative. @algdef ROCSOLVER_HouseholderQR """ - ROCSOLVER_QRIteration() + ROCSOLVER_QRIteration(; gaugefix::Bool = true) Algorithm type to denote the ROCSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the QR Iteration algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef ROCSOLVER_QRIteration """ - ROCSOLVER_Jacobi() + ROCSOLVER_Jacobi(; gaugefix::Bool = true) Algorithm type to denote the ROCSOLVER driver for computing the singular value decomposition of a general matrix using the Jacobi algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the singular +vectors, see also [`gaugefix!`](@ref). """ @algdef ROCSOLVER_Jacobi """ - ROCSOLVER_Bisection() + ROCSOLVER_Bisection(; gaugefix::Bool = true) Algorithm type to denote the ROCSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Bisection algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef ROCSOLVER_Bisection """ - ROCSOLVER_DivideAndConquer() + ROCSOLVER_DivideAndConquer(; gaugefix::Bool = true) Algorithm type to denote the ROCSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Divide and Conquer algorithm. +The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +singular vectors, see also [`gaugefix!`](@ref). """ @algdef ROCSOLVER_DivideAndConquer diff --git a/test/amd/svd.jl b/test/amd/svd.jl index 4f58f925..f8a48045 100644 --- a/test/amd/svd.jl +++ b/test/amd/svd.jl @@ -41,7 +41,7 @@ include(joinpath("..", "utilities.jl")) @test ROCArray(diagview(S)) ≈ Sd # ROCArray is necessary because norm of ROCArray view with non-unit step is broken if alg isa ROCSOLVER_QRIteration - @test_warn "GPU_QRIteration does not accept any keyword arguments" svd_compact!(copy!(Ac, A), (U, S, Vᴴ), ROCSOLVER_QRIteration(; bad = "bad")) + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_compact!(copy!(Ac, A), (U, S, Vᴴ), ROCSOLVER_QRIteration(; bad = "bad")) end end end @@ -83,8 +83,8 @@ end @test ROCArray(diagview(S)) ≈ Sc # ROCArray is necessary because norm of ROCArray view with non-unit step is broken if alg isa ROCSOLVER_QRIteration - @test_warn "GPU_QRIteration does not accept any keyword arguments" svd_full!(copy!(Ac, A), (U, S, Vᴴ), ROCSOLVER_QRIteration(; bad = "bad")) - @test_warn "GPU_QRIteration does not accept any keyword arguments" svd_vals!(copy!(Ac, A), Sc, ROCSOLVER_QRIteration(; bad = "bad")) + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_full!(copy!(Ac, A), (U, S, Vᴴ), ROCSOLVER_QRIteration(; bad = "bad")) + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_vals!(copy!(Ac, A), Sc, ROCSOLVER_QRIteration(; bad = "bad")) end end end diff --git a/test/cuda/svd.jl b/test/cuda/svd.jl index bfe993f8..e2173d03 100644 --- a/test/cuda/svd.jl +++ b/test/cuda/svd.jl @@ -41,7 +41,7 @@ include(joinpath("..", "utilities.jl")) @test CuArray(diagview(S)) ≈ Sd # CuArray is necessary because norm of CuArray view with non-unit step is broken if alg isa CUSOLVER_QRIteration - @test_warn "GPU_QRIteration does not accept any keyword arguments" svd_compact!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_compact!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) end end end @@ -84,8 +84,8 @@ end @test CuArray(diagview(S)) ≈ Sc # CuArray is necessary because norm of CuArray view with non-unit step is broken if alg isa CUSOLVER_QRIteration - @test_warn "GPU_QRIteration does not accept any keyword arguments" svd_full!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) - @test_warn "GPU_QRIteration does not accept any keyword arguments" svd_vals!(copy!(Ac, A), Sc, CUSOLVER_QRIteration(; bad = "bad")) + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_full!(copy!(Ac, A), (U, S, Vᴴ), CUSOLVER_QRIteration(; bad = "bad")) + @test_warn "invalid keyword arguments for GPU_QRIteration" svd_vals!(copy!(Ac, A), Sc, CUSOLVER_QRIteration(; bad = "bad")) end end end diff --git a/test/gen_eig.jl b/test/gen_eig.jl index eebfeee5..8f905b77 100644 --- a/test/gen_eig.jl +++ b/test/gen_eig.jl @@ -42,9 +42,9 @@ using LinearAlgebra: Diagonal A = randn(rng, T, m, m) B = randn(rng, T, m, m) @test_throws ArgumentError("LAPACK_Expert is not supported for ggev") gen_eig_full(A, B; alg = LAPACK_Expert()) - @test_throws ArgumentError("LAPACK_Simple (ggev) does not accept any keyword arguments") gen_eig_full(A, B; alg = LAPACK_Simple(bad = "sad")) + @test_throws ArgumentError("invalid keyword arguments for LAPACK_Simple") gen_eig_full(A, B; alg = LAPACK_Simple(bad = "sad")) @test_throws ArgumentError("LAPACK_Expert is not supported for ggev") gen_eig_vals(A, B; alg = LAPACK_Expert()) - @test_throws ArgumentError("LAPACK_Simple (ggev) does not accept any keyword arguments") gen_eig_vals(A, B; alg = LAPACK_Simple(bad = "sad")) + @test_throws ArgumentError("invalid keyword arguments for LAPACK_Simple") gen_eig_vals(A, B; alg = LAPACK_Simple(bad = "sad")) # a tuple of the input types is passed to `default_algorithm` @test_throws MethodError MatrixAlgebraKit.default_algorithm(gen_eig_full, A, B)