From cdcf46cacb3c6b46b38e7bb678af8cc81db42ab4 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 11 Mar 2022 10:08:16 -0500 Subject: [PATCH 01/19] adding nullspace --- src/ITensorTDVP.jl | 2 + src/nullspace.jl | 131 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 133 insertions(+) create mode 100644 src/nullspace.jl diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 765fc75..7b4a5b3 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -5,8 +5,10 @@ using KrylovKit: exponentiate, eigsolve using Printf using ITensors: AbstractMPS, @debug_check, @timeit_debug, check_hascommoninds, orthocenter +using ITensors.NDTensors include("tdvp.jl") +include("nullspace.jl") export tdvp diff --git a/src/nullspace.jl b/src/nullspace.jl new file mode 100644 index 0000000..b7bbd0b --- /dev/null +++ b/src/nullspace.jl @@ -0,0 +1,131 @@ + +ITensors.itensor(A::ITensor) = A + +# Insert missing diagonal blocks +function insert_diag_blocks!(T::Tensor) + for b in eachdiagblock(T) + blockT = blockview(T, b) + if isnothing(blockT) + # Block was not found in the list, insert it + insertblock!(T, b) + end + end +end +insert_diag_blocks!(T::ITensor) = insert_diag_blocks!(tensor(T)) + +# Reshape into an order-2 ITensor +matricize(T::ITensor, inds::Index...) = matricize(T, inds) + +function matricize(T::ITensor, inds) + left_inds = commoninds(T, inds) + right_inds = uniqueinds(T, inds) + return matricize(T, left_inds, right_inds) +end + +function matricize(T::ITensor, left_inds, right_inds) + CL = combiner(left_inds; dir=ITensors.Out, tags="CL") + CR = combiner(right_inds; dir=ITensors.In, tags="CR") + M = (T * CL) * CR + return M, CL, CR +end + +function setdims(t::NTuple{N,Pair{QN,Int}}, dims::NTuple{N,Int}) where {N} + return first.(t) .=> dims +end + +# XXX: generalize this function +function _getindex(T::DenseTensor{ElT,N}, I1::Colon, I2::UnitRange{Int64}) where {ElT,N} + A = array(T)[I1, I2] + return tensor(Dense(vec(A)), setdims(inds(T), size(A))) +end + +function getblock(i::Index, n::Integer) + return ITensors.space(i)[n] +end + +# Make `Pair{QN,Int}` act like a regular `dim` +NDTensors.dim(qnv::Pair{QN,Int}) = last(qnv) + +Base.:*(qnv::Pair{QN,Int}, d::ITensors.Arrow) = qn(qnv) * d => dim(qnv) + +function getblock_preserve_qns(T::Tensor, b::Block) + # TODO: make `T[b]` preserve QNs + Tb = T[b] + indsTb = getblock.(inds(T), Tuple(b)) .* dir.(inds(T)) + return ITensors.setinds(Tb, indsTb) +end + +function blocksparsetensor(blocks::Dict{B,TB}) where {B,TB} + b1, Tb1 = first(pairs(blocks)) + N = length(b1) + indstypes = typeof.(inds(Tb1)) + blocktype = eltype(Tb1) + indsT = getindex.(indstypes) + # Determine the indices from the blocks + for (b, Tb) in pairs(blocks) + indsTb = inds(Tb) + for n in 1:N + bn = b[n] + indsTn = indsT[n] + if bn > length(indsTn) + resize!(indsTn, bn) + end + indsTn[bn] = indsTb[n] + end + end + T = BlockSparseTensor(blocktype, indsT) + for (b, Tb) in pairs(blocks) + if !isempty(Tb) + T[b] = Tb + end + end + return T +end + +function _nullspace_hermitian(M::Tensor; atol::Real=0.0) + tol = atol + # Insert any missing diagonal blocks + insert_diag_blocks!(M) + #D, U = eigen(Hermitian(M)) + Dᵢₜ, Uᵢₜ = eigen(itensor(M); ishermitian=true) + D = tensor(Dᵢₜ) + U = tensor(Uᵢₜ) + nullspace_blocks = Dict() + for bU in nzblocks(U) + bM = Block(bU[1], bU[1]) + bD = Block(bU[2], bU[2]) + # Assume sorted from largest to smallest + indstart = sum(d -> abs(d) .> tol, storage(D[bD])) + 1 + Ub = getblock_preserve_qns(U, bU) + indstop = lastindex(Ub, 2) + Nb = _getindex(Ub, :, indstart:indstop) + nullspace_blocks[bU] = Nb + end + return blocksparsetensor(nullspace_blocks) +end + +function LinearAlgebra.nullspace(M::Hermitian{<:Number,<:Tensor}; kwargs...) + return _nullspace_hermitian(parent(M); kwargs...) +end + +function LinearAlgebra.nullspace(::Order{2}, M::ITensor, left_inds, right_inds; kwargs...) + @assert order(M) == 2 + M² = prime(dag(M), right_inds) * M + M² = permute(M², right_inds'..., right_inds...) + M²ₜ = tensor(M²) + Nₜ = nullspace(Hermitian(M²ₜ); kwargs...) + indsN = (Index(ind(Nₜ, 1); dir=ITensors.In), Index(ind(Nₜ, 2); dir=ITensors.In)) + N = dag(itensor(ITensors.setinds(Nₜ, indsN))) + # Make the index match the input index + Ñ = replaceinds(N, (ind(N, 1),) => right_inds) + return Ñ +end + +function LinearAlgebra.nullspace(T::ITensor, is...; kwargs...) + M, CL, CR = matricize(T, is...) + @assert order(M) == 2 + cL = commoninds(M, CL) + cR = commoninds(M, CR) + N₂ = nullspace(Order(2), M, cL, cR; kwargs...) + return N₂ * CR +end From 6fab5a3b0112e3c8c07e3c57143c9de6ada242f7 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 11 Mar 2022 16:52:29 -0500 Subject: [PATCH 02/19] Initial commit for subspace expansion --- src/subspace_expansion.jl | 101 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 src/subspace_expansion.jl diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl new file mode 100644 index 0000000..95cf8d4 --- /dev/null +++ b/src/subspace_expansion.jl @@ -0,0 +1,101 @@ + +function replaceind_indval(IV::Tuple, iĩ::Pair) + i, ĩ = iĩ + return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV)) +end + +# atol controls the tolerance cutoff for determining which eigenvectors are in the null +# space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep +# the eigenvectors corresponding to eigenvalues of approximately 1. + +###the most general implementation may be passing MPS, lims, and bondtensor (defaulting to the Id-matrix) +###then we can tensorsum both left and right tensor, and just apply the bondtensor to restore input gauge +###ideally we would also be able to apply tensorsum to the bond tensor instead of that loop below +function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, cutoff, atol=1e-2, kwargs... + ) + N = length(psi) + if !isortho(psi) || orthocenter(psi) != 1 + orthogonalize!(psi, 1) + end + PH.nsite=2 + position!(PH, psi, 1) + for (b, ha) in sweepnext(N; ncenter=2) + orthogonalize!(psi,b) + position!(PH, psi, b) + + if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) + b1 = (ha == 1 ? b + 1 : b) + Δ = (ha == 1 ? +1 : -1) + _=subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),(b,b+Δ);maxdim, cutoff, atol=1e-2, kwargs... + ) + end + end +end + +function subspace_expansion!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum},lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs... + ) + ##this should only work for the case where rlim-llim > 1 + ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) + llim,rlim = lims + n1, n2 = b + @assert n1 + 1 == n2 + PH.nsite=2 + position!(PH,psi,n1) + linkind=commoninds(ψ[n1],ψ[n2]) + NL=nullspace(ψ[n1],linkind) + NR=nullspace(ψ[n2],linkind) + ϕ=ψ[n1]*ψ[n2] + ψ[n1]=prime(ψ[n1],linkind) + if isnothing(bondtensor) + bondtensor=diagITensor(1.0,prime(linkind,1),linkind) + end + + newL,S,newR=_subspace_expand_core(ϕ,PH,NL,NR,;maxdim, cutoff, atol=1e-2, kwargs...) + nullbond=diagITensor(0.0,uniqueinds(newL, ψ[n1]),uniqueinds(newR, ψ[n2])) + + + ALⁿ¹, newl = ITensors.directsum( + ψ[n1], dag(newL), uniqueinds(ψ[n1], newL), uniqueinds(newL, ψ[n1]); tags=("Left",) + ) + ARⁿ², newr = ITensors.directsum( + ψ.AR[n2], dag(newR), uniqueinds(ψ[n2], newR), uniqueinds(newR, ψ[n2]); tags=("Right",) + ) + C,newlr=ITensors.directsum(bondtensor => (prime(linkind,1), linkind), nullbond => (uniqueinds(newL, ψ[n1]), uniqueinds(newR, ψ[n2])),tags=("Left","Right")) + + if rlim==n2 + ψ[n2]=ARⁿ² + elseif rlim>n2 + ψ[n2]=noprime(ARⁿ²*C) + end + + if llim==n1 + ψ[n1]=ALⁿ¹ + elseif llim Date: Fri, 11 Mar 2022 18:10:09 -0500 Subject: [PATCH 03/19] Setup minimal example for debugging subspace expansion --- examples/tdvp.jl | 12 +++++----- src/ITensorTDVP.jl | 7 ++++-- src/subspace_expansion.jl | 50 +++++++++++++++++++-------------------- src/tdvp.jl | 16 +++++++++---- 4 files changed, 48 insertions(+), 37 deletions(-) diff --git a/examples/tdvp.jl b/examples/tdvp.jl index 6b7bc26..b89b347 100644 --- a/examples/tdvp.jl +++ b/examples/tdvp.jl @@ -1,8 +1,7 @@ using ITensors using ITensorTDVP - n = 10 -s = siteinds("S=1/2", n) +s = siteinds("S=1/2", n,conserve_qns=true) function heisenberg(n) os = OpSum() @@ -15,16 +14,17 @@ function heisenberg(n) end H = MPO(heisenberg(n), s) -ψ = randomMPS(s, "↑"; linkdims=10) +ψ = randomMPS(s,i -> isodd(i) ? "Up" : "Dn"; linkdims=10) @show inner(ψ', H, ψ) / inner(ψ, ψ) ϕ = tdvp( H, ψ, - -1; - nsweeps=10, - reverse_step=false, + 0.05im; + nsweeps=1, + nsite=1, + reverse_step=true, normalize=true, maxdim=20, cutoff=1e-10, diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 7b4a5b3..9b80157 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -3,13 +3,16 @@ module ITensorTDVP using ITensors using KrylovKit: exponentiate, eigsolve using Printf +using LinearAlgebra using ITensors: AbstractMPS, @debug_check, @timeit_debug, check_hascommoninds, orthocenter using ITensors.NDTensors +using ITensors.NDTensors: eachdiagblock, blockview + -include("tdvp.jl") include("nullspace.jl") +include("subspace_expansion.jl") +include("tdvp.jl") export tdvp - end diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 95cf8d4..3fb081e 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -11,39 +11,45 @@ end ###the most general implementation may be passing MPS, lims, and bondtensor (defaulting to the Id-matrix) ###then we can tensorsum both left and right tensor, and just apply the bondtensor to restore input gauge ###ideally we would also be able to apply tensorsum to the bond tensor instead of that loop below -function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, cutoff, atol=1e-2, kwargs... - ) - N = length(psi) - if !isortho(psi) || orthocenter(psi) != 1 - orthogonalize!(psi, 1) +""" +expand subspace (2-site like) in a sweep +""" +function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, cutoff, atol=1e-2, kwargs...) + N = length(ψ) + if !isortho(ψ) || orthocenter(ψ) != 1 + orthogonalize!(ψ, 1) end PH.nsite=2 - position!(PH, psi, 1) + nsite=2 + position!(PH, ψ, 1) for (b, ha) in sweepnext(N; ncenter=2) - orthogonalize!(psi,b) - position!(PH, psi, b) + println(b) + orthogonalize!(ψ,b) + position!(PH, ψ, b) if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) b1 = (ha == 1 ? b + 1 : b) Δ = (ha == 1 ? +1 : -1) - _=subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),(b,b+Δ);maxdim, cutoff, atol=1e-2, kwargs... + println("") + _=subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),(b,b+Δ);maxdim, cutoff, atol=atol, kwargs... ) end end + return nothing end -function subspace_expansion!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum},lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs... - ) +function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) ##this should only work for the case where rlim-llim > 1 ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) llim,rlim = lims n1, n2 = b @assert n1 + 1 == n2 PH.nsite=2 - position!(PH,psi,n1) - linkind=commoninds(ψ[n1],ψ[n2]) - NL=nullspace(ψ[n1],linkind) - NR=nullspace(ψ[n2],linkind) + position!(PH,ψ,n1) + linkind=commonind(ψ[n1],ψ[n2]) + @show ψ[n1] + NL=nullspace(ψ[n1],linkind;atol=1e-6) + NR=nullspace(ψ[n2],linkind;atol=1e-6) ϕ=ψ[n1]*ψ[n2] ψ[n1]=prime(ψ[n1],linkind) if isnothing(bondtensor) @@ -77,22 +83,16 @@ function subspace_expansion!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum},lims::Tuple{I return C end - - -function _subspace_expand_core(centerwf::Vector{ITensor}, env,NL,NR;maxdim, cutoff, atol=1e-2, kwargs... - ) +function _subspace_expand_core(centerwf::Vector{ITensor}, env,NL,NR;maxdim, cutoff, atol=1e-2, kwargs...) ϕ = ITensor(1.0) for atensor in centerwf ϕ *= atensor end - - return _subspace_expand_core(ϕ, env, NL, NR;maxdim, cutoff, atol=1e-2, kwargs...) - + return _subspace_expand_core(ϕ, env, NL, NR;maxdim, cutoff, atol=1e-2, kwargs...) end -function _subspace_expand_core(ϕ::ITensor, env,NL,NR;maxdim, cutoff, atol=1e-2, kwargs... - ) - ϕH = env * ϕ #add noprime? +function _subspace_expand_core(ϕ::ITensor, env,NL,NR;maxdim, cutoff, atol=1e-2, kwargs...) + ϕH = env*ϕ #add noprime? ϕH = NL * ϕH * NR U,S,V=svd(ϕH,commoninds(ϕH,NL);maxdim=maxdim, cutoff=cutoff, kwargs...) NL *= dag(U) diff --git a/src/tdvp.jl b/src/tdvp.jl index 6da1573..d097692 100644 --- a/src/tdvp.jl +++ b/src/tdvp.jl @@ -19,7 +19,8 @@ function tdvp(solver, PH, psi0::MPS, t::Number, sweeps; kwargs...)::MPS which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing) svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer") obs = get(kwargs, :observer, NoObserver()) - + atol = get(kwargs, :atol,1e-6) + maxdim=get(kwargs,:maxdim,100) write_when_maxdim_exceeds::Union{Int,Nothing} = get( kwargs, :write_when_maxdim_exceeds, nothing ) @@ -31,10 +32,15 @@ function tdvp(solver, PH, psi0::MPS, t::Number, sweeps; kwargs...)::MPS orthogonalize!(psi, 1) end @assert isortho(psi) && orthocenter(psi) == 1 + #PH.nsite=1 + #position!(PH, psi, 1) + #psid=PH * psi - position!(PH, psi, 1) - + #println("did a successful env itensor mult") for sw in 1:nsweep(sweeps) + if nsite==1 + subspace_expansion_sweep!(psi,PH;maxdim, cutoff, atol=atol, kwargs...) + end sw_time = @elapsed begin maxtruncerr = 0.0 @@ -57,7 +63,8 @@ function tdvp(solver, PH, psi0::MPS, t::Number, sweeps; kwargs...)::MPS elseif nsite == 2 phi1 = psi[b] * psi[b + 1] end - + phiprime=PH*phi1 + println("did a successful env itensor mult") phi1, info = solver(PH, t / 2, phi1) ## if info.converged == 0 @@ -248,6 +255,7 @@ function tdvp(PH, psi0::MPS, t::Number; reverse_step=true, kwargs...) verbosity=get(kwargs, :exponentiate_verbosity, 0), ) function solver(PH, t, psi0) + #psidot=PH*psi0 psi, info = exponentiate(PH, t, psi0; solver_kwargs...) return psi, info end From 478cdb38c7e53437beffe12aa6e8681216e30e69 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 25 Mar 2022 17:15:12 -0400 Subject: [PATCH 04/19] Partial fix to subspace expansion. Works in sweep in single tdvp iteration, breaks for nsteps>1. Print debugging output. --- examples/tdvp.jl | 35 +++++++++++--- src/subspace_expansion.jl | 98 ++++++++++++++++++++++++++++++--------- src/tdvp.jl | 9 ++-- 3 files changed, 111 insertions(+), 31 deletions(-) diff --git a/examples/tdvp.jl b/examples/tdvp.jl index b89b347..897e91b 100644 --- a/examples/tdvp.jl +++ b/examples/tdvp.jl @@ -16,27 +16,50 @@ end H = MPO(heisenberg(n), s) ψ = randomMPS(s,i -> isodd(i) ? "Up" : "Dn"; linkdims=10) -@show inner(ψ', H, ψ) / inner(ψ, ψ) +E0= inner(ψ', H, ψ) / inner(ψ, ψ) ϕ = tdvp( H, ψ, - 0.05im; + 0.2im; nsweeps=1, nsite=1, reverse_step=true, normalize=true, - maxdim=20, + maxdim=50, cutoff=1e-10, outputlevel=1, ) - -@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ) +ϕ = tdvp( + H, + ϕ, + 0.2im; + nsweeps=1, + nsite=1, + reverse_step=true, + normalize=true, + maxdim=50, + cutoff=1e-10, + outputlevel=1, +) +ϕ = tdvp( + H, + ϕ, + 0.2im; + nsweeps=2, + nsite=1, + reverse_step=true, + normalize=true, + maxdim=100, + cutoff=1e-10, + outputlevel=1, +) +@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ) , E0 e2, ϕ2 = dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10) @show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2) -ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1) +ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10) @show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3) diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 3fb081e..1547219 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -30,8 +30,8 @@ function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) b1 = (ha == 1 ? b + 1 : b) Δ = (ha == 1 ? +1 : -1) - println("") - _=subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),(b,b+Δ);maxdim, cutoff, atol=atol, kwargs... + inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) + _=subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... ) end end @@ -43,31 +43,72 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) llim,rlim = lims n1, n2 = b - @assert n1 + 1 == n2 + #@assert n1 + 1 == n2 || n1 -1 == n2 PH.nsite=2 - position!(PH,ψ,n1) - linkind=commonind(ψ[n1],ψ[n2]) - @show ψ[n1] - NL=nullspace(ψ[n1],linkind;atol=1e-6) - NR=nullspace(ψ[n2],linkind;atol=1e-6) - ϕ=ψ[n1]*ψ[n2] - ψ[n1]=prime(ψ[n1],linkind) - if isnothing(bondtensor) - bondtensor=diagITensor(1.0,prime(linkind,1),linkind) + if llim == n1 + @assert rlim == n2+1 + U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=0., kwargs...) ##lookup svd interface again + ϕ_1=ψ[n1] + ϕ_2=U + bondtensor = S * V + + elseif rlim==n2 + @assert llim == n1-1 + U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=0., kwargs...) + ϕ_1=U + ϕ_2=ψ[n2] + bondtensor = S * V end + position!(PH,ψ,min(n1,n2)) + @show PH.lpos + @show PH.rpos + @show ψ.llim + @show ψ.rlim + @show llim, rlim + @show n1,n2 - newL,S,newR=_subspace_expand_core(ϕ,PH,NL,NR,;maxdim, cutoff, atol=1e-2, kwargs...) - nullbond=diagITensor(0.0,uniqueinds(newL, ψ[n1]),uniqueinds(newR, ψ[n2])) - + + #orthogonalize(ψ,n1) + linkind_l=commonind(ϕ_1,bondtensor) + linkind_r=commonind(ϕ_2,bondtensor) + + NL=nullspace(ϕ_1,linkind_l;atol=1e-9) + NR=nullspace(ϕ_2,linkind_r;atol=1e-9) + #@show norm(NL) + #@show norm(NR) + ###this is a crucial decision, justified for MPS, but will fail for rank-1 trees with physical DOFs on leafs + ###vanishing norm should trigger a one-sided subspace expansion + if norm(NL)==0.0 || norm(NR)==0. + return bondtensor + end + + ϕ=ϕ_1 * bondtensor * ϕ_2 + #ψ[n1]=prime(ψ[n1],linkind) + + + newL,S,newR,success=_subspace_expand_core(ϕ,PH,NL,NR,;maxdim, cutoff, atol=1e-2, kwargs...) + if success == false + return nothing + end ALⁿ¹, newl = ITensors.directsum( - ψ[n1], dag(newL), uniqueinds(ψ[n1], newL), uniqueinds(newL, ψ[n1]); tags=("Left",) + ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) ) ARⁿ², newr = ITensors.directsum( - ψ.AR[n2], dag(newR), uniqueinds(ψ[n2], newR), uniqueinds(newR, ψ[n2]); tags=("Right",) + ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) ) - C,newlr=ITensors.directsum(bondtensor => (prime(linkind,1), linkind), nullbond => (uniqueinds(newL, ψ[n1]), uniqueinds(newR, ψ[n2])),tags=("Left","Right")) - + #C,newlr=ITensors.directsum(bondtensor => (prime(linkind_l,1), linkind_r), nullbond => (uniqueinds(newL, ϕ_1), uniqueinds(newR, ϕ_2)),tags=("Left","Right")) + C = ITensor(dag(newl)..., dag(newr)...) + ψC = permute(bondtensor, linkind_l, linkind_r) + for I in eachindex(ψC) + v = ψC[I] + if !iszero(v) + C[I] = ψC[I] + end + end + println("before") + @show inds(ψ[n1]) + @show inds(ψ[n2]) if rlim==n2 ψ[n2]=ARⁿ² elseif rlim>n2 @@ -79,6 +120,10 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b elseif llim Date: Tue, 29 Mar 2022 16:52:59 -0400 Subject: [PATCH 05/19] Fixed subspace expansion on backward sweep, sets orthogonality limits manually for 1-site. 2-site TDVP currently broken. --- examples/tdvp.jl | 42 +++------------- src/subspace_expansion.jl | 102 +++++++++++++++++++++----------------- src/tdvp.jl | 12 ++++- 3 files changed, 75 insertions(+), 81 deletions(-) diff --git a/examples/tdvp.jl b/examples/tdvp.jl index 897e91b..98c9d27 100644 --- a/examples/tdvp.jl +++ b/examples/tdvp.jl @@ -1,6 +1,6 @@ using ITensors using ITensorTDVP -n = 10 +n = 25 s = siteinds("S=1/2", n,conserve_qns=true) function heisenberg(n) @@ -18,48 +18,22 @@ H = MPO(heisenberg(n), s) E0= inner(ψ', H, ψ) / inner(ψ, ψ) + ϕ = tdvp( H, ψ, - 0.2im; - nsweeps=1, - nsite=1, - reverse_step=true, - normalize=true, - maxdim=50, - cutoff=1e-10, - outputlevel=1, -) -ϕ = tdvp( - H, - ϕ, - 0.2im; - nsweeps=1, + -0.1im; + nsweeps=20, nsite=1, reverse_step=true, - normalize=true, + normalize=false, maxdim=50, - cutoff=1e-10, - outputlevel=1, -) -ϕ = tdvp( - H, - ϕ, - 0.2im; - nsweeps=2, - nsite=1, - reverse_step=true, - normalize=true, - maxdim=100, - cutoff=1e-10, + cutoff=5e-2, + atol=1e-11, outputlevel=1, ) @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ) , E0 -e2, ϕ2 = dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10) - -@show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2) - -ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10) +ϕ3 = ITensors.dmrg(H, ψ; nsweeps=10, maxdim=50, cutoff=1e-10) @show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3) diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 1547219..b39872d 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -23,7 +23,8 @@ function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, nsite=2 position!(PH, ψ, 1) for (b, ha) in sweepnext(N; ncenter=2) - println(b) + + ##TODO: figure out whether these calls should be here or inside subspace expansion, currently we do both? orthogonalize!(ψ,b) position!(PH, ψ, b) @@ -31,7 +32,7 @@ function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, b1 = (ha == 1 ? b + 1 : b) Δ = (ha == 1 ? +1 : -1) inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) - _=subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... + subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... ) end end @@ -43,72 +44,89 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) llim,rlim = lims n1, n2 = b - #@assert n1 + 1 == n2 || n1 -1 == n2 PH.nsite=2 + old_linkdim=dim(commonind(ψ[n1],ψ[n2])) + + ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly + ###the cutoff should be scaled with timestep, otherwise one runs into problems with non-monotonic error behaviour like in TEBD approaches if llim == n1 @assert rlim == n2+1 - U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=0., kwargs...) ##lookup svd interface again - ϕ_1=ψ[n1] + U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=1e-17, kwargs...) ##lookup svd interface again + ϕ_1=ψ[n1] * V ϕ_2=U - bondtensor = S * V + old_linkdim=dim(commonind(U,S)) + bondtensor = S elseif rlim==n2 @assert llim == n1-1 - U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=0., kwargs...) + U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=1e-17, kwargs...) ϕ_1=U - ϕ_2=ψ[n2] - bondtensor = S * V + ϕ_2=ψ[n2] * V + old_linkdim=dim(commonind(U,S)) + bondtensor = S + end + + ###don't expand if we are already at maxdim + if old_linkdim>=maxdim + println("not expanding") + return nothing end position!(PH,ψ,min(n1,n2)) - @show PH.lpos - @show PH.rpos - @show ψ.llim - @show ψ.rlim - @show llim, rlim - @show n1,n2 #orthogonalize(ψ,n1) linkind_l=commonind(ϕ_1,bondtensor) linkind_r=commonind(ϕ_2,bondtensor) - NL=nullspace(ϕ_1,linkind_l;atol=1e-9) - NR=nullspace(ϕ_2,linkind_r;atol=1e-9) + NL=nullspace(ϕ_1,linkind_l;atol=atol) + NR=nullspace(ϕ_2,linkind_r;atol=atol) - #@show norm(NL) - #@show norm(NR) - ###this is a crucial decision, justified for MPS, but will fail for rank-1 trees with physical DOFs on leafs - ###vanishing norm should trigger a one-sided subspace expansion + ###NOTE: This will fail for rank-1 trees with physical DOFs on leafs + ###NOTE: one-sided subspace expansion seems to not work well at least for trees according to Lachlan Lindoy if norm(NL)==0.0 || norm(NR)==0. return bondtensor end + ###form 2site wavefunction ϕ=ϕ_1 * bondtensor * ϕ_2 - #ψ[n1]=prime(ψ[n1],linkind) - - newL,S,newR,success=_subspace_expand_core(ϕ,PH,NL,NR,;maxdim, cutoff, atol=1e-2, kwargs...) + ###get subspace expansion + newL,S,newR,success=_subspace_expand_core(ϕ,PH,NL,NR,;maxdim=maxdim-old_linkdim, cutoff, kwargs...) if success == false return nothing end + + ###add expansion direction to current site tensors ALⁿ¹, newl = ITensors.directsum( ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) ) ARⁿ², newr = ITensors.directsum( ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) ) - #C,newlr=ITensors.directsum(bondtensor => (prime(linkind_l,1), linkind_r), nullbond => (uniqueinds(newL, ϕ_1), uniqueinds(newR, ϕ_2)),tags=("Left","Right")) + + ###TODO remove assertions regarding expansion not exceeding maxdim + @assert (dim(commonind(newL,S))+old_linkdim) <=maxdim + @assert dim(commonind(newL,S))==dim(commonind(newR,S)) + @assert(dim(uniqueind(ϕ_1, newL))+dim(uniqueind(newL, ϕ_1))==dim(newl)) + @assert(dim(uniqueind(ϕ_2, newR))+dim(uniqueind(newR, ϕ_2))==dim(newr)) + @assert (old_linkdim + dim(commonind(newL,S))) <=maxdim + @assert (old_linkdim + dim(commonind(newR,S))) <=maxdim + @assert dim(newl)<=maxdim + @assert dim(newr)<=maxdim + + ###zero-pad bond-tensor (the orthogonality center) C = ITensor(dag(newl)..., dag(newr)...) - ψC = permute(bondtensor, linkind_l, linkind_r) + ψC = bondtensor + ### FIXME: the permute below fails, maybe because this already the layout of bondtensor --- in any case it shouldn't fail? + #ψC = permute(bondtensor, linkind_l, linkind_r) for I in eachindex(ψC) v = ψC[I] if !iszero(v) C[I] = ψC[I] end end - println("before") - @show inds(ψ[n1]) - @show inds(ψ[n2]) + + ###move orthogonality center back to site (should restore input orthogonality limits) if rlim==n2 ψ[n2]=ARⁿ² elseif rlim>n2 @@ -120,36 +138,28 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b elseif llim Date: Tue, 29 Mar 2022 21:20:37 -0400 Subject: [PATCH 06/19] Start implementation of krylov/truncated SVD for subspace expansion, not working yet. --- src/ITensorTDVP.jl | 4 +- src/subspace_expansion.jl | 153 +++++++++++++++++++++++++++++++++++++- 2 files changed, 153 insertions(+), 4 deletions(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 9b80157..4f398f1 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,14 +1,14 @@ module ITensorTDVP using ITensors -using KrylovKit: exponentiate, eigsolve +using KrylovKit: exponentiate, eigsolve,svdsolve using Printf using LinearAlgebra using ITensors: AbstractMPS, @debug_check, @timeit_debug, check_hascommoninds, orthocenter using ITensors.NDTensors using ITensors.NDTensors: eachdiagblock, blockview - +using ITensors.ITensorNetworkMaps include("nullspace.jl") include("subspace_expansion.jl") diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index b39872d..028b26e 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -1,4 +1,3 @@ - function replaceind_indval(IV::Tuple, iĩ::Pair) i, ĩ = iĩ return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV)) @@ -32,18 +31,168 @@ function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, b1 = (ha == 1 ? b + 1 : b) Δ = (ha == 1 ? +1 : -1) inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) - subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... + subspace_expansion_krylov!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... ) end end return nothing end +function subspace_expansion_krylov!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) + ##this should only work for the case where rlim-llim > 1 + ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) + llim,rlim = lims + n1, n2 = b + PH.nsite=2 + old_linkdim=dim(commonind(ψ[n1],ψ[n2])) + + ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly + ###the cutoff should be scaled with timestep, otherwise one runs into problems with non-monotonic error behaviour like in TEBD approaches + if llim == n1 + @assert rlim == n2+1 + U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=1e-17, kwargs...) ##lookup svd interface again + ϕ_1=ψ[n1] * V + ϕ_2=U + old_linkdim=dim(commonind(U,S)) + bondtensor = S + + elseif rlim==n2 + @assert llim == n1-1 + U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=1e-17, kwargs...) + ϕ_1=U + ϕ_2=ψ[n2] * V + old_linkdim=dim(commonind(U,S)) + bondtensor = S + end + + ###don't expand if we are already at maxdim + if old_linkdim>=maxdim + println("not expanding") + return nothing + end + position!(PH,ψ,min(n1,n2)) + + + #orthogonalize(ψ,n1) + linkind_l=commonind(ϕ_1,bondtensor) + linkind_r=commonind(ϕ_2,bondtensor) + + NL=nullspace(ϕ_1,linkind_l;atol=atol) + NR=nullspace(ϕ_2,linkind_r;atol=atol) + + ###NOTE: This will fail for rank-1 trees with physical DOFs on leafs + ###NOTE: one-sided subspace expansion seems to not work well at least for trees according to Lachlan Lindoy + if norm(NL)==0.0 || norm(NR)==0. + return bondtensor + end + + ###form 2site wavefunction + ϕ=[ϕ_1,bondtensor, ϕ_2] + + ###get subspace expansion + newL,S,newR,success=_subspace_expand_core_krylov(ϕ,PH,NL,NR,;maxdim=maxdim-old_linkdim, cutoff, kwargs...) + if success == false + return nothing + end + + ###add expansion direction to current site tensors + ALⁿ¹, newl = ITensors.directsum( + ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) + ) + ARⁿ², newr = ITensors.directsum( + ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) + ) + + ###TODO remove assertions regarding expansion not exceeding maxdim + @assert (dim(commonind(newL,S))+old_linkdim) <=maxdim + @assert dim(commonind(newL,S))==dim(commonind(newR,S)) + @assert(dim(uniqueind(ϕ_1, newL))+dim(uniqueind(newL, ϕ_1))==dim(newl)) + @assert(dim(uniqueind(ϕ_2, newR))+dim(uniqueind(newR, ϕ_2))==dim(newr)) + @assert (old_linkdim + dim(commonind(newL,S))) <=maxdim + @assert (old_linkdim + dim(commonind(newR,S))) <=maxdim + @assert dim(newl)<=maxdim + @assert dim(newr)<=maxdim + + ###zero-pad bond-tensor (the orthogonality center) + C = ITensor(dag(newl)..., dag(newr)...) + ψC = bondtensor + ### FIXME: the permute below fails, maybe because this already the layout of bondtensor --- in any case it shouldn't fail? + #ψC = permute(bondtensor, linkind_l, linkind_r) + for I in eachindex(ψC) + v = ψC[I] + if !iszero(v) + C[I] = ψC[I] + end + end + + ###move orthogonality center back to site (should restore input orthogonality limits) + if rlim==n2 + ψ[n2]=ARⁿ² + elseif rlim>n2 + ψ[n2]=noprime(ARⁿ²*C) + end + + if llim==n1 + ψ[n1]=ALⁿ¹ + elseif llim noprime(B2(x)),y -> noprime(B2dag(y))),trial,maxdim) + #TO DO construct U,S,V objects, using only the vals > cutoff, and at most maxdim + @show vals + @show lvecs + @show rvecs + +end + +function _subspace_expand_core(ϕ::ITensor, env,NL,NR;maxdim, cutoff, kwargs...) + ϕH = noprime(env(ϕ)) #add noprime? + ϕH = NL * ϕH * NR + if norm(ϕH) == 0.0 + return false,false,false,false + end + U,S,V=svd(ϕH,commoninds(ϕH,NL);maxdim=maxdim, cutoff=cutoff, kwargs...) + + @assert dim(commonind(U,S))<=maxdim + + + NL *= dag(U) + NR *= dag(V) + return NL,S,NR,true +end + function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) ##this should only work for the case where rlim-llim > 1 ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) llim,rlim = lims n1, n2 = b + #@show n1+1==n2 PH.nsite=2 old_linkdim=dim(commonind(ψ[n1],ψ[n2])) From 061b8d5224385b9c26380567e96f6adf468ad3a1 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 30 Mar 2022 09:57:58 -0400 Subject: [PATCH 07/19] Input to svdsolve added, reconversion back to site-tensors missing. --- src/subspace_expansion.jl | 44 ++++++++++++++++++--------------------- 1 file changed, 20 insertions(+), 24 deletions(-) diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 028b26e..fb61b84 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -3,6 +3,7 @@ function replaceind_indval(IV::Tuple, iĩ::Pair) return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV)) end + # atol controls the tolerance cutoff for determining which eigenvectors are in the null # space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep # the eigenvectors corresponding to eigenvalues of approximately 1. @@ -158,33 +159,28 @@ function _subspace_expand_core_krylov(centerwf::Vector{ITensor}, env,NL,NR;maxdi outinds=uniqueinds(NL,NL) ininds=uniqueinds(NR,NR) B2=ITensorNetworkMap([NR,R,L,NL],ininds,outinds) - B2dag=adjoint(B2) - trial=randomITensor(uniqueinds(NL,L)) #columnspace of B2, i.e. NL - trial=noprime(B2(noprime(B2dag(trial)))) - + B2dag=adjoint(B2) + trial=randomITensor(eltype(L),uniqueinds(NL,L)) + trialadj=randomITensor(eltype(R),uniqueinds(NR,R)) #columnspace of B2, i.e. NL + #trial=noprime(B2(noprime(B2dag(trial)))) + #vals, lvecs, rvecs, info = svdsolve(trial, maxdim) do (x, flag) + # if flag + # y = B2dag * copy(x)# y = compute action of adjoint map on x + # else + # y = B2 * copy(x)# y = compute action of linear map on x + # end + # return y + #end + vals, lvecs, rvecs, info=svdsolve((x -> noprime(B2*x),y -> noprime(B2dag*y)),trial,maxdim) ###seems to work - vals, lvecs, rvecs, info=svdsolve((x -> noprime(B2(x)),y -> noprime(B2dag(y))),trial,maxdim) #TO DO construct U,S,V objects, using only the vals > cutoff, and at most maxdim + # @show vals - @show lvecs - @show rvecs - -end - -function _subspace_expand_core(ϕ::ITensor, env,NL,NR;maxdim, cutoff, kwargs...) - ϕH = noprime(env(ϕ)) #add noprime? - ϕH = NL * ϕH * NR - if norm(ϕH) == 0.0 - return false,false,false,false - end - U,S,V=svd(ϕH,commoninds(ϕH,NL);maxdim=maxdim, cutoff=cutoff, kwargs...) - - @assert dim(commonind(U,S))<=maxdim - - - NL *= dag(U) - NR *= dag(V) - return NL,S,NR,true + @show uniqueinds(NL,L) + @show uniqueinds(NR,R) + @show inds(lvecs[1]) + @show inds(rvecs[1]) + println("success!!") end function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) From ccfd7e58671526563fcd4ad2c74da406b629d878 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Tue, 12 Apr 2022 10:51:38 -0400 Subject: [PATCH 08/19] Moving to new tdvp interface, abandoning implementation in old interface. --- src/tdvp.jl | 323 +++++++--------------------------------------------- 1 file changed, 43 insertions(+), 280 deletions(-) diff --git a/src/tdvp.jl b/src/tdvp.jl index 8e2d550..2bca32f 100644 --- a/src/tdvp.jl +++ b/src/tdvp.jl @@ -1,293 +1,56 @@ -function tdvp(solver, PH, psi0::MPS, t::Number, sweeps; kwargs...)::MPS - if length(psi0) == 1 - error( - "`tdvp` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.exponentiate`, etc.", +function exponentiate_solver(; kwargs...) + function solver(H, t, psi0; kws...) + solver_kwargs = (; + ishermitian=get(kwargs, :ishermitian, true), + issymmetric=get(kwargs, :issymmetric, true), + tol=get(kwargs, :solver_tol, 1E-12), + krylovdim=get(kwargs, :solver_krylovdim, 30), + maxiter=get(kwargs, :solver_maxiter, 100), + verbosity=get(kwargs, :solver_outputlevel, 0), + eager=true, ) + psi, info = exponentiate(H, t, psi0; solver_kwargs..., kws...) + return psi, info end - - @debug_check begin - # Debug level checks - # Enable with ITensors.enable_debug_checks() - checkflux(psi0) - checkflux(PH) - end - - nsite::Int = get(kwargs, :nsite, 2) - reverse_step::Bool = get(kwargs, :reverse_step, true) - normalize::Bool = get(kwargs, :normalize, false) - outputlevel::Int = get(kwargs, :outputlevel, 0) - which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing) - svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer") - obs = get(kwargs, :observer, NoObserver()) - atol = get(kwargs, :atol,1e-6) - maxdim=get(kwargs,:maxdim,100) - write_when_maxdim_exceeds::Union{Int,Nothing} = get( - kwargs, :write_when_maxdim_exceeds, nothing - ) - - psi = copy(psi0) - N = length(psi) - - if !isortho(psi) || orthocenter(psi) != 1 - orthogonalize!(psi, 1) - end - @assert isortho(psi) && orthocenter(psi) == 1 - #PH.nsite=1 - #position!(PH, psi, 1) - #psid=PH * psi - - #println("did a successful env itensor mult") - for sw in 1:nsweep(sweeps) - orthogonalize!(psi, 1) - if nsite==1 - @time subspace_expansion_sweep!(psi,PH;maxdim, cutoff, atol=atol, kwargs...) - end - sw_time = @elapsed begin - maxtruncerr = 0.0 - - if !isnothing(write_when_maxdim_exceeds) && - maxdim(sweeps, sw) > write_when_maxdim_exceeds - if outputlevel >= 2 - println( - "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim(sweeps, sw) = $(maxdim(sweeps, sw)), writing environment tensors to disk", - ) - end - PH = disk(PH) - end - - for (b, ha) in sweepnext(N; ncenter=nsite) - PH.nsite = nsite - position!(PH, psi, b) - - if nsite == 1 - phi1 = psi[b] - elseif nsite == 2 - phi1 = psi[b] * psi[b + 1] - end - #phiprime=PH(phi1) - #println("did a successful env itensor mult") - phi1, info = solver(PH, t / 2, phi1) - - ## if info.converged == 0 - ## println("exponentiate not converged (b,ha)=($b,$ha)") - ## ITensors.pause() - ## end - - normalize && (phi1 /= norm(phi1)) - - spec = nothing - if nsite == 1 - psi[b] = phi1 - elseif nsite == 2 - ortho = ha == 1 ? "left" : "right" - - drho = nothing - if noise(sweeps, sw) > 0.0 && ha == 1 - drho = noise(sweeps, sw) * noiseterm(PH, phi, ortho) - end - - spec = replacebond!( - psi, - b, - phi1; - maxdim=maxdim(sweeps, sw), - mindim=mindim(sweeps, sw), - cutoff=cutoff(sweeps, sw), - eigen_perturbation=drho, - ortho=ortho, - normalize, - which_decomp, - svd_alg, - ) - maxtruncerr = max(maxtruncerr, spec.truncerr) - end - - # - # Do backwards evolution step - # - if reverse_step && (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) - b1 = (ha == 1 ? b + 1 : b) - Δ = (ha == 1 ? +1 : -1) - if nsite == 2 - phi0 = psi[b1] - elseif nsite == 1 - uinds = uniqueinds(phi1, psi[b + Δ]) - U, S, V = svd(phi1, uinds) - psi[b] = U - if ha==1 - ITensors.setleftlim!(psi,b) - elseif ha==2 - ITensors.setrightlim!(psi,b) - end - phi0 = S * V - end - - PH.nsite = nsite - 1 - position!(PH, psi, b1) - - phi0, info = solver(PH, -t / 2, phi0) - - normalize && (phi0 ./= norm(phi0)) - - if nsite == 2 - psi[b1] = phi0 - elseif nsite == 1 - psi[b + Δ] = phi0 * psi[b + Δ] - if ha==1 - ITensors.setrightlim!(psi,b + Δ + 1) - elseif ha==2 - ITensors.setleftlim!(psi,b + Δ - 1) - end - end - PH.nsite = nsite - end - - if outputlevel >= 2 - @printf("Sweep %d, half %d, bond (%d,%d) \n", sw, ha, b, b + 1) - @printf( - " Truncated using cutoff=%.1E maxdim=%d mindim=%d\n", - cutoff(sweeps, sw), - maxdim(sweeps, sw), - mindim(sweeps, sw) - ) - @printf( - " Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(psi, b)) - ) - flush(stdout) - end - - sweep_is_done = (b == 1 && ha == 2) - measure!( - obs; psi, bond=b, sweep=sw, half_sweep=ha, spec, outputlevel, sweep_is_done - ) - end - end #@elapsed for sw_time - - if outputlevel >= 1 - @printf( - "After sweep %d maxlinkdim=%d maxerr=%.2E time=%.3f\n", - sw, - maxlinkdim(psi), - maxtruncerr, - sw_time - ) - flush(stdout) - end - isdone = checkdone!(obs; psi, sweep=sw, outputlevel) - - isdone && break - end - - # Just to be sure: - normalize && normalize!(psi) - - return psi -end - -function _tdvp_sweeps(; - nsweeps=1, maxdim=typemax(Int), mindim=1, cutoff=1E-8, noise=0.0, kwargs... -) - sweeps = Sweeps(nsweeps) - setmaxdim!(sweeps, maxdim...) - setmindim!(sweeps, mindim...) - setcutoff!(sweeps, cutoff...) - setnoise!(sweeps, noise...) - return sweeps -end - -function tdvp(solver, x1, x2, psi0::MPS, t::Number; kwargs...) - return tdvp(solver, x1, x2, psi0, t, _tdvp_sweeps(; kwargs...); kwargs...) + return solver end -function tdvp(solver, x1, psi0::MPS, t::Number; kwargs...) - return tdvp(solver, x1, psi0, t, _tdvp_sweeps(; kwargs...); kwargs...) +function applyexp_solver(; kwargs...) + function solver(H, t, psi0; kws...) + tol_per_unit_time = get(kwargs, :solver_tol, 1E-8) + solver_kwargs = (; + maxiter=get(kwargs, :solver_krylovdim, 30), + outputlevel=get(kwargs, :solver_outputlevel, 0), + ) + #applyexp tol is absolute, compute from tol_per_unit_time: + tol = abs(t) * tol_per_unit_time + psi, info = applyexp(H, t, psi0; tol, solver_kwargs..., kws...) + return psi, info + end + return solver end -""" - tdvp(H::MPO,psi0::MPS,t::Number,sweeps::Sweeps; kwargs...) - -Use the time dependent variational principle (TDVP) algorithm -to compute `exp(t*H)*psi0` using an efficient algorithm based -on alternating optimization of the MPS tensors and local Krylov -exponentiation of H. - -Returns: -* `psi::MPS` - time-evolved MPS - -Optional keyword arguments: -* `outputlevel::Int = 1` - larger outputlevel values resulting in printing more information and 0 means no output -* `observer` - object implementing the [Observer](@ref observer) interface which can perform measurements and stop early -* `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations -""" -function tdvp(solver, H::MPO, psi0::MPS, t::Number, sweeps::Sweeps; kwargs...)::MPS - check_hascommoninds(siteinds, H, psi0) - check_hascommoninds(siteinds, H, psi0') - # Permute the indices to have a better memory layout - # and minimize permutations - H = ITensors.permute(H, (linkind, siteinds, linkind)) - PH = ProjMPO(H) - return tdvp(solver, PH, psi0, t, sweeps; kwargs...) +function tdvp_solver(; kwargs...) + solver_backend = get(kwargs, :solver_backend, "applyexp") + if solver_backend == "applyexp" + return applyexp_solver(; kwargs...) + elseif solver_backend == "exponentiate" + return exponentiate_solver(; kwargs...) + else + error( + "solver_backend=$solver_backend not recognized (options are \"applyexp\" or \"exponentiate\")", + ) + end end -""" - tdvp(Hs::Vector{MPO},psi0::MPS,t::Number,sweeps::Sweeps;kwargs...) - -Use the time dependent variational principle (TDVP) algorithm -to compute `exp(t*H)*psi0` using an efficient algorithm based -on alternating optimization of the MPS tensors and local Krylov -exponentiation of H. - -This version of `tdvp` accepts a representation of H as a -Vector of MPOs, Hs = [H1,H2,H3,...] such that H is defined -as H = H1+H2+H3+... -Note that this sum of MPOs is not actually computed; rather -the set of MPOs [H1,H2,H3,..] is efficiently looped over at -each step of the algorithm when optimizing the MPS. - -Returns: -* `psi::MPS` - time-evolved MPS -""" -function tdvp(solver, Hs::Vector{MPO}, psi0::MPS, t::Number, sweeps::Sweeps; kwargs...)::MPS - for H in Hs - check_hascommoninds(siteinds, H, psi0) - check_hascommoninds(siteinds, H, psi0') - end - Hs .= ITensors.permute.(Hs, Ref((linkind, siteinds, linkind))) - PHS = ProjMPOSum(Hs) - return tdvp(solver, PHS, psi0, t, sweeps; kwargs...) +function tdvp(H, t::Number, psi0::MPS; kwargs...) + return tdvp(tdvp_solver(; kwargs...), H, t, psi0; kwargs...) end -function tdvp(PH, psi0::MPS, t::Number; reverse_step=true, kwargs...) - solver_kwargs = (; - ishermitian=get(kwargs, :ishermitian, true), - tol=get(kwargs, :exponentiate_tol, 1e-12), - krylovdim=get(kwargs, :exponentiate_krylovdim, 30), - maxiter=get(kwargs, :exponentiate_maxiter, 100), - verbosity=get(kwargs, :exponentiate_verbosity, 0), - ) - function solver(PH, t, psi0) - #psidot=PH(psi0) - psi, info = exponentiate(PH, t, psi0; solver_kwargs...) - return psi, info - end - return tdvp(solver, PH, psi0, t; reverse_step, kwargs...) +function tdvp(t::Number, H, psi0::MPS; kwargs...) + return tdvp(H, t, psi0; kwargs...) end -function dmrg(PH, psi0::MPS; reverse_step=false, kwargs...) - t = Inf # DMRG is TDVP with an infinite timestep and no reverse step - howmany = 1 - which = get(kwargs, :eigsolve_which_eigenvalue, :SR) - solver_kwargs = (; - ishermitian=get(kwargs, :ishermitian, true), - tol=get(kwargs, :eigsolve_tol, 1e-14), - krylovdim=get(kwargs, :eigsolve_krylovdim, 3), - maxiter=get(kwargs, :eigsolve_maxiter, 1), - verbosity=get(kwargs, :eigsolve_verbosity, 0), - ) - function solver(PH, t, psi0) - vals, vecs, info = eigsolve(PH, psi0, howmany, which; solver_kwargs...) - psi = vecs[1] - return psi, info - end - return tdvp(solver, PH, psi0, t; reverse_step, kwargs...) +function tdvp(H, psi0::MPS, t::Number; kwargs...) + return tdvp(H, t, psi0; kwargs...) end From d633c6f2eaf3117b4c91aaee4eb8915904c2e963 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Tue, 26 Apr 2022 18:25:46 -0400 Subject: [PATCH 09/19] separate standard and krylov subspace expansion into different files --- src/subspace_expansion.jl | 2 +- src/subspace_expansion_krylov.jl | 185 +++++++++++++++++++++++++++++++ 2 files changed, 186 insertions(+), 1 deletion(-) create mode 100644 src/subspace_expansion_krylov.jl diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index fb61b84..f44672a 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -32,7 +32,7 @@ function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, b1 = (ha == 1 ? b + 1 : b) Δ = (ha == 1 ? +1 : -1) inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) - subspace_expansion_krylov!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... + subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... ) end end diff --git a/src/subspace_expansion_krylov.jl b/src/subspace_expansion_krylov.jl new file mode 100644 index 0000000..6d0b62f --- /dev/null +++ b/src/subspace_expansion_krylov.jl @@ -0,0 +1,185 @@ +function replaceind_indval(IV::Tuple, iĩ::Pair) + i, ĩ = iĩ + return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV)) + end + + + # atol controls the tolerance cutoff for determining which eigenvectors are in the null + # space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep + # the eigenvectors corresponding to eigenvalues of approximately 1. + + ###the most general implementation may be passing MPS, lims, and bondtensor (defaulting to the Id-matrix) + ###then we can tensorsum both left and right tensor, and just apply the bondtensor to restore input gauge + ###ideally we would also be able to apply tensorsum to the bond tensor instead of that loop below + """ + expand subspace (2-site like) in a sweep + """ + function subspace_expansion_krylov_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, cutoff, atol=1e-2, kwargs...) + N = length(ψ) + if !isortho(ψ) || orthocenter(ψ) != 1 + orthogonalize!(ψ, 1) + end + PH.nsite=2 + nsite=2 + position!(PH, ψ, 1) + for (b, ha) in sweepnext(N; ncenter=2) + + ##TODO: figure out whether these calls should be here or inside subspace expansion, currently we do both? + orthogonalize!(ψ,b) + position!(PH, ψ, b) + + if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) + b1 = (ha == 1 ? b + 1 : b) + Δ = (ha == 1 ? +1 : -1) + inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) + subspace_expansion_krylov!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... + ) + end + end + return nothing + end + + function subspace_expansion_krylov!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) + ##this should only work for the case where rlim-llim > 1 + ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) + llim,rlim = lims + n1, n2 = b + PH.nsite=2 + old_linkdim=dim(commonind(ψ[n1],ψ[n2])) + + ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly + ###the cutoff should be scaled with timestep, otherwise one runs into problems with non-monotonic error behaviour like in TEBD approaches + if llim == n1 + @assert rlim == n2+1 + U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=1e-17, kwargs...) ##lookup svd interface again + ϕ_1=ψ[n1] * V + ϕ_2=U + old_linkdim=dim(commonind(U,S)) + bondtensor = S + + elseif rlim==n2 + @assert llim == n1-1 + U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=1e-17, kwargs...) + ϕ_1=U + ϕ_2=ψ[n2] * V + old_linkdim=dim(commonind(U,S)) + bondtensor = S + end + + ###don't expand if we are already at maxdim + if old_linkdim>=maxdim + println("not expanding") + return nothing + end + position!(PH,ψ,min(n1,n2)) + + + #orthogonalize(ψ,n1) + linkind_l=commonind(ϕ_1,bondtensor) + linkind_r=commonind(ϕ_2,bondtensor) + + NL=nullspace(ϕ_1,linkind_l;atol=atol) + NR=nullspace(ϕ_2,linkind_r;atol=atol) + + ###NOTE: This will fail for rank-1 trees with physical DOFs on leafs + ###NOTE: one-sided subspace expansion seems to not work well at least for trees according to Lachlan Lindoy + if norm(NL)==0.0 || norm(NR)==0. + return bondtensor + end + + ###form 2site wavefunction + ϕ=[ϕ_1,bondtensor, ϕ_2] + + ###get subspace expansion + newL,S,newR,success=_subspace_expand_core_krylov(ϕ,PH,NL,NR,;maxdim=maxdim-old_linkdim, cutoff, kwargs...) + if success == false + return nothing + end + + ###add expansion direction to current site tensors + ALⁿ¹, newl = ITensors.directsum( + ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) + ) + ARⁿ², newr = ITensors.directsum( + ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) + ) + + ###TODO remove assertions regarding expansion not exceeding maxdim + @assert (dim(commonind(newL,S))+old_linkdim) <=maxdim + @assert dim(commonind(newL,S))==dim(commonind(newR,S)) + @assert(dim(uniqueind(ϕ_1, newL))+dim(uniqueind(newL, ϕ_1))==dim(newl)) + @assert(dim(uniqueind(ϕ_2, newR))+dim(uniqueind(newR, ϕ_2))==dim(newr)) + @assert (old_linkdim + dim(commonind(newL,S))) <=maxdim + @assert (old_linkdim + dim(commonind(newR,S))) <=maxdim + @assert dim(newl)<=maxdim + @assert dim(newr)<=maxdim + + ###zero-pad bond-tensor (the orthogonality center) + C = ITensor(dag(newl)..., dag(newr)...) + ψC = bondtensor + ### FIXME: the permute below fails, maybe because this already the layout of bondtensor --- in any case it shouldn't fail? + #ψC = permute(bondtensor, linkind_l, linkind_r) + for I in eachindex(ψC) + v = ψC[I] + if !iszero(v) + C[I] = ψC[I] + end + end + + ###move orthogonality center back to site (should restore input orthogonality limits) + if rlim==n2 + ψ[n2]=ARⁿ² + elseif rlim>n2 + ψ[n2]=noprime(ARⁿ²*C) + end + + if llim==n1 + ψ[n1]=ALⁿ¹ + elseif llim noprime(B2*x),y -> noprime(B2dag*y)),trial,maxdim) ###seems to work + + #TO DO construct U,S,V objects, using only the vals > cutoff, and at most maxdim + # + @show vals + @show uniqueinds(NL,L) + @show uniqueinds(NR,R) + @show inds(lvecs[1]) + @show inds(rvecs[1]) + println("success!!") + end + \ No newline at end of file From 6ba87a8d18c46ae847e7ad70d4c42683d73080a1 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 27 Apr 2022 12:20:37 -0400 Subject: [PATCH 10/19] Adding subspace expansion to every 1-site TDVP timestep and test. --- src/subspace_expansion.jl | 148 ++------------------------------------ src/tdvp_step.jl | 25 ++++++- test/test_tdvp.jl | 105 ++++++++++++++++++++++++++- 3 files changed, 132 insertions(+), 146 deletions(-) diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index f44672a..06b6a97 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -32,6 +32,8 @@ function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, b1 = (ha == 1 ? b + 1 : b) Δ = (ha == 1 ? +1 : -1) inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) + + #println("expanding", b) subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... ) end @@ -39,150 +41,6 @@ function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, return nothing end -function subspace_expansion_krylov!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) - ##this should only work for the case where rlim-llim > 1 - ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) - llim,rlim = lims - n1, n2 = b - PH.nsite=2 - old_linkdim=dim(commonind(ψ[n1],ψ[n2])) - - ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly - ###the cutoff should be scaled with timestep, otherwise one runs into problems with non-monotonic error behaviour like in TEBD approaches - if llim == n1 - @assert rlim == n2+1 - U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=1e-17, kwargs...) ##lookup svd interface again - ϕ_1=ψ[n1] * V - ϕ_2=U - old_linkdim=dim(commonind(U,S)) - bondtensor = S - - elseif rlim==n2 - @assert llim == n1-1 - U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=1e-17, kwargs...) - ϕ_1=U - ϕ_2=ψ[n2] * V - old_linkdim=dim(commonind(U,S)) - bondtensor = S - end - - ###don't expand if we are already at maxdim - if old_linkdim>=maxdim - println("not expanding") - return nothing - end - position!(PH,ψ,min(n1,n2)) - - - #orthogonalize(ψ,n1) - linkind_l=commonind(ϕ_1,bondtensor) - linkind_r=commonind(ϕ_2,bondtensor) - - NL=nullspace(ϕ_1,linkind_l;atol=atol) - NR=nullspace(ϕ_2,linkind_r;atol=atol) - - ###NOTE: This will fail for rank-1 trees with physical DOFs on leafs - ###NOTE: one-sided subspace expansion seems to not work well at least for trees according to Lachlan Lindoy - if norm(NL)==0.0 || norm(NR)==0. - return bondtensor - end - - ###form 2site wavefunction - ϕ=[ϕ_1,bondtensor, ϕ_2] - - ###get subspace expansion - newL,S,newR,success=_subspace_expand_core_krylov(ϕ,PH,NL,NR,;maxdim=maxdim-old_linkdim, cutoff, kwargs...) - if success == false - return nothing - end - - ###add expansion direction to current site tensors - ALⁿ¹, newl = ITensors.directsum( - ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) - ) - ARⁿ², newr = ITensors.directsum( - ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) - ) - - ###TODO remove assertions regarding expansion not exceeding maxdim - @assert (dim(commonind(newL,S))+old_linkdim) <=maxdim - @assert dim(commonind(newL,S))==dim(commonind(newR,S)) - @assert(dim(uniqueind(ϕ_1, newL))+dim(uniqueind(newL, ϕ_1))==dim(newl)) - @assert(dim(uniqueind(ϕ_2, newR))+dim(uniqueind(newR, ϕ_2))==dim(newr)) - @assert (old_linkdim + dim(commonind(newL,S))) <=maxdim - @assert (old_linkdim + dim(commonind(newR,S))) <=maxdim - @assert dim(newl)<=maxdim - @assert dim(newr)<=maxdim - - ###zero-pad bond-tensor (the orthogonality center) - C = ITensor(dag(newl)..., dag(newr)...) - ψC = bondtensor - ### FIXME: the permute below fails, maybe because this already the layout of bondtensor --- in any case it shouldn't fail? - #ψC = permute(bondtensor, linkind_l, linkind_r) - for I in eachindex(ψC) - v = ψC[I] - if !iszero(v) - C[I] = ψC[I] - end - end - - ###move orthogonality center back to site (should restore input orthogonality limits) - if rlim==n2 - ψ[n2]=ARⁿ² - elseif rlim>n2 - ψ[n2]=noprime(ARⁿ²*C) - end - - if llim==n1 - ψ[n1]=ALⁿ¹ - elseif llim noprime(B2*x),y -> noprime(B2dag*y)),trial,maxdim) ###seems to work - - #TO DO construct U,S,V objects, using only the vals > cutoff, and at most maxdim - # - @show vals - @show uniqueinds(NL,L) - @show uniqueinds(NR,R) - @show inds(lvecs[1]) - @show inds(rvecs[1]) - println("success!!") -end - function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) ##this should only work for the case where rlim-llim > 1 ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) @@ -237,6 +95,8 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b ###get subspace expansion newL,S,newR,success=_subspace_expand_core(ϕ,PH,NL,NR,;maxdim=maxdim-old_linkdim, cutoff, kwargs...) + @assert success + #@show expanS if success == false return nothing end diff --git a/src/tdvp_step.jl b/src/tdvp_step.jl index f9bc797..e589612 100644 --- a/src/tdvp_step.jl +++ b/src/tdvp_step.jl @@ -5,6 +5,29 @@ function tdvp( sub_time_steps = ITensorTDVP.sub_time_steps(order) sub_time_steps *= time_step global info + + #psi_evolved = copy(psi) #copy required, otherwise this tdvp function mutates psi in the loop (and subspace expansion if performed) + + nsite::Int = get(kwargs, :nsite, 2) + expand::Bool = get(kwargs, :expand, nsite == 2 ? false : (hasqns(psi) ? true : false )) + ###TODO: Pass accuracy and other control parameters for subspace expansion instead of hardcoding + if expand + ###TODO: check that shallow copy leaves psi invariant in the following 'in-place' subspace expansion + maxdim = get(kwargs, :maxdim, 100) + @show kwargs + @show expand + @show nsite + @show maxdim + cutoff = get(kwargs, :cutoff, 1e-12) + atol = get(kwargs, :cutoff, 1e-12) + println("running subspace expansion sweep") + println("maxlinkdim before") + @show maxlinkdim(psi) + ITensorTDVP.subspace_expansion_sweep!(psi, PH;maxdim, cutoff, atol=atol) + println("maxlinkdim after") + @show maxlinkdim(psi) + end + for substep in 1:length(sub_time_steps) psi, PH, info = tdvp( orderings[substep], solver, PH, sub_time_steps[substep], psi; current_time, kwargs... @@ -15,7 +38,7 @@ function tdvp( end function tdvp(direction::Base.Ordering, solver, PH, time_step::Number, psi::MPS; kwargs...) - PH = copy(PH) + #PH = copy(PH) ###TODO: implement copy for PH, otherwhise this function may mutate PH psi = copy(psi) if length(psi) == 1 diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl index d366ff0..c4b8490 100644 --- a/test/test_tdvp.jl +++ b/test/test_tdvp.jl @@ -8,7 +8,7 @@ using Test @testset "Basic TDVP" begin N = 10 cutoff = 1e-12 - + #maxdim = 100 s = siteinds("S=1/2", N) os = OpSum() @@ -236,6 +236,7 @@ end #normalize!(psi) nsite = (step <= 3 ? 2 : 1) + println("step: ", step, " doing nsite ", nsite) phi = tdvp( H, -tau * im, phi; nsweeps=1, cutoff, nsite, normalize=true, exponentiate_krylovdim=15 ) @@ -266,6 +267,7 @@ end H, -im * ttotal, phi; + nsite=2, time_step=-im * tau, cutoff, normalize=false, @@ -283,6 +285,107 @@ end @test norm(En1 - En2) < 1e-3 end +@testset "Subspace expansions vs 2-site evolution" begin + N = 10 + cutoff = 1e-12 + tau = 0.1 + ttotal = 1.0 + + s = siteinds("S=1/2", N; conserve_qns=true) + + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + + H = MPO(os, s) + + + psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn") + phi = copy(psi) + c = div(N, 2) + + + # + # Evolve using TDVP + # + struct TDVPObserver <: AbstractObserver end + + Nsteps = convert(Int, ceil(abs(ttotal / tau))) + Sz1 = zeros(Nsteps) + En1 = zeros(Nsteps) + Sz2 = zeros(Nsteps) + En2 = zeros(Nsteps) + + function ITensors.measure!(obs::TDVPObserver; sweep, bond, half_sweep, psi, kwargs...) + if bond == 1 && half_sweep == 2 + Sz1[sweep] = expect(psi, "Sz"; sites=c) + En1[sweep] = real(inner(psi', H, psi)) + end + end + + phi = productMPS(s, n -> isodd(n) ? "Up" : "Dn") + + phi = tdvp( + H, + -im * tau, + phi; + nsite=2, + time_step=-im * tau, + cutoff, + normalize=false, + ) + + psi=deepcopy(phi) + + phi = tdvp( + H, + -im * ttotal, + phi; + nsite=2, + time_step=-im * tau, + cutoff, + normalize=false, + (observer!)=TDVPObserver(), + ) + + function ITensors.measure!(obs::TDVPObserver; sweep, bond, half_sweep, psi, kwargs...) + if bond == 1 && half_sweep == 2 + Sz2[sweep] = expect(psi, "Sz"; sites=c) + En2[sweep] = real(inner(psi', H, psi)) + end + end + + psi = tdvp( + H, + -im * ttotal, + psi; + time_step=-im * tau, + #cutoff, + normalize=false, + nsite=1, + maxdim=100, + expand=true, + #cutoff=5e-2, + #atol=1e-11, + (observer!)=TDVPObserver(), + ) + @show maxlinkdim(psi) + #display(En1) + #display(En2) + #display(Sz1) + #display(Sz2) + #@show norm(Sz1 - Sz2) + #@show norm(En1 - En2) + @show(Sz1) + @show(Sz2) + + @test norm(Sz1 - Sz2) < 1e-3 + @test norm(En1 - En2) < 1e-3 +end + @testset "Imaginary Time Evolution" begin N = 10 cutoff = 1e-12 From 1425746d296329a53798453e8cf835097c3fc2c5 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 27 Apr 2022 15:26:59 -0400 Subject: [PATCH 11/19] Add rescaling of cutoff parameters for subspace expansion. Clean-up output and test. --- src/subspace_expansion.jl | 34 +++++++++++++++++++++------------- src/tdvp_step.jl | 21 +++++++-------------- test/test_tdvp.jl | 24 ++++++++++++++---------- 3 files changed, 42 insertions(+), 37 deletions(-) diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 06b6a97..c74032c 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -14,14 +14,17 @@ end """ expand subspace (2-site like) in a sweep """ -function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, cutoff, atol=1e-2, kwargs...) +function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, cutoff, atol=1e-10, kwargs...) N = length(ψ) + if !isortho(ψ) || orthocenter(ψ) != 1 orthogonalize!(ψ, 1) end + PH.nsite=2 nsite=2 position!(PH, ψ, 1) + for (b, ha) in sweepnext(N; ncenter=2) ##TODO: figure out whether these calls should be here or inside subspace expansion, currently we do both? @@ -32,29 +35,36 @@ function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, b1 = (ha == 1 ? b + 1 : b) Δ = (ha == 1 ? +1 : -1) inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) - - #println("expanding", b) - subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... + subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff=cutoff, atol=atol, kwargs... ) end end return nothing end -function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) +function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-10, kwargs...) + """ + expands in nullspace of site-tensors b + """ + #atol refers to the tolerance in nullspace determination (for finite MPS can probably be set rather small) + #cutoff refers to largest singular value of gradient (acceleration of population gain of expansion vectors) to keep + #this should be related to the cutoff in two-site TDVP: \Delta_rho_i = 0.5 * lambda_y * tau **2 + #note that in the initial SVD there is another cutoff parameter, that we set to roughly machine precision for the time being + #(allows to move bond-dimension between different partial QN sectors, if no such truncation occurs distribution of bond-dimensions + #between different QNs locally is static once bond dimension saturates maxdim.) + ##this should only work for the case where rlim-llim > 1 ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) llim,rlim = lims n1, n2 = b - #@show n1+1==n2 PH.nsite=2 old_linkdim=dim(commonind(ψ[n1],ψ[n2])) + cutoff_compress=get(kwargs,:cutoff_compress,1e-12) ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly - ###the cutoff should be scaled with timestep, otherwise one runs into problems with non-monotonic error behaviour like in TEBD approaches if llim == n1 @assert rlim == n2+1 - U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=1e-17, kwargs...) ##lookup svd interface again + U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=cutoff_compress, kwargs...) ##lookup svd interface again ϕ_1=ψ[n1] * V ϕ_2=U old_linkdim=dim(commonind(U,S)) @@ -62,7 +72,7 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b elseif rlim==n2 @assert llim == n1-1 - U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=1e-17, kwargs...) + U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=cutoff_compress, kwargs...) ϕ_1=U ϕ_2=ψ[n2] * V old_linkdim=dim(commonind(U,S)) @@ -109,7 +119,7 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) ) - ###TODO remove assertions regarding expansion not exceeding maxdim + ###TODO remove assertions regarding expansion not exceeding maxdim ? @assert (dim(commonind(newL,S))+old_linkdim) <=maxdim @assert dim(commonind(newL,S))==dim(commonind(newR,S)) @assert(dim(uniqueind(ϕ_1, newL))+dim(uniqueind(newL, ϕ_1))==dim(newl)) @@ -122,8 +132,6 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b ###zero-pad bond-tensor (the orthogonality center) C = ITensor(dag(newl)..., dag(newr)...) ψC = bondtensor - ### FIXME: the permute below fails, maybe because this already the layout of bondtensor --- in any case it shouldn't fail? - #ψC = permute(bondtensor, linkind_l, linkind_r) for I in eachindex(ψC) v = ψC[I] if !iszero(v) @@ -161,7 +169,7 @@ function _subspace_expand_core(ϕ::ITensor, env,NL,NR;maxdim, cutoff, kwargs...) if norm(ϕH) == 0.0 return false,false,false,false end - U,S,V=svd(ϕH,commoninds(ϕH,NL);maxdim=maxdim, cutoff=cutoff, kwargs...) + U,S,V=svd(ϕH,commoninds(ϕH,NL);maxdim, cutoff,use_relative_cutoff = false, use_absolute_cutoff = true, kwargs...) @assert dim(commonind(U,S))<=maxdim diff --git a/src/tdvp_step.jl b/src/tdvp_step.jl index e589612..8f03328 100644 --- a/src/tdvp_step.jl +++ b/src/tdvp_step.jl @@ -10,23 +10,15 @@ function tdvp( nsite::Int = get(kwargs, :nsite, 2) expand::Bool = get(kwargs, :expand, nsite == 2 ? false : (hasqns(psi) ? true : false )) - ###TODO: Pass accuracy and other control parameters for subspace expansion instead of hardcoding + if expand ###TODO: check that shallow copy leaves psi invariant in the following 'in-place' subspace expansion maxdim = get(kwargs, :maxdim, 100) - @show kwargs - @show expand - @show nsite - @show maxdim - cutoff = get(kwargs, :cutoff, 1e-12) - atol = get(kwargs, :cutoff, 1e-12) - println("running subspace expansion sweep") - println("maxlinkdim before") - @show maxlinkdim(psi) - ITensorTDVP.subspace_expansion_sweep!(psi, PH;maxdim, cutoff, atol=atol) - println("maxlinkdim after") - @show maxlinkdim(psi) - end + cutoff = get(kwargs, :cutoff, 1e-11) + atol = get(kwargs, :atol, 1e-12) + cutoff_trunc = 0.5 * abs(time_step)^2 * cutoff ### cutoff_trunc is min. acceleration of population gain of expansion vectors, thus rescaling with timestep + ITensorTDVP.subspace_expansion_sweep!(psi, PH;maxdim, cutoff=cutoff_trunc, atol=atol) + end for substep in 1:length(sub_time_steps) psi, PH, info = tdvp( @@ -34,6 +26,7 @@ function tdvp( ) current_time += sub_time_steps[substep] end + #@show maxlinkdim(psi) return psi, PH, info end diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl index c4b8490..531df1c 100644 --- a/test/test_tdvp.jl +++ b/test/test_tdvp.jl @@ -236,7 +236,6 @@ end #normalize!(psi) nsite = (step <= 3 ? 2 : 1) - println("step: ", step, " doing nsite ", nsite) phi = tdvp( H, -tau * im, phi; nsweeps=1, cutoff, nsite, normalize=true, exponentiate_krylovdim=15 ) @@ -286,10 +285,10 @@ end end @testset "Subspace expansions vs 2-site evolution" begin - N = 10 + N = 16 cutoff = 1e-12 - tau = 0.1 - ttotal = 1.0 + tau = 0.05 + ttotal = 2.0 s = siteinds("S=1/2", N; conserve_qns=true) @@ -314,6 +313,9 @@ end struct TDVPObserver <: AbstractObserver end Nsteps = convert(Int, ceil(abs(ttotal / tau))) + maxdims1 = zeros(Nsteps) + maxdims2 = zeros(Nsteps) + Sz1 = zeros(Nsteps) En1 = zeros(Nsteps) Sz2 = zeros(Nsteps) @@ -323,6 +325,7 @@ end if bond == 1 && half_sweep == 2 Sz1[sweep] = expect(psi, "Sz"; sites=c) En1[sweep] = real(inner(psi', H, psi)) + maxdims1[sweep] = maxlinkdim(psi) end end @@ -346,7 +349,7 @@ end phi; nsite=2, time_step=-im * tau, - cutoff, + cutoff = cutoff, normalize=false, (observer!)=TDVPObserver(), ) @@ -355,32 +358,33 @@ end if bond == 1 && half_sweep == 2 Sz2[sweep] = expect(psi, "Sz"; sites=c) En2[sweep] = real(inner(psi', H, psi)) + maxdims2[sweep] = maxlinkdim(psi) end end - psi = tdvp( H, -im * ttotal, psi; time_step=-im * tau, - #cutoff, + cutoff=cutoff*1e8, ##not clear why this is so much larger than the 2-site TDVP cutoff for comparable bond dimensions? + cutoff_compress=cutoff, normalize=false, nsite=1, maxdim=100, expand=true, + atol=1e-10, #cutoff=5e-2, #atol=1e-11, (observer!)=TDVPObserver(), ) - @show maxlinkdim(psi) #display(En1) #display(En2) #display(Sz1) #display(Sz2) #@show norm(Sz1 - Sz2) #@show norm(En1 - En2) - @show(Sz1) - @show(Sz2) + #@show(Sz1) + #@show(Sz2) @test norm(Sz1 - Sz2) < 1e-3 @test norm(En1 - En2) < 1e-3 From 13ca2f4303801a5c60a37021f0620acae0d73cf2 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 27 Apr 2022 17:36:31 -0400 Subject: [PATCH 12/19] Format. --- src/ITensorTDVP.jl | 6 +- src/subspace_expansion.jl | 172 ++++++++------ src/subspace_expansion_krylov.jl | 370 ++++++++++++++++--------------- src/tdvp_step.jl | 12 +- test/test_tdvp.jl | 26 +-- 5 files changed, 307 insertions(+), 279 deletions(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 50c3058..3a9f538 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,19 +1,15 @@ module ITensorTDVP using ITensors -using KrylovKit: exponentiate, eigsolve,svdsolve +using KrylovKit: exponentiate, eigsolve, svdsolve using Printf using LinearAlgebra - - - using TimerOutputs using Observers using ITensors: AbstractMPS, @debug_check, @timeit_debug, check_hascommoninds, orthocenter, set_nsite! -using ITensors: AbstractMPS, @debug_check, @timeit_debug, check_hascommoninds, orthocenter using ITensors.NDTensors using ITensors.NDTensors: eachdiagblock, blockview using ITensors.ITensorNetworkMaps diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index c74032c..b26fecc 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -3,7 +3,6 @@ function replaceind_indval(IV::Tuple, iĩ::Pair) return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV)) end - # atol controls the tolerance cutoff for determining which eigenvectors are in the null # space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep # the eigenvectors corresponding to eigenvalues of approximately 1. @@ -14,35 +13,48 @@ end """ expand subspace (2-site like) in a sweep """ -function subspace_expansion_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, cutoff, atol=1e-10, kwargs...) +function subspace_expansion_sweep!( + ψ::MPS, PH::Union{ProjMPO,ProjMPOSum}; maxdim, cutoff, atol=1e-10, kwargs... +) N = length(ψ) - + if !isortho(ψ) || orthocenter(ψ) != 1 orthogonalize!(ψ, 1) end - - PH.nsite=2 - nsite=2 + + PH.nsite = 2 + nsite = 2 position!(PH, ψ, 1) for (b, ha) in sweepnext(N; ncenter=2) - + ##TODO: figure out whether these calls should be here or inside subspace expansion, currently we do both? - orthogonalize!(ψ,b) + orthogonalize!(ψ, b) position!(PH, ψ, b) - + if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) b1 = (ha == 1 ? b + 1 : b) Δ = (ha == 1 ? +1 : -1) - inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) - subspace_expansion!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff=cutoff, atol=atol, kwargs... + inds = (ha == 1 ? (b, b + Δ) : (b + Δ, b)) + subspace_expansion!( + ψ, PH, (ψ.llim, ψ.rlim), inds; maxdim, cutoff=cutoff, atol=atol, kwargs... ) end end return nothing end -function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-10, kwargs...) +function subspace_expansion!( + ψ::MPS, + PH, + lims::Tuple{Int,Int}, + b::Tuple{Int,Int}; + bondtensor=nothing, + maxdim, + cutoff, + atol=1e-10, + kwargs..., +) """ expands in nullspace of site-tensors b """ @@ -52,65 +64,70 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b #note that in the initial SVD there is another cutoff parameter, that we set to roughly machine precision for the time being #(allows to move bond-dimension between different partial QN sectors, if no such truncation occurs distribution of bond-dimensions #between different QNs locally is static once bond dimension saturates maxdim.) - + ##this should only work for the case where rlim-llim > 1 ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) - llim,rlim = lims + llim, rlim = lims n1, n2 = b - PH.nsite=2 - old_linkdim=dim(commonind(ψ[n1],ψ[n2])) - cutoff_compress=get(kwargs,:cutoff_compress,1e-12) - + PH.nsite = 2 + old_linkdim = dim(commonind(ψ[n1], ψ[n2])) + cutoff_compress = get(kwargs, :cutoff_compress, 1e-12) + ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly if llim == n1 - @assert rlim == n2+1 - U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=cutoff_compress, kwargs...) ##lookup svd interface again - ϕ_1=ψ[n1] * V - ϕ_2=U - old_linkdim=dim(commonind(U,S)) - bondtensor = S - - elseif rlim==n2 - @assert llim == n1-1 - U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=cutoff_compress, kwargs...) - ϕ_1=U - ϕ_2=ψ[n2] * V - old_linkdim=dim(commonind(U,S)) + @assert rlim == n2 + 1 + U, S, V = svd( + ψ[n2], uniqueinds(ψ[n2], ψ[n1]); maxdim=maxdim, cutoff=cutoff_compress, kwargs... + ) ##lookup svd interface again + ϕ_1 = ψ[n1] * V + ϕ_2 = U + old_linkdim = dim(commonind(U, S)) + bondtensor = S + + elseif rlim == n2 + @assert llim == n1 - 1 + U, S, V = svd( + ψ[n1], uniqueinds(ψ[n1], ψ[n2]); maxdim=maxdim, cutoff=cutoff_compress, kwargs... + ) + ϕ_1 = U + ϕ_2 = ψ[n2] * V + old_linkdim = dim(commonind(U, S)) bondtensor = S end - + ###don't expand if we are already at maxdim - if old_linkdim>=maxdim + if old_linkdim >= maxdim println("not expanding") return nothing end - position!(PH,ψ,min(n1,n2)) - - + position!(PH, ψ, min(n1, n2)) + #orthogonalize(ψ,n1) - linkind_l=commonind(ϕ_1,bondtensor) - linkind_r=commonind(ϕ_2,bondtensor) - - NL=nullspace(ϕ_1,linkind_l;atol=atol) - NR=nullspace(ϕ_2,linkind_r;atol=atol) + linkind_l = commonind(ϕ_1, bondtensor) + linkind_r = commonind(ϕ_2, bondtensor) + + NL = nullspace(ϕ_1, linkind_l; atol=atol) + NR = nullspace(ϕ_2, linkind_r; atol=atol) ###NOTE: This will fail for rank-1 trees with physical DOFs on leafs ###NOTE: one-sided subspace expansion seems to not work well at least for trees according to Lachlan Lindoy - if norm(NL)==0.0 || norm(NR)==0. + if norm(NL) == 0.0 || norm(NR) == 0.0 return bondtensor end - + ###form 2site wavefunction - ϕ=ϕ_1 * bondtensor * ϕ_2 - + ϕ = ϕ_1 * bondtensor * ϕ_2 + ###get subspace expansion - newL,S,newR,success=_subspace_expand_core(ϕ,PH,NL,NR,;maxdim=maxdim-old_linkdim, cutoff, kwargs...) + newL, S, newR, success = _subspace_expand_core( + ϕ, PH, NL, NR, ; maxdim=maxdim - old_linkdim, cutoff, kwargs... + ) @assert success #@show expanS if success == false return nothing end - + ###add expansion direction to current site tensors ALⁿ¹, newl = ITensors.directsum( ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) @@ -120,15 +137,15 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b ) ###TODO remove assertions regarding expansion not exceeding maxdim ? - @assert (dim(commonind(newL,S))+old_linkdim) <=maxdim - @assert dim(commonind(newL,S))==dim(commonind(newR,S)) - @assert(dim(uniqueind(ϕ_1, newL))+dim(uniqueind(newL, ϕ_1))==dim(newl)) - @assert(dim(uniqueind(ϕ_2, newR))+dim(uniqueind(newR, ϕ_2))==dim(newr)) - @assert (old_linkdim + dim(commonind(newL,S))) <=maxdim - @assert (old_linkdim + dim(commonind(newR,S))) <=maxdim - @assert dim(newl)<=maxdim - @assert dim(newr)<=maxdim - + @assert (dim(commonind(newL, S)) + old_linkdim) <= maxdim + @assert dim(commonind(newL, S)) == dim(commonind(newR, S)) + @assert(dim(uniqueind(ϕ_1, newL)) + dim(uniqueind(newL, ϕ_1)) == dim(newl)) + @assert(dim(uniqueind(ϕ_2, newR)) + dim(uniqueind(newR, ϕ_2)) == dim(newr)) + @assert (old_linkdim + dim(commonind(newL, S))) <= maxdim + @assert (old_linkdim + dim(commonind(newR, S))) <= maxdim + @assert dim(newl) <= maxdim + @assert dim(newr) <= maxdim + ###zero-pad bond-tensor (the orthogonality center) C = ITensor(dag(newl)..., dag(newr)...) ψC = bondtensor @@ -140,41 +157,50 @@ function subspace_expansion!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};b end ###move orthogonality center back to site (should restore input orthogonality limits) - if rlim==n2 - ψ[n2]=ARⁿ² - elseif rlim>n2 - ψ[n2]=noprime(ARⁿ²*C) + if rlim == n2 + ψ[n2] = ARⁿ² + elseif rlim > n2 + ψ[n2] = noprime(ARⁿ² * C) end - if llim==n1 - ψ[n1]=ALⁿ¹ - elseif llim first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV)) + i, ĩ = iĩ + return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV)) +end + +# atol controls the tolerance cutoff for determining which eigenvectors are in the null +# space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep +# the eigenvectors corresponding to eigenvalues of approximately 1. + +###the most general implementation may be passing MPS, lims, and bondtensor (defaulting to the Id-matrix) +###then we can tensorsum both left and right tensor, and just apply the bondtensor to restore input gauge +###ideally we would also be able to apply tensorsum to the bond tensor instead of that loop below +""" +expand subspace (2-site like) in a sweep +""" +function subspace_expansion_krylov_sweep!( + ψ::MPS, PH::Union{ProjMPO,ProjMPOSum}; maxdim, cutoff, atol=1e-2, kwargs... +) + N = length(ψ) + if !isortho(ψ) || orthocenter(ψ) != 1 + orthogonalize!(ψ, 1) end - - - # atol controls the tolerance cutoff for determining which eigenvectors are in the null - # space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep - # the eigenvectors corresponding to eigenvalues of approximately 1. - - ###the most general implementation may be passing MPS, lims, and bondtensor (defaulting to the Id-matrix) - ###then we can tensorsum both left and right tensor, and just apply the bondtensor to restore input gauge - ###ideally we would also be able to apply tensorsum to the bond tensor instead of that loop below - """ - expand subspace (2-site like) in a sweep - """ - function subspace_expansion_krylov_sweep!(ψ::MPS,PH::Union{ProjMPO,ProjMPOSum};maxdim, cutoff, atol=1e-2, kwargs...) - N = length(ψ) - if !isortho(ψ) || orthocenter(ψ) != 1 - orthogonalize!(ψ, 1) - end - PH.nsite=2 - nsite=2 - position!(PH, ψ, 1) - for (b, ha) in sweepnext(N; ncenter=2) - - ##TODO: figure out whether these calls should be here or inside subspace expansion, currently we do both? - orthogonalize!(ψ,b) - position!(PH, ψ, b) - - if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) - b1 = (ha == 1 ? b + 1 : b) - Δ = (ha == 1 ? +1 : -1) - inds = (ha == 1 ? (b,b+Δ) : (b+Δ,b)) - subspace_expansion_krylov!(ψ,PH,(ψ.llim,ψ.rlim),inds;maxdim, cutoff, atol=atol, kwargs... - ) - end + PH.nsite = 2 + nsite = 2 + position!(PH, ψ, 1) + for (b, ha) in sweepnext(N; ncenter=2) + + ##TODO: figure out whether these calls should be here or inside subspace expansion, currently we do both? + orthogonalize!(ψ, b) + position!(PH, ψ, b) + + if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) + b1 = (ha == 1 ? b + 1 : b) + Δ = (ha == 1 ? +1 : -1) + inds = (ha == 1 ? (b, b + Δ) : (b + Δ, b)) + subspace_expansion_krylov!( + ψ, PH, (ψ.llim, ψ.rlim), inds; maxdim, cutoff, atol=atol, kwargs... + ) end + end + return nothing +end + +function subspace_expansion_krylov!( + ψ::MPS, + PH, + lims::Tuple{Int,Int}, + b::Tuple{Int,Int}; + bondtensor=nothing, + maxdim, + cutoff, + atol=1e-2, + kwargs..., +) + ##this should only work for the case where rlim-llim > 1 + ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) + llim, rlim = lims + n1, n2 = b + PH.nsite = 2 + old_linkdim = dim(commonind(ψ[n1], ψ[n2])) + + ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly + ###the cutoff should be scaled with timestep, otherwise one runs into problems with non-monotonic error behaviour like in TEBD approaches + if llim == n1 + @assert rlim == n2 + 1 + U, S, V = svd(ψ[n2], uniqueinds(ψ[n2], ψ[n1]); maxdim=maxdim, cutoff=1e-17, kwargs...) ##lookup svd interface again + ϕ_1 = ψ[n1] * V + ϕ_2 = U + old_linkdim = dim(commonind(U, S)) + bondtensor = S + + elseif rlim == n2 + @assert llim == n1 - 1 + U, S, V = svd(ψ[n1], uniqueinds(ψ[n1], ψ[n2]); maxdim=maxdim, cutoff=1e-17, kwargs...) + ϕ_1 = U + ϕ_2 = ψ[n2] * V + old_linkdim = dim(commonind(U, S)) + bondtensor = S + end + + ###don't expand if we are already at maxdim + if old_linkdim >= maxdim + println("not expanding") return nothing end - - function subspace_expansion_krylov!(ψ::MPS,PH,lims::Tuple{Int,Int},b::Tuple{Int,Int};bondtensor=nothing,maxdim, cutoff, atol=1e-2, kwargs...) - ##this should only work for the case where rlim-llim > 1 - ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) - llim,rlim = lims - n1, n2 = b - PH.nsite=2 - old_linkdim=dim(commonind(ψ[n1],ψ[n2])) - - ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly - ###the cutoff should be scaled with timestep, otherwise one runs into problems with non-monotonic error behaviour like in TEBD approaches - if llim == n1 - @assert rlim == n2+1 - U,S,V=svd(ψ[n2],uniqueinds(ψ[n2],ψ[n1]);maxdim=maxdim, cutoff=1e-17, kwargs...) ##lookup svd interface again - ϕ_1=ψ[n1] * V - ϕ_2=U - old_linkdim=dim(commonind(U,S)) - bondtensor = S - - elseif rlim==n2 - @assert llim == n1-1 - U,S,V=svd(ψ[n1],uniqueinds(ψ[n1],ψ[n2]);maxdim=maxdim, cutoff=1e-17, kwargs...) - ϕ_1=U - ϕ_2=ψ[n2] * V - old_linkdim=dim(commonind(U,S)) - bondtensor = S - end - - ###don't expand if we are already at maxdim - if old_linkdim>=maxdim - println("not expanding") - return nothing - end - position!(PH,ψ,min(n1,n2)) - - - #orthogonalize(ψ,n1) - linkind_l=commonind(ϕ_1,bondtensor) - linkind_r=commonind(ϕ_2,bondtensor) - - NL=nullspace(ϕ_1,linkind_l;atol=atol) - NR=nullspace(ϕ_2,linkind_r;atol=atol) - - ###NOTE: This will fail for rank-1 trees with physical DOFs on leafs - ###NOTE: one-sided subspace expansion seems to not work well at least for trees according to Lachlan Lindoy - if norm(NL)==0.0 || norm(NR)==0. - return bondtensor - end - - ###form 2site wavefunction - ϕ=[ϕ_1,bondtensor, ϕ_2] - - ###get subspace expansion - newL,S,newR,success=_subspace_expand_core_krylov(ϕ,PH,NL,NR,;maxdim=maxdim-old_linkdim, cutoff, kwargs...) - if success == false - return nothing - end - - ###add expansion direction to current site tensors - ALⁿ¹, newl = ITensors.directsum( - ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) - ) - ARⁿ², newr = ITensors.directsum( - ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) - ) - - ###TODO remove assertions regarding expansion not exceeding maxdim - @assert (dim(commonind(newL,S))+old_linkdim) <=maxdim - @assert dim(commonind(newL,S))==dim(commonind(newR,S)) - @assert(dim(uniqueind(ϕ_1, newL))+dim(uniqueind(newL, ϕ_1))==dim(newl)) - @assert(dim(uniqueind(ϕ_2, newR))+dim(uniqueind(newR, ϕ_2))==dim(newr)) - @assert (old_linkdim + dim(commonind(newL,S))) <=maxdim - @assert (old_linkdim + dim(commonind(newR,S))) <=maxdim - @assert dim(newl)<=maxdim - @assert dim(newr)<=maxdim - - ###zero-pad bond-tensor (the orthogonality center) - C = ITensor(dag(newl)..., dag(newr)...) - ψC = bondtensor - ### FIXME: the permute below fails, maybe because this already the layout of bondtensor --- in any case it shouldn't fail? - #ψC = permute(bondtensor, linkind_l, linkind_r) - for I in eachindex(ψC) - v = ψC[I] - if !iszero(v) - C[I] = ψC[I] - end - end - - ###move orthogonality center back to site (should restore input orthogonality limits) - if rlim==n2 - ψ[n2]=ARⁿ² - elseif rlim>n2 - ψ[n2]=noprime(ARⁿ²*C) - end - - if llim==n1 - ψ[n1]=ALⁿ¹ - elseif llim noprime(B2*x),y -> noprime(B2dag*y)),trial,maxdim) ###seems to work - - #TO DO construct U,S,V objects, using only the vals > cutoff, and at most maxdim - # - @show vals - @show uniqueinds(NL,L) - @show uniqueinds(NR,R) - @show inds(lvecs[1]) - @show inds(rvecs[1]) - println("success!!") + + ###add expansion direction to current site tensors + ALⁿ¹, newl = ITensors.directsum( + ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) + ) + ARⁿ², newr = ITensors.directsum( + ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) + ) + + ###TODO remove assertions regarding expansion not exceeding maxdim + @assert (dim(commonind(newL, S)) + old_linkdim) <= maxdim + @assert dim(commonind(newL, S)) == dim(commonind(newR, S)) + @assert(dim(uniqueind(ϕ_1, newL)) + dim(uniqueind(newL, ϕ_1)) == dim(newl)) + @assert(dim(uniqueind(ϕ_2, newR)) + dim(uniqueind(newR, ϕ_2)) == dim(newr)) + @assert (old_linkdim + dim(commonind(newL, S))) <= maxdim + @assert (old_linkdim + dim(commonind(newR, S))) <= maxdim + @assert dim(newl) <= maxdim + @assert dim(newr) <= maxdim + + ###zero-pad bond-tensor (the orthogonality center) + C = ITensor(dag(newl)..., dag(newr)...) + ψC = bondtensor + ### FIXME: the permute below fails, maybe because this already the layout of bondtensor --- in any case it shouldn't fail? + #ψC = permute(bondtensor, linkind_l, linkind_r) + for I in eachindex(ψC) + v = ψC[I] + if !iszero(v) + C[I] = ψC[I] + end + end + + ###move orthogonality center back to site (should restore input orthogonality limits) + if rlim == n2 + ψ[n2] = ARⁿ² + elseif rlim > n2 + ψ[n2] = noprime(ARⁿ² * C) + end + + if llim == n1 + ψ[n1] = ALⁿ¹ + elseif llim < n1 + ψ[n1] = noprime(ALⁿ¹ * C) end - \ No newline at end of file + + return nothing +end + +function _subspace_expand_core_krylov( + centerwf::Vector{ITensor}, env, NL, NR; maxdim, cutoff, kwargs... +) + maxdim = min(maxdim, 15) + envL = lproj(env) + envR = rproj(env) + envsiteMPOs = env.H[ITensors.site_range(env)] + ##NOTE: requires n1+1=n2, otherwise assignment of left and right is of + L = envL * centerwf[1] + L *= envsiteMPOs[1] + R = envR * centerwf[3] + R *= envsiteMPOs[2] + C = centerwf[2] + R *= C ##contract C into one of them because otherwise the application is more costly? + R = noprime(R) + L = noprime(L) + outinds = uniqueinds(NL, NL) + ininds = uniqueinds(NR, NR) + B2 = ITensorNetworkMap([NR, R, L, NL], ininds, outinds) + B2dag = adjoint(B2) + trial = randomITensor(eltype(L), uniqueinds(NL, L)) + trialadj = randomITensor(eltype(R), uniqueinds(NR, R)) #columnspace of B2, i.e. NL + #trial=noprime(B2(noprime(B2dag(trial)))) + #vals, lvecs, rvecs, info = svdsolve(trial, maxdim) do (x, flag) + # if flag + # y = B2dag * copy(x)# y = compute action of adjoint map on x + # else + # y = B2 * copy(x)# y = compute action of linear map on x + # end + # return y + #end + vals, lvecs, rvecs, info = svdsolve( + (x -> noprime(B2 * x), y -> noprime(B2dag * y)), trial, maxdim + ) ###seems to work + + #TO DO construct U,S,V objects, using only the vals > cutoff, and at most maxdim + # + @show vals + @show uniqueinds(NL, L) + @show uniqueinds(NR, R) + @show inds(lvecs[1]) + @show inds(rvecs[1]) + return println("success!!") +end diff --git a/src/tdvp_step.jl b/src/tdvp_step.jl index 8f03328..793b06e 100644 --- a/src/tdvp_step.jl +++ b/src/tdvp_step.jl @@ -7,19 +7,19 @@ function tdvp( global info #psi_evolved = copy(psi) #copy required, otherwise this tdvp function mutates psi in the loop (and subspace expansion if performed) - + nsite::Int = get(kwargs, :nsite, 2) - expand::Bool = get(kwargs, :expand, nsite == 2 ? false : (hasqns(psi) ? true : false )) - + expand::Bool = get(kwargs, :expand, nsite == 2 ? false : (hasqns(psi) ? true : false)) + if expand ###TODO: check that shallow copy leaves psi invariant in the following 'in-place' subspace expansion maxdim = get(kwargs, :maxdim, 100) cutoff = get(kwargs, :cutoff, 1e-11) atol = get(kwargs, :atol, 1e-12) cutoff_trunc = 0.5 * abs(time_step)^2 * cutoff ### cutoff_trunc is min. acceleration of population gain of expansion vectors, thus rescaling with timestep - ITensorTDVP.subspace_expansion_sweep!(psi, PH;maxdim, cutoff=cutoff_trunc, atol=atol) - end - + ITensorTDVP.subspace_expansion_sweep!(psi, PH; maxdim, cutoff=cutoff_trunc, atol=atol) + end + for substep in 1:length(sub_time_steps) psi, PH, info = tdvp( orderings[substep], solver, PH, sub_time_steps[substep], psi; current_time, kwargs... diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl index 531df1c..a3279fd 100644 --- a/test/test_tdvp.jl +++ b/test/test_tdvp.jl @@ -301,12 +301,10 @@ end H = MPO(os, s) - psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn") phi = copy(psi) c = div(N, 2) - # # Evolve using TDVP # @@ -315,7 +313,7 @@ end Nsteps = convert(Int, ceil(abs(ttotal / tau))) maxdims1 = zeros(Nsteps) maxdims2 = zeros(Nsteps) - + Sz1 = zeros(Nsteps) En1 = zeros(Nsteps) Sz2 = zeros(Nsteps) @@ -331,25 +329,17 @@ end phi = productMPS(s, n -> isodd(n) ? "Up" : "Dn") - phi = tdvp( - H, - -im * tau, - phi; - nsite=2, - time_step=-im * tau, - cutoff, - normalize=false, - ) - - psi=deepcopy(phi) - + phi = tdvp(H, -im * tau, phi; nsite=2, time_step=-im * tau, cutoff, normalize=false) + + psi = deepcopy(phi) + phi = tdvp( H, -im * ttotal, phi; nsite=2, time_step=-im * tau, - cutoff = cutoff, + cutoff=cutoff, normalize=false, (observer!)=TDVPObserver(), ) @@ -366,7 +356,7 @@ end -im * ttotal, psi; time_step=-im * tau, - cutoff=cutoff*1e8, ##not clear why this is so much larger than the 2-site TDVP cutoff for comparable bond dimensions? + cutoff=cutoff * 1e8, ##not clear why this is so much larger than the 2-site TDVP cutoff for comparable bond dimensions? cutoff_compress=cutoff, normalize=false, nsite=1, @@ -385,7 +375,7 @@ end #@show norm(En1 - En2) #@show(Sz1) #@show(Sz2) - + @test norm(Sz1 - Sz2) < 1e-3 @test norm(En1 - En2) < 1e-3 end From 01c6891c23b8a9b584a51d896583c3f95a0f74fb Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Wed, 8 Jun 2022 17:02:04 -0400 Subject: [PATCH 13/19] Format src/tdvp_step.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/tdvp_step.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/tdvp_step.jl b/src/tdvp_step.jl index de5ff7b..584eb45 100644 --- a/src/tdvp_step.jl +++ b/src/tdvp_step.jl @@ -30,7 +30,6 @@ function tdvp_step( return psi, PH, info end - isforward(direction::Base.ForwardOrdering) = true isforward(direction::Base.ReverseOrdering) = false isreverse(direction) = !isforward(direction) From de099804420cc71eee9121d8a77aa14993039f79 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 7 Jul 2022 16:54:31 -0400 Subject: [PATCH 14/19] Switch to using nullspace.jl from ITensors --- Project.toml | 2 +- src/ITensorTDVP.jl | 6 +-- src/nullspace.jl | 131 --------------------------------------------- 3 files changed, 4 insertions(+), 135 deletions(-) delete mode 100644 src/nullspace.jl diff --git a/Project.toml b/Project.toml index 155e965..c69a6de 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] -ITensors = "0.3.3" +ITensors = "0.3.16" KrylovKit = "0.5" Observers = "0.0.8" TimerOutputs = "0.5" diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 3a9f538..a4a7a31 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -3,13 +3,13 @@ module ITensorTDVP using ITensors using KrylovKit: exponentiate, eigsolve, svdsolve using Printf -using LinearAlgebra +#using LinearAlgebra using TimerOutputs using Observers using ITensors: - AbstractMPS, @debug_check, @timeit_debug, check_hascommoninds, orthocenter, set_nsite! + AbstractMPS,AbstractProjMPO, @debug_check, @timeit_debug, check_hascommoninds, orthocenter, nullspace,set_nsite! using ITensors.NDTensors using ITensors.NDTensors: eachdiagblock, blockview using ITensors.ITensorNetworkMaps @@ -29,7 +29,7 @@ include("tdvp_generic.jl") include("tdvp.jl") include("dmrg.jl") include("dmrg_x.jl") -include("nullspace.jl") +#include("nullspace.jl") include("subspace_expansion.jl") export tdvp, dmrg_x, to_vec, TimeDependentSum diff --git a/src/nullspace.jl b/src/nullspace.jl deleted file mode 100644 index b7bbd0b..0000000 --- a/src/nullspace.jl +++ /dev/null @@ -1,131 +0,0 @@ - -ITensors.itensor(A::ITensor) = A - -# Insert missing diagonal blocks -function insert_diag_blocks!(T::Tensor) - for b in eachdiagblock(T) - blockT = blockview(T, b) - if isnothing(blockT) - # Block was not found in the list, insert it - insertblock!(T, b) - end - end -end -insert_diag_blocks!(T::ITensor) = insert_diag_blocks!(tensor(T)) - -# Reshape into an order-2 ITensor -matricize(T::ITensor, inds::Index...) = matricize(T, inds) - -function matricize(T::ITensor, inds) - left_inds = commoninds(T, inds) - right_inds = uniqueinds(T, inds) - return matricize(T, left_inds, right_inds) -end - -function matricize(T::ITensor, left_inds, right_inds) - CL = combiner(left_inds; dir=ITensors.Out, tags="CL") - CR = combiner(right_inds; dir=ITensors.In, tags="CR") - M = (T * CL) * CR - return M, CL, CR -end - -function setdims(t::NTuple{N,Pair{QN,Int}}, dims::NTuple{N,Int}) where {N} - return first.(t) .=> dims -end - -# XXX: generalize this function -function _getindex(T::DenseTensor{ElT,N}, I1::Colon, I2::UnitRange{Int64}) where {ElT,N} - A = array(T)[I1, I2] - return tensor(Dense(vec(A)), setdims(inds(T), size(A))) -end - -function getblock(i::Index, n::Integer) - return ITensors.space(i)[n] -end - -# Make `Pair{QN,Int}` act like a regular `dim` -NDTensors.dim(qnv::Pair{QN,Int}) = last(qnv) - -Base.:*(qnv::Pair{QN,Int}, d::ITensors.Arrow) = qn(qnv) * d => dim(qnv) - -function getblock_preserve_qns(T::Tensor, b::Block) - # TODO: make `T[b]` preserve QNs - Tb = T[b] - indsTb = getblock.(inds(T), Tuple(b)) .* dir.(inds(T)) - return ITensors.setinds(Tb, indsTb) -end - -function blocksparsetensor(blocks::Dict{B,TB}) where {B,TB} - b1, Tb1 = first(pairs(blocks)) - N = length(b1) - indstypes = typeof.(inds(Tb1)) - blocktype = eltype(Tb1) - indsT = getindex.(indstypes) - # Determine the indices from the blocks - for (b, Tb) in pairs(blocks) - indsTb = inds(Tb) - for n in 1:N - bn = b[n] - indsTn = indsT[n] - if bn > length(indsTn) - resize!(indsTn, bn) - end - indsTn[bn] = indsTb[n] - end - end - T = BlockSparseTensor(blocktype, indsT) - for (b, Tb) in pairs(blocks) - if !isempty(Tb) - T[b] = Tb - end - end - return T -end - -function _nullspace_hermitian(M::Tensor; atol::Real=0.0) - tol = atol - # Insert any missing diagonal blocks - insert_diag_blocks!(M) - #D, U = eigen(Hermitian(M)) - Dᵢₜ, Uᵢₜ = eigen(itensor(M); ishermitian=true) - D = tensor(Dᵢₜ) - U = tensor(Uᵢₜ) - nullspace_blocks = Dict() - for bU in nzblocks(U) - bM = Block(bU[1], bU[1]) - bD = Block(bU[2], bU[2]) - # Assume sorted from largest to smallest - indstart = sum(d -> abs(d) .> tol, storage(D[bD])) + 1 - Ub = getblock_preserve_qns(U, bU) - indstop = lastindex(Ub, 2) - Nb = _getindex(Ub, :, indstart:indstop) - nullspace_blocks[bU] = Nb - end - return blocksparsetensor(nullspace_blocks) -end - -function LinearAlgebra.nullspace(M::Hermitian{<:Number,<:Tensor}; kwargs...) - return _nullspace_hermitian(parent(M); kwargs...) -end - -function LinearAlgebra.nullspace(::Order{2}, M::ITensor, left_inds, right_inds; kwargs...) - @assert order(M) == 2 - M² = prime(dag(M), right_inds) * M - M² = permute(M², right_inds'..., right_inds...) - M²ₜ = tensor(M²) - Nₜ = nullspace(Hermitian(M²ₜ); kwargs...) - indsN = (Index(ind(Nₜ, 1); dir=ITensors.In), Index(ind(Nₜ, 2); dir=ITensors.In)) - N = dag(itensor(ITensors.setinds(Nₜ, indsN))) - # Make the index match the input index - Ñ = replaceinds(N, (ind(N, 1),) => right_inds) - return Ñ -end - -function LinearAlgebra.nullspace(T::ITensor, is...; kwargs...) - M, CL, CR = matricize(T, is...) - @assert order(M) == 2 - cL = commoninds(M, CL) - cR = commoninds(M, CR) - N₂ = nullspace(Order(2), M, cL, cR; kwargs...) - return N₂ * CR -end From ada52a38d121162b1627407c5ba398598a568154 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 7 Jul 2022 16:56:14 -0400 Subject: [PATCH 15/19] Implement and use non-mutating version of subspace expansion. Use updated interface to ITensors.directsum --- src/subspace_expansion.jl | 16 ++++++++++++++-- src/tdvp_step.jl | 6 ++---- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index b26fecc..64d6344 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -44,6 +44,18 @@ function subspace_expansion_sweep!( return nothing end +function subspace_expansion_sweep( + ψ::MPS, PH::Union{ProjMPO,ProjMPOSum}; maxdim, cutoff, atol=1e-10, kwargs... +) +ψc=copy(ψ) +PHc=copy(PH) +subspace_expansion_sweep!( + ψc, PHc; maxdim, cutoff, atol=1e-10, kwargs... +) +return ψc,PHc +end + + function subspace_expansion!( ψ::MPS, PH, @@ -130,10 +142,10 @@ function subspace_expansion!( ###add expansion direction to current site tensors ALⁿ¹, newl = ITensors.directsum( - ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) + ϕ_1 => uniqueinds(ϕ_1, newL), dag(newL) => uniqueinds(newL, ϕ_1); tags=("Left",) ) ARⁿ², newr = ITensors.directsum( - ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) + ϕ_2 => uniqueinds(ϕ_2, newR), dag(newR) => uniqueinds(newR, ϕ_2); tags=("Right",) ) ###TODO remove assertions regarding expansion not exceeding maxdim ? diff --git a/src/tdvp_step.jl b/src/tdvp_step.jl index 584eb45..4d3e51c 100644 --- a/src/tdvp_step.jl +++ b/src/tdvp_step.jl @@ -9,17 +9,15 @@ function tdvp_step( #psi_evolved = copy(psi) #copy required, otherwise this tdvp function mutates psi in the loop (and subspace expansion if performed) nsite::Int = get(kwargs, :nsite, 2) - expand::Bool = get(kwargs, :expand, nsite == 2 ? false : (hasqns(psi) ? true : false)) + expand::Bool = get(kwargs, :expand, nsite == 2 ? false : (hasqns(psi) ? false : false)) if expand - ###TODO: check that shallow copy leaves psi invariant in the following 'in-place' subspace expansion maxdim = get(kwargs, :maxdim, 100) cutoff = get(kwargs, :cutoff, 1e-11) atol = get(kwargs, :atol, 1e-12) cutoff_trunc = 0.5 * abs(time_step)^2 * cutoff ### cutoff_trunc is min. acceleration of population gain of expansion vectors, thus rescaling with timestep - ITensorTDVP.subspace_expansion_sweep!(psi, PH; maxdim, cutoff=cutoff_trunc, atol=atol) + psi,PH=ITensorTDVP.subspace_expansion_sweep(psi, PH; maxdim, cutoff=cutoff_trunc, atol=atol) end - for substep in 1:length(sub_time_steps) psi, PH, info = tdvp_sweep( orderings[substep], solver, PH, sub_time_steps[substep], psi; current_time, kwargs... From dbd3aa11aa484e1323f0e7ea58e7d2739c7a29e6 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 7 Jul 2022 17:07:45 -0400 Subject: [PATCH 16/19] Format. --- src/ITensorTDVP.jl | 9 ++++++++- src/subspace_expansion.jl | 11 ++++------- src/tdvp_step.jl | 4 +++- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index a4a7a31..692e11d 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -9,7 +9,14 @@ using TimerOutputs using Observers using ITensors: - AbstractMPS,AbstractProjMPO, @debug_check, @timeit_debug, check_hascommoninds, orthocenter, nullspace,set_nsite! + AbstractMPS, + AbstractProjMPO, + @debug_check, + @timeit_debug, + check_hascommoninds, + orthocenter, + nullspace, + set_nsite! using ITensors.NDTensors using ITensors.NDTensors: eachdiagblock, blockview using ITensors.ITensorNetworkMaps diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 64d6344..37abcb6 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -47,15 +47,12 @@ end function subspace_expansion_sweep( ψ::MPS, PH::Union{ProjMPO,ProjMPOSum}; maxdim, cutoff, atol=1e-10, kwargs... ) -ψc=copy(ψ) -PHc=copy(PH) -subspace_expansion_sweep!( - ψc, PHc; maxdim, cutoff, atol=1e-10, kwargs... -) -return ψc,PHc + ψc = copy(ψ) + PHc = copy(PH) + subspace_expansion_sweep!(ψc, PHc; maxdim, cutoff, atol=1e-10, kwargs...) + return ψc, PHc end - function subspace_expansion!( ψ::MPS, PH, diff --git a/src/tdvp_step.jl b/src/tdvp_step.jl index 4d3e51c..3c3068b 100644 --- a/src/tdvp_step.jl +++ b/src/tdvp_step.jl @@ -16,7 +16,9 @@ function tdvp_step( cutoff = get(kwargs, :cutoff, 1e-11) atol = get(kwargs, :atol, 1e-12) cutoff_trunc = 0.5 * abs(time_step)^2 * cutoff ### cutoff_trunc is min. acceleration of population gain of expansion vectors, thus rescaling with timestep - psi,PH=ITensorTDVP.subspace_expansion_sweep(psi, PH; maxdim, cutoff=cutoff_trunc, atol=atol) + psi, PH = ITensorTDVP.subspace_expansion_sweep( + psi, PH; maxdim, cutoff=cutoff_trunc, atol=atol + ) end for substep in 1:length(sub_time_steps) psi, PH, info = tdvp_sweep( From b4d77a29ce313772e2829f6435d5e6fbacc7d84b Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 8 Jul 2022 14:31:11 -0400 Subject: [PATCH 17/19] Formatting and cleaning up subspace expansion. --- src/subspace_expansion.jl | 58 +++++++++++++++++++++------------------ src/tdvp_step.jl | 9 +++--- 2 files changed, 36 insertions(+), 31 deletions(-) diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl index 37abcb6..8019069 100644 --- a/src/subspace_expansion.jl +++ b/src/subspace_expansion.jl @@ -67,32 +67,32 @@ function subspace_expansion!( """ expands in nullspace of site-tensors b """ - #atol refers to the tolerance in nullspace determination (for finite MPS can probably be set rather small) - #cutoff refers to largest singular value of gradient (acceleration of population gain of expansion vectors) to keep - #this should be related to the cutoff in two-site TDVP: \Delta_rho_i = 0.5 * lambda_y * tau **2 - #note that in the initial SVD there is another cutoff parameter, that we set to roughly machine precision for the time being - #(allows to move bond-dimension between different partial QN sectors, if no such truncation occurs distribution of bond-dimensions - #between different QNs locally is static once bond dimension saturates maxdim.) - - ##this should only work for the case where rlim-llim > 1 - ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) + # atol refers to the tolerance in nullspace determination (for finite MPS can probably be set rather small) + # cutoff refers to largest singular value of gradient (acceleration of population gain of expansion vectors) to keep + # this should be related to the cutoff in two-site TDVP: \Delta_rho_i = 0.5 * lambda_y * tau **2 + # note that in the initial SVD there is another cutoff parameter `cutoff_compress`, that we set to roughly machine precision for the time being + # (allows to move bond-dimension between different partial QN sectors, if no such truncation occurs distribution of bond-dimensions + # between different QNs locally is static once bond dimension saturates maxdim.) + llim, rlim = lims n1, n2 = b PH.nsite = 2 old_linkdim = dim(commonind(ψ[n1], ψ[n2])) cutoff_compress = get(kwargs, :cutoff_compress, 1e-12) - ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly + # move orthogonality center to bond + # and check whether there are vanishing contributions to the wavefunctions + # truncate accordingly + # at the cost of unitarity but introducing flexibility of redistributing bond dimensions among QN sectors) if llim == n1 @assert rlim == n2 + 1 U, S, V = svd( ψ[n2], uniqueinds(ψ[n2], ψ[n1]); maxdim=maxdim, cutoff=cutoff_compress, kwargs... - ) ##lookup svd interface again + ) ϕ_1 = ψ[n1] * V ϕ_2 = U old_linkdim = dim(commonind(U, S)) bondtensor = S - elseif rlim == n2 @assert llim == n1 - 1 U, S, V = svd( @@ -104,40 +104,45 @@ function subspace_expansion!( bondtensor = S end - ###don't expand if we are already at maxdim + # don't expand if we are already at maxdim if old_linkdim >= maxdim println("not expanding") return nothing end position!(PH, ψ, min(n1, n2)) - #orthogonalize(ψ,n1) + # orthogonalize(ψ,n1) linkind_l = commonind(ϕ_1, bondtensor) linkind_r = commonind(ϕ_2, bondtensor) + # compute nullspace to the left and right NL = nullspace(ϕ_1, linkind_l; atol=atol) NR = nullspace(ϕ_2, linkind_r; atol=atol) - ###NOTE: This will fail for rank-1 trees with physical DOFs on leafs - ###NOTE: one-sided subspace expansion seems to not work well at least for trees according to Lachlan Lindoy + # if nullspace is empty (happen's for product states with QNs) if norm(NL) == 0.0 || norm(NR) == 0.0 return bondtensor end - ###form 2site wavefunction + # form 2-site wavefunction ϕ = ϕ_1 * bondtensor * ϕ_2 - ###get subspace expansion + # get subspace expansion newL, S, newR, success = _subspace_expand_core( ϕ, PH, NL, NR, ; maxdim=maxdim - old_linkdim, cutoff, kwargs... ) - @assert success - #@show expanS + + # ToDo: Handle failed subspace expansion properly. Maybe pass an argument whether to raise an error if subspace expansion fails + # @assert success + if success == false + println( + "Subspace expansion not successful. This may indicate that 2-site TDVP also fails for the given state and Hamiltonian.", + ) return nothing end - ###add expansion direction to current site tensors + # expand current site tensors ALⁿ¹, newl = ITensors.directsum( ϕ_1 => uniqueinds(ϕ_1, newL), dag(newL) => uniqueinds(newL, ϕ_1); tags=("Left",) ) @@ -145,7 +150,9 @@ function subspace_expansion!( ϕ_2 => uniqueinds(ϕ_2, newR), dag(newR) => uniqueinds(newR, ϕ_2); tags=("Right",) ) - ###TODO remove assertions regarding expansion not exceeding maxdim ? + # Some checks + # ToDo: clean up + # ToDo: remove assertions regarding expansion not exceeding maxdim ? @assert (dim(commonind(newL, S)) + old_linkdim) <= maxdim @assert dim(commonind(newL, S)) == dim(commonind(newR, S)) @assert(dim(uniqueind(ϕ_1, newL)) + dim(uniqueind(newL, ϕ_1)) == dim(newl)) @@ -155,7 +162,7 @@ function subspace_expansion!( @assert dim(newl) <= maxdim @assert dim(newr) <= maxdim - ###zero-pad bond-tensor (the orthogonality center) + # zero-pad bond-tensor (the orthogonality center) C = ITensor(dag(newl)..., dag(newr)...) ψC = bondtensor for I in eachindex(ψC) @@ -165,13 +172,12 @@ function subspace_expansion!( end end - ###move orthogonality center back to site (should restore input orthogonality limits) + # move orthogonality center back to site (should restore input orthogonality limits) if rlim == n2 ψ[n2] = ARⁿ² elseif rlim > n2 ψ[n2] = noprime(ARⁿ² * C) end - if llim == n1 ψ[n1] = ALⁿ¹ elseif llim < n1 @@ -197,6 +203,7 @@ function _subspace_expand_core(ϕ::ITensor, env, NL, NR; maxdim, cutoff, kwargs. if norm(ϕH) == 0.0 return false, false, false, false end + U, S, V = svd( ϕH, commoninds(ϕH, NL); @@ -206,7 +213,6 @@ function _subspace_expand_core(ϕ::ITensor, env, NL, NR; maxdim, cutoff, kwargs. use_absolute_cutoff=true, kwargs..., ) - @assert dim(commonind(U, S)) <= maxdim NL *= dag(U) diff --git a/src/tdvp_step.jl b/src/tdvp_step.jl index 3c3068b..5af5dc3 100644 --- a/src/tdvp_step.jl +++ b/src/tdvp_step.jl @@ -5,17 +5,16 @@ function tdvp_step( sub_time_steps = ITensorTDVP.sub_time_steps(order) sub_time_steps *= time_step global info - - #psi_evolved = copy(psi) #copy required, otherwise this tdvp function mutates psi in the loop (and subspace expansion if performed) - nsite::Int = get(kwargs, :nsite, 2) expand::Bool = get(kwargs, :expand, nsite == 2 ? false : (hasqns(psi) ? false : false)) if expand + # do subspace expansion maxdim = get(kwargs, :maxdim, 100) cutoff = get(kwargs, :cutoff, 1e-11) atol = get(kwargs, :atol, 1e-12) - cutoff_trunc = 0.5 * abs(time_step)^2 * cutoff ### cutoff_trunc is min. acceleration of population gain of expansion vectors, thus rescaling with timestep + # cutoff_trunc is should scale with acceleration of population gain of expansion vectors, thus rescaling with timestep + cutoff_trunc = 0.5 * abs(time_step)^2 * cutoff psi, PH = ITensorTDVP.subspace_expansion_sweep( psi, PH; maxdim, cutoff=cutoff_trunc, atol=atol ) @@ -26,7 +25,7 @@ function tdvp_step( ) current_time += sub_time_steps[substep] end - #@show maxlinkdim(psi) + return psi, PH, info end From 4fdbe91f3f97e93740ce946b21723cc7147a5cb8 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 8 Jul 2022 14:31:47 -0400 Subject: [PATCH 18/19] Remove LinearAlgebra from dependencies. --- src/ITensorTDVP.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 692e11d..e4d2f68 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -3,7 +3,6 @@ module ITensorTDVP using ITensors using KrylovKit: exponentiate, eigsolve, svdsolve using Printf -#using LinearAlgebra using TimerOutputs using Observers From 45641a4b8cf8c68df6f2a72cd3833beac5009618 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Mon, 11 Jul 2022 08:26:49 -0400 Subject: [PATCH 19/19] remove Krylov subspace expansion --- src/subspace_expansion_krylov.jl | 201 ------------------------------- 1 file changed, 201 deletions(-) delete mode 100644 src/subspace_expansion_krylov.jl diff --git a/src/subspace_expansion_krylov.jl b/src/subspace_expansion_krylov.jl deleted file mode 100644 index bb2bb28..0000000 --- a/src/subspace_expansion_krylov.jl +++ /dev/null @@ -1,201 +0,0 @@ -function replaceind_indval(IV::Tuple, iĩ::Pair) - i, ĩ = iĩ - return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV)) -end - -# atol controls the tolerance cutoff for determining which eigenvectors are in the null -# space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep -# the eigenvectors corresponding to eigenvalues of approximately 1. - -###the most general implementation may be passing MPS, lims, and bondtensor (defaulting to the Id-matrix) -###then we can tensorsum both left and right tensor, and just apply the bondtensor to restore input gauge -###ideally we would also be able to apply tensorsum to the bond tensor instead of that loop below -""" -expand subspace (2-site like) in a sweep -""" -function subspace_expansion_krylov_sweep!( - ψ::MPS, PH::Union{ProjMPO,ProjMPOSum}; maxdim, cutoff, atol=1e-2, kwargs... -) - N = length(ψ) - if !isortho(ψ) || orthocenter(ψ) != 1 - orthogonalize!(ψ, 1) - end - PH.nsite = 2 - nsite = 2 - position!(PH, ψ, 1) - for (b, ha) in sweepnext(N; ncenter=2) - - ##TODO: figure out whether these calls should be here or inside subspace expansion, currently we do both? - orthogonalize!(ψ, b) - position!(PH, ψ, b) - - if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1) - b1 = (ha == 1 ? b + 1 : b) - Δ = (ha == 1 ? +1 : -1) - inds = (ha == 1 ? (b, b + Δ) : (b + Δ, b)) - subspace_expansion_krylov!( - ψ, PH, (ψ.llim, ψ.rlim), inds; maxdim, cutoff, atol=atol, kwargs... - ) - end - end - return nothing -end - -function subspace_expansion_krylov!( - ψ::MPS, - PH, - lims::Tuple{Int,Int}, - b::Tuple{Int,Int}; - bondtensor=nothing, - maxdim, - cutoff, - atol=1e-2, - kwargs..., -) - ##this should only work for the case where rlim-llim > 1 - ##not a valid MPS otherwise anyway (since bond matrix not part of MPS unless so defined like in VidalMPS struct) - llim, rlim = lims - n1, n2 = b - PH.nsite = 2 - old_linkdim = dim(commonind(ψ[n1], ψ[n2])) - - ###move orthogonality center to bond, check whether there are vanishing contributions to the wavefunctions and truncate accordingly - ###the cutoff should be scaled with timestep, otherwise one runs into problems with non-monotonic error behaviour like in TEBD approaches - if llim == n1 - @assert rlim == n2 + 1 - U, S, V = svd(ψ[n2], uniqueinds(ψ[n2], ψ[n1]); maxdim=maxdim, cutoff=1e-17, kwargs...) ##lookup svd interface again - ϕ_1 = ψ[n1] * V - ϕ_2 = U - old_linkdim = dim(commonind(U, S)) - bondtensor = S - - elseif rlim == n2 - @assert llim == n1 - 1 - U, S, V = svd(ψ[n1], uniqueinds(ψ[n1], ψ[n2]); maxdim=maxdim, cutoff=1e-17, kwargs...) - ϕ_1 = U - ϕ_2 = ψ[n2] * V - old_linkdim = dim(commonind(U, S)) - bondtensor = S - end - - ###don't expand if we are already at maxdim - if old_linkdim >= maxdim - println("not expanding") - return nothing - end - position!(PH, ψ, min(n1, n2)) - - #orthogonalize(ψ,n1) - linkind_l = commonind(ϕ_1, bondtensor) - linkind_r = commonind(ϕ_2, bondtensor) - - NL = nullspace(ϕ_1, linkind_l; atol=atol) - NR = nullspace(ϕ_2, linkind_r; atol=atol) - - ###NOTE: This will fail for rank-1 trees with physical DOFs on leafs - ###NOTE: one-sided subspace expansion seems to not work well at least for trees according to Lachlan Lindoy - if norm(NL) == 0.0 || norm(NR) == 0.0 - return bondtensor - end - - ###form 2site wavefunction - ϕ = [ϕ_1, bondtensor, ϕ_2] - - ###get subspace expansion - newL, S, newR, success = _subspace_expand_core_krylov( - ϕ, PH, NL, NR, ; maxdim=maxdim - old_linkdim, cutoff, kwargs... - ) - if success == false - return nothing - end - - ###add expansion direction to current site tensors - ALⁿ¹, newl = ITensors.directsum( - ϕ_1, dag(newL), uniqueinds(ϕ_1, newL), uniqueinds(newL, ϕ_1); tags=("Left",) - ) - ARⁿ², newr = ITensors.directsum( - ϕ_2, dag(newR), uniqueinds(ϕ_2, newR), uniqueinds(newR, ϕ_2); tags=("Right",) - ) - - ###TODO remove assertions regarding expansion not exceeding maxdim - @assert (dim(commonind(newL, S)) + old_linkdim) <= maxdim - @assert dim(commonind(newL, S)) == dim(commonind(newR, S)) - @assert(dim(uniqueind(ϕ_1, newL)) + dim(uniqueind(newL, ϕ_1)) == dim(newl)) - @assert(dim(uniqueind(ϕ_2, newR)) + dim(uniqueind(newR, ϕ_2)) == dim(newr)) - @assert (old_linkdim + dim(commonind(newL, S))) <= maxdim - @assert (old_linkdim + dim(commonind(newR, S))) <= maxdim - @assert dim(newl) <= maxdim - @assert dim(newr) <= maxdim - - ###zero-pad bond-tensor (the orthogonality center) - C = ITensor(dag(newl)..., dag(newr)...) - ψC = bondtensor - ### FIXME: the permute below fails, maybe because this already the layout of bondtensor --- in any case it shouldn't fail? - #ψC = permute(bondtensor, linkind_l, linkind_r) - for I in eachindex(ψC) - v = ψC[I] - if !iszero(v) - C[I] = ψC[I] - end - end - - ###move orthogonality center back to site (should restore input orthogonality limits) - if rlim == n2 - ψ[n2] = ARⁿ² - elseif rlim > n2 - ψ[n2] = noprime(ARⁿ² * C) - end - - if llim == n1 - ψ[n1] = ALⁿ¹ - elseif llim < n1 - ψ[n1] = noprime(ALⁿ¹ * C) - end - - return nothing -end - -function _subspace_expand_core_krylov( - centerwf::Vector{ITensor}, env, NL, NR; maxdim, cutoff, kwargs... -) - maxdim = min(maxdim, 15) - envL = lproj(env) - envR = rproj(env) - envsiteMPOs = env.H[ITensors.site_range(env)] - ##NOTE: requires n1+1=n2, otherwise assignment of left and right is of - L = envL * centerwf[1] - L *= envsiteMPOs[1] - R = envR * centerwf[3] - R *= envsiteMPOs[2] - C = centerwf[2] - R *= C ##contract C into one of them because otherwise the application is more costly? - R = noprime(R) - L = noprime(L) - outinds = uniqueinds(NL, NL) - ininds = uniqueinds(NR, NR) - B2 = ITensorNetworkMap([NR, R, L, NL], ininds, outinds) - B2dag = adjoint(B2) - trial = randomITensor(eltype(L), uniqueinds(NL, L)) - trialadj = randomITensor(eltype(R), uniqueinds(NR, R)) #columnspace of B2, i.e. NL - #trial=noprime(B2(noprime(B2dag(trial)))) - #vals, lvecs, rvecs, info = svdsolve(trial, maxdim) do (x, flag) - # if flag - # y = B2dag * copy(x)# y = compute action of adjoint map on x - # else - # y = B2 * copy(x)# y = compute action of linear map on x - # end - # return y - #end - vals, lvecs, rvecs, info = svdsolve( - (x -> noprime(B2 * x), y -> noprime(B2dag * y)), trial, maxdim - ) ###seems to work - - #TO DO construct U,S,V objects, using only the vals > cutoff, and at most maxdim - # - @show vals - @show uniqueinds(NL, L) - @show uniqueinds(NR, R) - @show inds(lvecs[1]) - @show inds(rvecs[1]) - return println("success!!") -end