update states

rework gauging of InfiniteMPS

Author: lkdvos

src/states/abstractmps.jl

@@ -109,20 +109,20 @@ function isfullrank(A::GenericMPSTensor; side = :both)
 end
 
 """
-    makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg=QRpos())
+    makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg=Defalts.alg_qr())
 
 Make the set of MPS tensors full rank by performing a series of orthogonalizations.
 """
-function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg = QRpos())
+function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg = Defaults.alg_qr())
     while true
         i = findfirst(!isfullrank, A)
         isnothing(i) && break
         if !isfullrank(A[i]; side = :left)
-            L, Q = rightorth!(_transpose_tail(A[i]); alg = alg')
+            L, Q = _right_orth!(_transpose_tail(A[i]); alg)
             A[i] = _transpose_front(Q)
             A[i - 1] = A[i - 1] * L
         else
-            A[i], R = leftorth!(A[i]; alg)
+            A[i], R = _left_orth!(A[i]; alg)
             A[i + 1] = _transpose_front(R * _transpose_tail(A[i + 1]))
         end
     end

src/states/finitemps.jl

@@ -211,15 +211,14 @@ function FiniteMPS(As::Vector{<:GenericMPSTensor}; normalize = false, overwrite
     # vectors anyways, maybe deprecate `overwrite`.
     As = overwrite ? As : copy(As)
     N = length(As)
-    for i in 1:(N - 1)
-        As[i], C = leftorth(As[i]; alg = QRpos())
+    As[1] = MatrixAlgebraKit.copy_input(qr_compact, As[1])
+    local C
+    for i in eachindex(As)
+        As[i], C = qr_compact!(As[i]; positive = true)
         normalize && normalize!(C)
-        As[i + 1] = _transpose_front(C * _transpose_tail(As[i + 1]))
+        i == N || (As[i + 1] = _transpose_front(C * _transpose_tail(As[i + 1])))
     end
 
-    As[end], C = leftorth(As[end]; alg = QRpos())
-    normalize && normalize!(C)
-
     A = eltype(As)
     B = typeof(C)
 
@@ -533,11 +532,11 @@ function Base.:+(ψ₁::MPS, ψ₂::MPS) where {MPS <: FiniteMPS}
     F₁ = isometry(
         storagetype(ψ), (_lastspace(ψ₁.AL[1]) ⊕ _lastspace(ψ₂.AL[1]))', _lastspace(ψ₁.AL[1])'
     )
-    F₂ = leftnull(F₁)
+    F₂ = left_null(F₁)
     @assert _lastspace(F₂) == _lastspace(ψ₂.AL[1])
 
     AL = ψ₁.AL[1] * F₁' + ψ₂.AL[1] * F₂'
-    ψ.ALs[1], R = leftorth!(AL)
+    ψ.ALs[1], R = left_orth!(AL)
 
     for i in 2:halfN
         AL₁ = _transpose_front(F₁ * _transpose_tail(ψ₁.AL[i]))
@@ -546,11 +545,11 @@ function Base.:+(ψ₁::MPS, ψ₂::MPS) where {MPS <: FiniteMPS}
         F₁ = isometry(
             storagetype(ψ), (_lastspace(AL₁) ⊕ _lastspace(ψ₂.AL[i]))', _lastspace(AL₁)'
         )
-        F₂ = leftnull(F₁)
+        F₂ = left_null(F₁)
         @assert _lastspace(F₂) == _lastspace(ψ₂.AL[i])
 
         AL = _transpose_front(R * _transpose_tail(AL₁ * F₁' + AL₂ * F₂'))
-        ψ.ALs[i], R = leftorth!(AL)
+        ψ.ALs[i], R = left_orth!(AL)
     end
 
     C₁ = F₁ * ψ₁.C[halfN]
@@ -560,11 +559,11 @@ function Base.:+(ψ₁::MPS, ψ₂::MPS) where {MPS <: FiniteMPS}
     F₁ = isometry(
         storagetype(ψ), _firstspace(ψ₁.AR[end]) ⊕ _firstspace(ψ₂.AR[end]), _firstspace(ψ₁.AR[end])
     )
-    F₂ = leftnull(F₁)
+    F₂ = left_null(F₁)
     @assert _lastspace(F₂) == _firstspace(ψ₂.AR[end])'
 
     AR = F₁ * _transpose_tail(ψ₁.AR[end]) + F₂ * _transpose_tail(ψ₂.AR[end])
-    L, AR′ = rightorth!(AR)
+    L, AR′ = right_orth!(AR)
     ψ.ARs[end] = _transpose_front(AR′)
 
     for i in Iterators.reverse((halfN + 1):(length(ψ) - 1))
@@ -574,11 +573,11 @@ function Base.:+(ψ₁::MPS, ψ₂::MPS) where {MPS <: FiniteMPS}
         F₁ = isometry(
             storagetype(ψ), _firstspace(ψ₁.AR[i]) ⊕ _firstspace(AR₂), _firstspace(ψ₁.AR[i])
         )
-        F₂ = leftnull(F₁)
+        F₂ = left_null(F₁)
         @assert _lastspace(F₂) == _firstspace(AR₂)'
 
         AR = _transpose_tail(_transpose_front(F₁ * AR₁ + F₂ * AR₂) * L)
-        L, AR′ = rightorth!(AR)
+        L, AR′ = right_orth!(AR)
         ψ.ARs[i] = _transpose_front(AR′)
     end
 

src/states/ortho.jl

@@ -21,7 +21,7 @@ $(TYPEDFIELDS)
     verbosity::Int = VERBOSE_WARN
 
     "algorithm used for orthogonalization of the tensors"
-    alg_orth = QRpos()
+    alg_orth = LAPACK_HouseholderQR(; positive = true)
     "algorithm used for the eigensolver"
     alg_eigsolve = _GAUGE_ALG_EIGSOLVE
     "minimal amount of iterations before using the eigensolver steps"
@@ -46,7 +46,7 @@ $(TYPEDFIELDS)
     verbosity::Int = VERBOSE_WARN
 
     "algorithm used for orthogonalization of the tensors"
-    alg_orth = LQpos()
+    alg_orth = LAPACK_HouseholderLQ(; positive = true)
     "algorithm used for the eigensolver"
     alg_eigsolve = _GAUGE_ALG_EIGSOLVE
     "minimal amount of iterations before using the eigensolver steps"
@@ -73,18 +73,18 @@ end
 
 function MixedCanonical(;
         tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter,
-        verbosity::Int = VERBOSE_WARN, alg_orth = QRpos(),
+        verbosity::Int = VERBOSE_WARN, alg_orth = LAPACK_HouseholderQR(; positive = true),
         alg_eigsolve = _GAUGE_ALG_EIGSOLVE,
         eig_miniter::Int = 10, order::Symbol = :LR
     )
-    if alg_orth isa QR || alg_orth isa QRpos
+    if alg_orth isa LAPACK_HouseholderQR
         alg_leftorth = alg_orth
-        alg_rightorth = alg_orth'
-    elseif alg_orth isa LQ || alg_orth isa LQpos
-        alg_leftorth = alg_orth'
+        alg_rightorth = LAPACK_HouseholderLQ(; alg_orth.kwargs...)
+    elseif alg_orth isa LAPACK_HouseholderLQ
+        alg_leftorth = LAPACK_HouseholderQR(; alg_orth.kwargs...)
         alg_rightorth = alg_orth
     else
-        throw(ArgumentError("Invalid orthogonalization algorithm: $(typeof(alg_orth))"))
+        alg_leftorth = alg_rightorth = alg_orth
     end
 
     left = LeftCanonical(;
@@ -145,45 +145,54 @@ function gaugefix!(ψ::InfiniteMPS, A, C₀, alg::RightCanonical)
 end
 
 @doc """
-    regauge!(AC::GenericMPSTensor, C::MPSBondTensor; alg=QRpos()) -> AL
-    regauge!(CL::MPSBondTensor, AC::GenericMPSTensor; alg=LQpos()) -> AR
+    regauge!(AC::GenericMPSTensor, C::MPSBondTensor; alg) -> AL
+    regauge!(CL::MPSBondTensor, AC::GenericMPSTensor; alg) -> AR
 
 Bring updated `AC` and `C` tensors back into a consistent set of left or right canonical
-tensors. This minimizes `∥AC_i - AL_i * C_i∥` or `∥AC_i - C_{i-1} * AR_i∥`. The optimal algorithm uses
-`Polar()` decompositions, but `QR`-based algorithms are typically more performant.
+tensors. This minimizes `∥AC_i - AL_i * C_i∥` or `∥AC_i - C_{i-1} * AR_i∥`.
+
+The `alg` is passed on to [`left_orth!`](@extref MatrixAlgebraKit) and
+[`right_orth!`](@extref MatrixAlgebraKit), and can be used to control the kind of 
+factorization used. By default, this is set to a (positive) QR/LQ, even though the
+optimal algorithm would use a polar decompositions instead, sacrificing a bit of
+performance for accuracy.
 
 !!! note
     Computing `AL` is slightly faster than `AR`, as it avoids an intermediate transposition.
 """
 regauge!
 
-function regauge!(AC::GenericMPSTensor, C::MPSBondTensor; alg = QRpos())
-    Q_AC, _ = leftorth!(AC; alg)
-    Q_C, _ = leftorth!(C; alg)
+function regauge!(
+        AC::GenericMPSTensor, C::MPSBondTensor; alg = Defaults.alg_qr()
+    )
+    Q_AC, _ = _left_orth!(AC; alg)
+    Q_C, _ = _left_orth!(C; alg)
     return mul!(AC, Q_AC, Q_C')
 end
-function regauge!(AC::Vector{<:GenericMPSTensor}, C::Vector{<:MPSBondTensor}; alg = QRpos())
+function regauge!(AC::Vector{<:GenericMPSTensor}, C::Vector{<:MPSBondTensor}; kwargs...)
     for i in 1:length(AC)
-        regauge!(AC[i], C[i]; alg)
+        regauge!(AC[i], C[i]; kwargs...)
     end
     return AC
 end
-function regauge!(CL::MPSBondTensor, AC::GenericMPSTensor; alg = LQpos())
+function regauge!(
+        CL::MPSBondTensor, AC::GenericMPSTensor; alg = Defaults.alg_lq()
+    )
     AC_tail = _transpose_tail(AC)
-    _, Q_AC = rightorth!(AC_tail; alg)
-    _, Q_C = rightorth!(CL; alg)
+    _, Q_AC = _right_orth!(AC_tail; alg)
+    _, Q_C = _right_orth!(CL; alg)
     AR_tail = mul!(AC_tail, Q_C', Q_AC)
     return repartition!(AC, AR_tail)
 end
-function regauge!(CL::Vector{<:MPSBondTensor}, AC::Vector{<:GenericMPSTensor}; alg = LQpos())
+function regauge!(CL::Vector{<:MPSBondTensor}, AC::Vector{<:GenericMPSTensor}; kwargs...)
     for i in length(CL):-1:1
-        regauge!(CL[i], AC[i]; alg)
+        regauge!(CL[i], AC[i]; kwargs...)
     end
     return CL
 end
 # fix ambiguity + error
-regauge!(::MPSBondTensor, ::MPSBondTensor; alg = QRpos()) = error("method ambiguity")
-function regauge!(::Vector{<:MPSBondTensor}, ::Vector{<:MPSBondTensor}; alg = QRpos())
+regauge!(::MPSBondTensor, ::MPSBondTensor; kwargs...) = error("method ambiguity")
+function regauge!(::Vector{<:MPSBondTensor}, ::Vector{<:MPSBondTensor}; kwargs...)
     return error("method ambiguity")
 end
 
@@ -232,17 +241,18 @@ function gauge_eigsolve_step!(it::IterativeSolver{LeftCanonical}, state)
     if iter ≥ it.eig_miniter
         alg_eigsolve = updatetol(it.alg_eigsolve, 1, ϵ^2)
         _, vec = fixedpoint(flip(TransferMatrix(A, AL)), C[end], :LM, alg_eigsolve)
-        _, C[end] = leftorth!(vec; alg = it.alg_orth)
+        _, C[end] = _left_orth!(vec; alg = it.alg_orth)
     end
     return C[end]
 end
 
 function gauge_orth_step!(it::IterativeSolver{LeftCanonical}, state)
     (; AL, C, A_tail, CA_tail) = state
     for i in 1:length(AL)
+        # repartition!(A_tail[i], AL[i])
         mul!(CA_tail[i], C[i - 1], A_tail[i])
         repartition!(AL[i], CA_tail[i])
-        AL[i], C[i] = leftorth!(AL[i]; alg = it.alg_orth)
+        AL[i], C[i] = _left_orth!(AL[i]; alg = it.alg_orth)
     end
     normalize!(C[end])
     return C[end]
@@ -289,7 +299,7 @@ function gauge_eigsolve_step!(it::IterativeSolver{RightCanonical}, state)
     if iter ≥ it.eig_miniter
         alg_eigsolve = updatetol(it.alg_eigsolve, 1, ϵ^2)
         _, vec = fixedpoint(TransferMatrix(A, AR), C[end], :LM, alg_eigsolve)
-        C[end], _ = rightorth!(vec; alg = it.alg_orth)
+        C[end], _ = _right_orth!(vec; alg = it.alg_orth)
     end
     return C[end]
 end
@@ -299,7 +309,7 @@ function gauge_orth_step!(it::IterativeSolver{RightCanonical}, state)
     for i in length(AR):-1:1
         AC = mul!(AR[i], A[i], C[i])   # use AR as temporary storage for A * C
         tmp = repartition!(AC_tail[i], AC)
-        C[i - 1], tmp = rightorth!(tmp; alg = it.alg_orth)
+        C[i - 1], tmp = _right_orth!(tmp; alg = it.alg_orth)
         repartition!(AR[i], tmp)       # TODO: avoid doing this every iteration
     end
     normalize!(C[end])

src/states/orthoview.jl

@@ -60,7 +60,7 @@ function Base.getindex(v::CView{<:FiniteMPS, E}, i::Int)::E where {E}
         end
 
         for j in Iterators.reverse((i + 1):center)
-            v.parent.Cs[j], tmp = rightorth(_transpose_tail(v.parent.ACs[j]); alg = LQpos())
+            v.parent.Cs[j], tmp = lq_compact!(_transpose_tail(v.parent.ACs[j]); positive = true)
             v.parent.ARs[j] = _transpose_front(tmp)
             if j != i + 1 # last AC not needed
                 v.parent.ACs[j - 1] = _mul_tail(v.parent.ALs[j - 1], v.parent.Cs[j])
@@ -76,7 +76,7 @@ function Base.getindex(v::CView{<:FiniteMPS, E}, i::Int)::E where {E}
         end
 
         for j in center:i
-            v.parent.ALs[j], v.parent.Cs[j + 1] = leftorth(v.parent.ACs[j]; alg = QRpos())
+            v.parent.ALs[j], v.parent.Cs[j + 1] = qr_compact(v.parent.ACs[j]; positive = true)
             if j != i # last AC not needed
                 v.parent.ACs[j + 1] = _mul_front(v.parent.Cs[j + 1], v.parent.ARs[j + 1])
             end
@@ -89,10 +89,10 @@ end
 function Base.setindex!(v::CView{<:FiniteMPS}, vec, i::Int)
     if ismissing(v.parent.Cs[i + 1])
         if !ismissing(v.parent.ALs[i])
-            v.parent.Cs[i + 1], temp = rightorth(_transpose_tail(v.parent.AC[i + 1]); alg = LQpos())
+            v.parent.Cs[i + 1], temp = lq_compact!(_transpose_tail(v.parent.AC[i + 1]); positive = true)
             v.parent.ARs[i + 1] = _transpose_front(temp)
         else
-            v.parent.ALs[i], v.parent.Cs[i + 1] = leftorth(v.parent.AC[i]; alg = QRpos())
+            v.parent.ALs[i], v.parent.Cs[i + 1] = qr_compact(v.parent.AC[i]; positive = true)
         end
     end
 

src/utility/utility.jl

@@ -149,3 +149,31 @@ function check_unambiguous_braiding(V::VectorSpace)
     return check_unambiguous_braiding(Bool, V) ||
         throw(ArgumentError("cannot unambiguously braid $V"))
 end
+
+# temporary workaround for the fact that left_orth and right_orth are poorly designed:
+function _left_orth!(t; alg::MatrixAlgebraKit.AbstractAlgorithm)
+    if alg isa LAPACK_HouseholderQR
+        return left_orth!(t; kind = :qr, alg_qr = alg)
+    elseif alg isa LAPACK_HouseholderLQ
+        return left_orth!(t; kind = :qr, alg_qr = LAPACK_HouseholderQR(; alg.kwargs...))
+    elseif alg isa PolarViaSVD
+        return left_orth!(t; kind = :polar, alg_polar = alg)
+    elseif alg isa LAPACK_SVDAlgorithm
+        return left_orth!(t; kind = :svd, alg_svd = alg)
+    else
+        error(lazy"unkown algorithm $alg")
+    end
+end
+function _right_orth!(t; alg::MatrixAlgebraKit.AbstractAlgorithm)
+    if alg isa LAPACK_HouseholderLQ
+        return right_orth!(t; kind = :lq, alg_lq = alg)
+    elseif alg isa LAPACK_HouseholderQr
+        return right_orth!(t; kind = :lq, alg_lq = LAPACK_HouseholderLQ(; alg.kwargs...))
+    elseif alg isa PolarViaSVD
+        return right_orth!(t; kind = :polar, alg_polar = alg)
+    elseif alg isa LAPACK_SVDAlgorithm
+        return right_orth!(t; kind = :svd, alg_svd = alg)
+    else
+        error(lazy"unkown algorithm $alg")
+    end
+end

test/states.jl

@@ -207,7 +207,7 @@ module TestStates
 
         @test real(e2) ≤ real(e1)
 
-        window, envs = timestep(window, ham, 0.1, 0.0, TDVP2(; trscheme = truncdim(20)), envs)
+        window, envs = timestep(window, ham, 0.1, 0.0, TDVP2(; trscheme = truncrank(20)), envs)
         window, envs = timestep(window, ham, 0.1, 0.0, TDVP(), envs)
 
         e3 = expectation_value(window, (2, 3) => O)