Make gaugefix! optional (#95)

* move gaugefix code to gaugefix section

* add gaugefix keyword argument for enabling/disabling `svd` gauge fixing

* add gaugefix keyword argument for enabling/disabling `eig` gauge fixing

* add gaugefix kewyord argument for enabling/disabling `eigh` gauge fixing

* add gaugefix keyword argument for enabling/disabling `gen_eig` gauge fixing

* update docstrings

* centralize gaugefix function

* remove unnecessary size arguments

* bring eig gaugefix in the same shape

* move include order to ensure things are defined

* update error messages in tests

* update docs

* add gaugefix keyword to GenericLinearAlgebra extension

* rename `dogaugefix` to `do_gauge_fix`

* add global gaugefix toggle

* implement review comments

* more careful about gaugefixing svd_trunc

* rename `lapack_kwargs` to `alg_kwargs`

Author: lkdvos

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
+```
+

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)

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")

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

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

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