Skip to content

Commit

Permalink
update left and right_virtualspace (#205)
Browse files Browse the repository at this point in the history
  • Loading branch information
@lkdvos
lkdvos authored Dec 12, 2024
1 parent 19fae27 commit 283a62c
Show file tree
Hide file tree
Showing 17 changed files with 113 additions and 164 deletions.
2 changes: 0 additions & 2 deletions src/MPSKit.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,6 @@ export exact_diagonalization
export TransferMatrix
export transfer_left, transfer_right

@deprecate virtualspace left_virtualspace # there is a possible ambiguity when C isn't square, necessitating specifying left or right virtualspace

# Abstract type defs
abstract type Algorithm end

Expand Down
6 changes: 3 additions & 3 deletions src/algorithms/expval.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,8 +66,8 @@ function expectation_value(ψ::AbstractMPS, (inds, O)::Pair)
# left side
T = storagetype(site_type(ψ))
@plansor Vl[-1 -2; -3] := isomorphism(T,
left_virtualspace(ψ, sites[1] - 1),
left_virtualspace(ψ, sites[1] - 1))[-1; -3] *
left_virtualspace(ψ, sites[1]),
left_virtualspace(ψ, sites[1]))[-1; -3] *
conj(Ut[-2])

# middle
Expand Down Expand Up @@ -104,7 +104,7 @@ function expectation_value(ψ::InfiniteMPS, H::InfiniteMPOHamiltonian,
envs::AbstractMPSEnvironments=environments(ψ, H))
return sum(1:length(ψ)) do i
util = fill_data!(similar(ψ.AL[1], right_virtualspace(H, i)[end]), one)
@plansor GR[-1 -2; -3] := r_LL(ψ, i)[-1; -3] * conj(util[-2])
@plansor GR[-1 -2; -3] := r_LL(ψ, i)[-1; -3] * util[-2]
return contract_mpo_expval(ψ.AL[i], leftenv(envs, i, ψ), H[i][:, 1, 1, end], GR)
end
end
Expand Down
2 changes: 1 addition & 1 deletion src/algorithms/fidelity_susceptibility.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ function fidelity_susceptibility(state::Union{FiniteMPS,InfiniteMPS}, H₀::T,
Tos = LeftGaugedQP(rand, state)
for (i, ac) in enumerate(state.AC)
temp = ∂∂AC(i, state, V, venvs) * ac
help = fill_data!(similar(ac, utilleg(Tos)), one)
help = fill_data!(similar(ac, auxiliaryspace(Tos)), one)
@plansor Tos[i][-1 -2; -3 -4] := temp[-1 -2; -4] * help[-3]
end

Expand Down
2 changes: 1 addition & 1 deletion src/algorithms/timestep/timeevmpo.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -327,7 +327,7 @@ function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T}
Vᵣ = right_virtualspace(H, i)[1:(end - 1)]
P = physicalspace(H, i)

h′ = similar(H[i], Vₗ ⊗ P ← P ⊗ Vᵣ')
h′ = similar(H[i], Vₗ ⊗ P ← P ⊗ Vᵣ)
h′[2:end, 1, 1, 2:end] = WA[i]
h′[2:end, 1, 1, 1] = WB[i]
h′[1, 1, 1, 2:end] = WC[i]
Expand Down
22 changes: 11 additions & 11 deletions src/algorithms/toolbox.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,10 +63,10 @@ domain of each eigenvector. The `tol` and `num_vals` keyword arguments are passe
"""
function transfer_spectrum(above::InfiniteMPS; below=above, tol=Defaults.tol, num_vals=20,
sector=first(sectors(oneunit(left_virtualspace(above, 1)))))
init = randomize!(similar(above.AL[1], left_virtualspace(below, 0),
ℂ[typeof(sector)](sector => 1)' * left_virtualspace(above, 0)))
init = randomize!(similar(above.AL[1], left_virtualspace(below, 1),
ℂ[typeof(sector)](sector => 1)' * left_virtualspace(above, 1)))

transferspace = fuse(left_virtualspace(above, 0) * left_virtualspace(below, 0)')
transferspace = fuse(left_virtualspace(above, 1) * left_virtualspace(below, 1)')
num_vals = min(dim(transferspace, sector), num_vals) # we can ask at most this many values
eigenvals, eigenvecs, convhist = eigsolve(flip(TransferMatrix(above.AL, below.AL)),
init, num_vals, :LM; tol=tol)
Expand Down Expand Up @@ -239,26 +239,26 @@ function periodic_boundary_conditions(mpo::InfiniteMPO{O},

# allocate output
output = Vector{O}(undef, L)
V_wrap = left_virtualspace(mpo, 1)'
V_wrap = left_virtualspace(mpo, 1)
ST = storagetype(O)

util = fill!(similar(mpo[1], oneunit(V_wrap)), one(scalartype(O)))
@plansor cup[-1; -2 -3] := id(ST, V_wrap)[-3; -2] * util[-1]
@plansor cup[-1; -2 -3] := id(ST, V_wrap)[-2; -3] * util[-1]

local F_right
for i in 1:L
V_left = i == 1 ? oneunit(V_wrap) : fuse(V_wrap ⊗ left_virtualspace(mpo, i))
V_right = i == L ? oneunit(V_wrap) : fuse(V_wrap' ⊗ right_virtualspace(mpo, i)')
V_right = i == L ? oneunit(V_wrap) : fuse(V_wrap ⊗ right_virtualspace(mpo, i))
output[i] = similar(mpo[i],
V_left * physicalspace(mpo, i) ←
physicalspace(mpo, i) * V_right)
F_left = i == 1 ? cup : F_right
F_right = i == L ? cup :
isomorphism(ST, V_right ← V_wrap * right_virtualspace(mpo, i)')
@plansor output[i][-1 -2; -3 -4] = F_left[-1; 1 2] *
τ[-3 1; 4 3] *
mpo[i][2 -2; 3 5] *
conj(F_right[-4; 4 5])
isomorphism(ST, V_right ← V_wrap' * right_virtualspace(mpo, i))
@plansor contractcheck = true output[i][-1 -2; -3 -4] = F_left[-1; 1 2] *
τ[-3 1; 4 3] *
mpo[i][2 -2; 3 5] *
conj(F_right[-4; 4 5])
end

mpo isa SparseMPO && dropzeros!.(output) # the above process fills sparse mpos with zeros.
Expand Down
16 changes: 7 additions & 9 deletions src/environments/abstract_envs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,10 @@ Base.unlock(envs::AbstractMPSEnvironments) = unlock(envs.lock);
# Allocating tensors
# ------------------

# TODO: fix the fucking left/right virtualspace bullshit
# TODO: storagetype stuff
function allocate_GL(bra::AbstractMPS, mpo::AbstractMPO, ket::AbstractMPS, i::Int)
T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
V = left_virtualspace(bra, i - 1) ⊗ left_virtualspace(mpo, i)' ←
left_virtualspace(ket, i - 1)
V = left_virtualspace(bra, i) ⊗ left_virtualspace(mpo, i)' ←
left_virtualspace(ket, i)
if V isa BlockTensorKit.TensorMapSumSpace
TT = blocktensormaptype(spacetype(bra), numout(V), numin(V), T)
else
Expand All @@ -30,7 +28,7 @@ end

function allocate_GR(bra::AbstractMPS, mpo::AbstractMPO, ket::AbstractMPS, i::Int)
T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
V = right_virtualspace(ket, i) ⊗ right_virtualspace(mpo, i)' ←
V = right_virtualspace(ket, i) ⊗ right_virtualspace(mpo, i) ←
right_virtualspace(bra, i)
if V isa BlockTensorKit.TensorMapSumSpace
TT = blocktensormaptype(spacetype(bra), numout(V), numin(V), T)
Expand All @@ -42,8 +40,8 @@ end

function allocate_GBL(bra::QP, mpo::AbstractMPO, ket::QP, i::Int)
T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
V = left_virtualspace(bra.left_gs, i - 1) ⊗ left_virtualspace(mpo, i)' ←
auxiliaryspace(ket)' ⊗ left_virtualspace(ket.left_gs, i - 1)
V = left_virtualspace(bra, i) ⊗ left_virtualspace(mpo, i)' ←
auxiliaryspace(ket)' ⊗ left_virtualspace(ket, i)
if V isa BlockTensorKit.TensorMapSumSpace
TT = blocktensormaptype(spacetype(bra), numout(V), numin(V), T)
else
Expand All @@ -54,8 +52,8 @@ end

function allocate_GBR(bra::QP, mpo::AbstractMPO, ket::QP, i::Int)
T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
V = right_virtualspace(ket.right_gs, i) ⊗ right_virtualspace(mpo, i)' ←
auxiliaryspace(ket)' ⊗ right_virtualspace(bra.right_gs, i)
V = right_virtualspace(ket, i) ⊗ right_virtualspace(mpo, i) ←
auxiliaryspace(ket)' ⊗ right_virtualspace(bra, i)
if V isa BlockTensorKit.TensorMapSumSpace
TT = blocktensormaptype(spacetype(bra), numout(V), numin(V), T)
else
Expand Down
18 changes: 3 additions & 15 deletions src/environments/finite_envs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,30 +34,18 @@ function environments(below, operator, above, leftstart, rightstart)
rightenvs)
end

function environments(below::FiniteMPS{S}, O::DenseMPO, above=nothing) where {S}
N = length(below)
leftstart = isomorphism(storagetype(S),
left_virtualspace(below, 0) ⊗ space(O[1], 1)' ←
left_virtualspace(something(above, below), 0))
rightstart = isomorphism(storagetype(S),
right_virtualspace(something(above, below), N) ⊗
space(O[N], 4)' ←
right_virtualspace(below, length(below)))
return environments(below, O, above, leftstart, rightstart)
end

function environments(below::FiniteMPS{S}, O::Union{FiniteMPO,FiniteMPOHamiltonian},
above=nothing) where {S}
Vl_bot = left_virtualspace(below, 0)
Vl_bot = left_virtualspace(below, 1)
Vl_mid = left_virtualspace(O, 1)
Vl_top = isnothing(above) ? left_virtualspace(below, 0) : left_virtualspace(above, 0)
Vl_top = isnothing(above) ? left_virtualspace(below, 1) : left_virtualspace(above, 1)
leftstart = isomorphism(storagetype(S), Vl_bot ⊗ Vl_mid' ← Vl_top)

N = length(below)
Vr_bot = right_virtualspace(below, N)
Vr_mid = right_virtualspace(O, N)
Vr_top = isnothing(above) ? right_virtualspace(below, N) : right_virtualspace(above, N)
rightstart = isomorphism(storagetype(S), Vr_top ⊗ Vr_mid' ← Vr_bot)
rightstart = isomorphism(storagetype(S), Vr_top ⊗ Vr_mid ← Vr_bot)

return environments(below, O, above, leftstart, rightstart)
end
Expand Down
6 changes: 3 additions & 3 deletions src/operators/abstractmpo.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,7 @@ Compute the mpo tensor that arises from multiplying MPOs.
function fuse_mul_mpo(O1::MPOTensor, O2::MPOTensor)
T = promote_type(scalartype(O1), scalartype(O2))
F_left = fuser(T, left_virtualspace(O2), left_virtualspace(O1))
F_right = fuser(T, right_virtualspace(O2)', right_virtualspace(O1)')
F_right = fuser(T, right_virtualspace(O2), right_virtualspace(O1))
@plansor O[-1 -2; -3 -4] := F_left[-1; 1 2] *
O2[1 5; -3 3] *
O1[2 -2; 5 4] *
Expand Down Expand Up @@ -203,7 +203,7 @@ function add_physical_charge(O::BraidingTensor, charge::Sector)
sectortype(O) === typeof(charge) || throw(SectorMismatch())
auxspace = Vect[typeof(charge)](charge => 1)
V = left_virtualspace(O) ⊗ fuse(physicalspace(O), auxspace) ←
fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)'
fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)
return BraidingTensor{scalartype(O)}(V)
end
function add_physical_charge(O::AbstractBlockTensorMap{<:Any,<:Any,2,2}, charge::Sector)
Expand All @@ -213,7 +213,7 @@ function add_physical_charge(O::AbstractBlockTensorMap{<:Any,<:Any,2,2}, charge:

Odst = similar(O,
left_virtualspace(O) ⊗ fuse(physicalspace(O), auxspace) ←
fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)')
fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O))
for (I, v) in nonzero_pairs(O)
Odst[I] = add_physical_charge(v, charge)
end
Expand Down
55 changes: 23 additions & 32 deletions src/operators/mpo.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ const FiniteMPO{O<:MPOTensor} = MPO{O,Vector{O}}

function FiniteMPO(Os::AbstractVector{O}) where {O<:MPOTensor}
for i in eachindex(Os)[1:(end - 1)]
dual(right_virtualspace(Os[i])) == left_virtualspace(Os[i + 1]) ||
right_virtualspace(Os[i]) == left_virtualspace(Os[i + 1]) ||
throw(SpaceMismatch("unmatching virtual spaces at site $i"))
end
return FiniteMPO{O}(Os)
Expand All @@ -38,7 +38,7 @@ const InfiniteMPO{O<:MPOTensor} = MPO{O,PeriodicVector{O}}

function InfiniteMPO(Os::AbstractVector{O}) where {O<:MPOTensor}
for i in eachindex(Os)
dual(right_virtualspace(Os[i])) == left_virtualspace(Os[mod1(i + 1, end)]) ||
right_virtualspace(Os[i]) == left_virtualspace(Os[mod1(i + 1, end)]) ||
throw(SpaceMismatch("umatching virtual spaces at site $i"))
end
return InfiniteMPO{O}(Os)
Expand Down Expand Up @@ -139,8 +139,8 @@ function Base.convert(::Type{TensorMap}, mpo::FiniteMPO)
U_left = ones(scalartype(mpo), V_left)'

V_right = right_virtualspace(mpo, length(mpo))
@assert V_right == oneunit(V_right)'
U_right = ones(scalartype(mpo), V_right')
@assert V_right == oneunit(V_right)
U_right = ones(scalartype(mpo), V_right)

tensors = vcat(U_left, parent(mpo), U_right)
indices = [[i, -i, -(2N - i + 1), i + 1] for i in 1:length(mpo)]
Expand All @@ -166,10 +166,10 @@ function Base.:+(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
A = storagetype(TO)

# left half
F₁ = isometry(A, (right_virtualspace(mpo1, 1) ⊕ right_virtualspace(mpo2, 1))',
right_virtualspace(mpo1, 1)')
F₁ = isometry(A, (right_virtualspace(mpo1, 1) ⊕ right_virtualspace(mpo2, 1)),
right_virtualspace(mpo1, 1))
F₂ = leftnull(F₁)
@assert _lastspace(F₂) == right_virtualspace(mpo2, 1)
@assert _lastspace(F₂) == right_virtualspace(mpo2, 1)'

@plansor O[-3 -1 -2; -4] := mpo1[1][-1 -2; -3 1] * conj(F₁[-4; 1]) +
mpo2[1][-1 -2; -3 1] * conj(F₂[-4; 1])
Expand All @@ -184,10 +184,10 @@ function Base.:+(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
@plansor O₂[-1 -2; -3 -4] := R[-1; 1] * F₂[1; 2] * mpo2[i][2 -2; -3 -4]

# incorporate fusers from right side
F₁ = isometry(A, (right_virtualspace(mpo1, i) ⊕ right_virtualspace(mpo2, i))',
right_virtualspace(mpo1, i)')
F₁ = isometry(A, (right_virtualspace(mpo1, i) ⊕ right_virtualspace(mpo2, i)),
right_virtualspace(mpo1, i))
F₂ = leftnull(F₁)
@assert _lastspace(F₂) == right_virtualspace(mpo2, i)
@assert _lastspace(F₂) == right_virtualspace(mpo2, i)'
@plansor O[-3 -1 -2; -4] := O₁[-1 -2; -3 1] * conj(F₁[-4; 1]) +
O₂[-1 -2; -3 1] * conj(F₂[-4; 1])

Expand Down Expand Up @@ -248,8 +248,8 @@ function Base.:*(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
S = spacetype(TO)
if (left_virtualspace(mpo1, 1) != oneunit(S) ||
left_virtualspace(mpo2, 1) != oneunit(S)) ||
(right_virtualspace(mpo1, N)' != oneunit(S) ||
right_virtualspace(mpo2, N)' != oneunit(S))
(right_virtualspace(mpo1, N) != oneunit(S) ||
right_virtualspace(mpo2, N) != oneunit(S))
@warn "left/right virtual space is not trivial, fusion may not be unique"
# this is a warning because technically any isomorphism that fuses the left/right
# would work and for now I dont feel like figuring out if this is important
Expand All @@ -261,13 +261,8 @@ function Base.:*(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
# note order of mpos: mpo1 * mpo2 * state -> mpo2 on top of mpo1
local Fᵣ # trick to make Fᵣ defined in the loop
for i in 1:N
Fₗ = i != 1 ? Fᵣ :
isomorphism(A, fuse(left_virtualspace(mpo2, i), left_virtualspace(mpo1, i)),
left_virtualspace(mpo2, i) * left_virtualspace(mpo1, i))
Fᵣ = isomorphism(A,
fuse(right_virtualspace(mpo2, i)', right_virtualspace(mpo1, i)'),
right_virtualspace(mpo2, i)' * right_virtualspace(mpo1, i)')

Fₗ = i != 1 ? Fᵣ : fuser(A, left_virtualspace(mpo2, i), left_virtualspace(mpo1, i))
Fᵣ = fuser(A, right_virtualspace(mpo2, i), right_virtualspace(mpo1, i))
@plansor O[i][-1 -2; -3 -4] := Fₗ[-1; 1 4] * mpo2[i][1 2; -3 3] *
mpo1[i][4 -2; 2 5] *
conj(Fᵣ[-4; 3 5])
Expand All @@ -284,11 +279,8 @@ function Base.:*(mpo::FiniteMPO, mps::FiniteMPS)

local Fᵣ # trick to make Fᵣ defined in the loop
for i in 1:length(mps)
Fₗ = i != 1 ? Fᵣ :
isomorphism(TT, fuse(left_virtualspace(A[i]), left_virtualspace(mpo, i)),
left_virtualspace(A[i]) * left_virtualspace(mpo, i))
Fᵣ = isomorphism(TT, fuse(right_virtualspace(A[i])', right_virtualspace(mpo, i)'),
right_virtualspace(A[i])' * right_virtualspace(mpo, i)')
Fₗ = i != 1 ? Fᵣ : fuser(TT, left_virtualspace(mps, i), left_virtualspace(mpo, i))
Fᵣ = fuser(TT, right_virtualspace(mps, i), right_virtualspace(mpo, i))
A[i] = _fuse_mpo_mps(mpo[i], A[i], Fₗ, Fᵣ)
end

Expand All @@ -298,12 +290,11 @@ function Base.:*(mpo::FiniteMPO, mps::FiniteMPS)
end

function Base.:*(mpo::InfiniteMPO, mps::InfiniteMPS)
check_length(mpo, mps)
L = check_length(mpo, mps)
T = promote_type(scalartype(mpo), scalartype(mps))
fusers = PeriodicArray(map(mps.AL, mpo) do al, mp
return fuser(T, _firstspace(al), _firstspace(mp))
end)
As = map(1:length(mps)) do i
fusers = PeriodicArray(fuser.(T, left_virtualspace.(Ref(mps), 1:L),
left_virtualspace.(Ref(mpo), 1:L)))
As = map(1:L) do i
return _fuse_mpo_mps(mpo[i], mps.AL[i], fusers[i], fusers[i + 1])
end
return changebonds(InfiniteMPS(As), SvdCut(; trscheme=notrunc()))
Expand Down Expand Up @@ -346,14 +337,14 @@ function TensorKit.dot(bra::FiniteMPS{T}, mpo::FiniteMPO, ket::FiniteMPS{T}) whe
Nhalf = N ÷ 2
# left half
ρ_left = isomorphism(storagetype(T),
left_virtualspace(bra, 0) ⊗ left_virtualspace(mpo, 1)',
left_virtualspace(ket, 0))
left_virtualspace(bra, 1) ⊗ left_virtualspace(mpo, 1)',
left_virtualspace(ket, 1))
T_left = TransferMatrix(ket.AL[1:Nhalf], mpo[1:Nhalf], bra.AL[1:Nhalf])
ρ_left = ρ_left * T_left

# right half
ρ_right = isomorphism(storagetype(T),
right_virtualspace(ket, N) ⊗ right_virtualspace(mpo, N)',
right_virtualspace(ket, N) ⊗ right_virtualspace(mpo, N),
right_virtualspace(ket, length(ket)))
T_right = TransferMatrix(ket.AR[(Nhalf + 1):end], mpo[(Nhalf + 1):end],
bra.AR[(Nhalf + 1):end])
Expand Down
Loading

0 comments on commit 283a62c

Please sign in to comment.