From 0d504a1a8feca16cd2b0fcee77913598e87c0e15 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Mon, 1 Apr 2024 16:47:55 -0400 Subject: [PATCH 01/19] Initial commit, checking out files from other branch. --- .../local_solvers/expand_and_exponentiate.jl | 94 +++++ .../subspace_expansions/linalg/rsvd_aux.jl | 251 +++++++++++ .../subspace_expansions/linalg/rsvd_linalg.jl | 259 ++++++++++++ .../linalg/standard_svd.jl | 12 + .../subspace_expansions/subspace_expansion.jl | 398 ++++++++++++++++++ 5 files changed, 1014 insertions(+) create mode 100644 src/solvers/local_solvers/expand_and_exponentiate.jl create mode 100644 src/solvers/subspace_expansions/linalg/rsvd_aux.jl create mode 100644 src/solvers/subspace_expansions/linalg/rsvd_linalg.jl create mode 100644 src/solvers/subspace_expansions/linalg/standard_svd.jl create mode 100644 src/solvers/subspace_expansions/subspace_expansion.jl diff --git a/src/solvers/local_solvers/expand_and_exponentiate.jl b/src/solvers/local_solvers/expand_and_exponentiate.jl new file mode 100644 index 00000000..6df90efe --- /dev/null +++ b/src/solvers/local_solvers/expand_and_exponentiate.jl @@ -0,0 +1,94 @@ +function local_expand_and_exponentiate_updater( + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs, +) + # expansion + expand_kwargs = updater_kwargs.expand_kwargs # defaults are set in the expansion_updater + expanded_init, _ = two_site_expansion_updater( + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs=expand_kwargs, + ) + + # exponentiate + # ToDo: also call exponentiate_updater, instead of reimplementing it here + # + default_exponentiate_kwargs = (; + krylovdim=30, + maxiter=100, + verbosity=0, + tol=1E-8, + ishermitian=true, + issymmetric=true, + eager=true, + ) + exponentiate_kwargs = updater_kwargs.exponentiate_kwargs + exponentiate_kwargs = merge(default_exponentiate_kwargs, exponentiate_kwargs) #last collection has precedence + result, exp_info = exponentiate( + projected_operator![], region_kwargs.time_step, expanded_init; exponentiate_kwargs... + ) + + # + #ToDo: return truncation error and append to info + # truncate + # ToDo: refactor + # ToDo: either remove do_truncate, or make accessible as kwarg + # ToDo: think about more elaborate ways to trigger truncation, i.e. truncate only once maxdim is hit? + # ToDo: add truncation info + do_truncate = true + # simply return if we are not truncating + !do_truncate && return result, (; info=exp_info) + # check for cases where we don't want to try truncating + region = first(sweep_plan[which_region_update]) + typeof(region) <: NamedEdge && return result, (; info=exp_info) + (; maxdim, cutoff) = region_kwargs + region = only(region) + next_region = if which_region_update == length(sweep_plan) + nothing + else + first(sweep_plan[which_region_update + 1]) + end + previous_region = + which_region_update == 1 ? nothing : first(sweep_plan[which_region_update - 1]) + isnothing(next_region) && return result, (; info=exp_info) + !(typeof(next_region) <: NamedEdge) && return result, (; info=exp_info) + !(region == src(next_region) || region == dst(next_region)) && + return result, (; info=exp_info) + # actually truncate + left_inds = uniqueinds(state![], next_region) + U, S, V = svd( + result, + left_inds; + lefttags=tags(state![], next_region), + righttags=tags(state![], next_region), + maxdim, + cutoff, + ) + next_vertex = src(next_region) == region ? dst(next_region) : src(next_region) + state = copy(state![]) + #@show inds(V) + state[next_vertex] = state[next_vertex] * V + _nsites = (region isa AbstractEdge) ? 0 : length(region) #should be 1 + #@show _nsites + PH = copy(projected_operator![]) + PH = set_nsite(PH, 2) + PH = position(PH, state, [region, next_vertex]) + PH = set_nsite(PH, _nsites) + PH = position(PH, state, first(sweep_plan[which_region_update])) + state![] = state + projected_operator![] = PH + return U * S, (; info=exp_info) +end diff --git a/src/solvers/subspace_expansions/linalg/rsvd_aux.jl b/src/solvers/subspace_expansions/linalg/rsvd_aux.jl new file mode 100644 index 00000000..cedc11c1 --- /dev/null +++ b/src/solvers/subspace_expansions/linalg/rsvd_aux.jl @@ -0,0 +1,251 @@ +function get_column_space(A::Vector{<:ITensor}, lc::Index,rc::Index) + #gets non-zero blocks in rc by sticking in lc and contracting through + viable_sectors=Vector{Pair{QN,Int64}} + for s in space(lc) + qn = first(s) + trial=randomITensor(flux(dag(qn)),lc) + adtrial=foldl(A;init=trial) + nnzblocks(adtrial)==0 && continue + thesector=only(space(only(inds(adtrial)))) + push!(viable_sectors, thesector) + end + return viable_sectors +end + + +function build_guess_matrix( + eltype::Type{<:Number}, ind, sectors::Union{Nothing,Vector{Pair{QN,Int64}}}, n::Int, p::Int + ) + if hasqns(ind) + aux_spaces = Pair{QN,Int64}[] + for s in sectors + thedim = last(s) + qn = first(s) + en = min(n + p, thedim) + push!(aux_spaces, Pair(qn, en)) + end + aux_ind = Index(aux_spaces; dir=dir(ind)) + try + M = randomITensor(eltype, dag(ind), aux_ind) #defaults to zero flux + # @show theflux, flux(M) + catch e + @show e + @show aux_ind + @show ind + error("stopping here something is wrong") + end + @assert nnzblocks(M) != 0 + else + thedim = dim(ind) + en = min(n + p, thedim) + #ep = max(0,en-n) + aux_ind = Index(en) + M = randomITensor(eltype, ind, aux_ind) + end + #@show M + return M +end +qns(t::ITensor) = qns(collect(eachnzblock(t)), t) +function qns(bs::Vector{Block{n}}, t::ITensor) where {n} + return [qns(b, t) for b in bs] +end + +function qns(b::Block{n}, t::ITensor) where {n} + theqns = QN[] + for i in 1:order(t) + theb = getindex(b, i) + theind = inds(t)[i] + push!(theqns, ITensors.qn(space(theind), Int(theb))) + end + return theqns +end + +function build_guess_matrix( + eltype::Type{<:Number}, ind, sectors::Union{Nothing,Vector{Pair{QN,Int64}}}, ndict::Dict; theflux=nothing, auxdir=dir(ind) +) + if hasqns(ind) + #translate_qns = get_qn_dict(ind, theflux; auxdir) + aux_spaces = Pair{QN,Int64}[] + #@show first.(space(ind)) + for s in sectors + thedim = last(s) + qn = first(s) + en = min(ndict[qn], thedim) + push!(aux_spaces, Pair(qn, en)) + end + aux_ind = Index(aux_spaces; dir=dir(ind)) + #M=randomITensor(eltype,-theflux,dag(ind),aux_ind) + #@show aux_ind + try + M = randomITensor(eltype, dag(ind), aux_ind) + catch e + @show e + @show aux_ind + @show ind + error("stopping here something is wrong") + end + @assert nnzblocks(M) != 0 + else + thedim = dim(ind) + en = min(ndict[space(ind)], thedim) + aux_ind = Index(en) + M = randomITensor(eltype, ind, aux_ind) + end + return M +end + +function init_guess_sizes(cind, sectors::Union{Nothing,Vector{Pair{QN,Int64}}}, n::Int, rule; theflux=nothing, auxdir=dir(cind)) + if hasqns(cind) + ndict = Dict{QN,Int64}() + for s in sectors + thedim = last(s) + qn = first(s) + ndict[qn] = min(rule(n), thedim) + end + else + @assert sectors==nothing + thedim = dim(cind) + ndict = Dict{Int64,Int64}() + ndict[thedim] = min(thedim, rule(n)) + end + return ndict +end + +function increment_guess_sizes(ndict::Dict{QN,Int64}, n_inc::Int, rule) + for key in keys(ndict) + #thedim=last(key) + ndict[key] = ndict[key] + rule(n_inc) + end + return ndict +end + +function increment_guess_sizes(ndict::Dict{Int64,Int64}, n_inc::Int, rule) + for key in keys(ndict) + #thedim=key + ndict[key] = ndict[key] + rule(n_inc) + end + return ndict +end + +function approx_the_same(o, n; abs_eps=1e-12, rel_eps=1e-6) + absdev = abs.(o .- n) + reldev = abs.((o .- n) ./ n) + abs_conv = absdev .< abs_eps + rel_conv = reldev .< rel_eps + return all(abs_conv .|| rel_conv) +end + +function is_converged_block(o, n; svd_kwargs...) + maxdim = get(svd_kwargs, :maxdim, Inf) + if length(o) != length(n) + return false + else + r = min(maxdim, length(o)) + #ToDo: Pass kwargs? + return approx_the_same(o[1:r], n[1:r]) + end +end + +function is_converged!(ndict, old_fact, new_fact; n_inc=1, has_qns=true, svd_kwargs...) + oS = old_fact.S + nS = new_fact.S + theflux = flux(nS) + oldflux = flux(oS) + if has_qns + if oldflux == nothing || theflux == nothing + if norm(oS) == 0.0 && norm(nS) == 0.0 + return true + else + return false + end + else + try + @assert theflux == flux(oS) + catch e + @show e + @show theflux + @show oldflux + error("Somehow the fluxes are not matching here! Exiting") + end + end + else + ###not entirely sure if this is legal for empty factorization + if norm(oS) == 0.0 + if norm(nS) == 0.0 + return true + else + return false + end + end + end + + maxdim = get(svd_kwargs, :maxdim, Inf) + os = sort(storage(oS); rev=true) + ns = sort(storage(nS); rev=true) + if length(os) >= maxdim && length(ns) >= maxdim + conv_global = approx_the_same(os[1:maxdim], ns[1:maxdim]) + elseif length(os) != length(ns) + conv_global = false + else + r = length(ns) + conv_global = approx_the_same(os[1:r], ns[1:r]) + end + if !hasqns(oS) + if conv_global == false + #ndict has only one key + ndict[first(keys(ndict))] *= 2 + end + return conv_global + end + conv_bool_total = true + ##a lot of this would be more convenient with ITensor internal_inds_space + ##ToDo: refactor, but not entirely trivial because functionality is implemented on the level of QNIndex, not a set of QNIndices + ##e.g. it is cumbersome to query the collection of QNs associated with a Block{n} of an ITensor with n>1 + soS = space(inds(oS)[1]) + snS = space(inds(nS)[1]) + qns = union(ITensors.qn.(soS), ITensors.qn.(snS)) + + oblocks = eachnzblock(oS) + oblockdict = Int.(getindex.(oblocks, 1)) + oqnindtoblock = Dict(collect(values(oblockdict)) .=> collect(keys(oblockdict))) + + nblocks = eachnzblock(nS) + nblockdict = Int.(getindex.(nblocks, 1)) + nqnindtoblock = Dict(collect(values(nblockdict)) .=> collect(keys(nblockdict))) + + for qn in qns + if qn in ITensors.qn.(snS) && qn in ITensors.qn.(soS) + oqnind = findfirst((first.(soS)) .== [qn]) + nqnind = findfirst((first.(snS)) .== (qn,)) + oblock = oqnindtoblock[oqnind] + nblock = nqnindtoblock[nqnind] + #oblock=ITensors.block(first,inds(oS)[1],qn) + #nblock=ITensors.block(first,inds(nS)[1],qn) + + #make sure blocks are the same QN when we compare them + #@assert first(soS[oqnind])==first(snS[nqnind])# + ovals = diag(oS[oblock]) + nvals = diag(nS[nblock]) + conv_bool = is_converged_block(collect(ovals), collect(nvals); svd_kwargs...) + else + conv_bool = false + end + if conv_bool == false + ndict[qn] *= 2 + end + conv_bool_total *= conv_bool + end + if conv_bool_total == true + @assert conv_global == true + else + if conv_global == true + println( + "Subspace expansion, rand. svd: singular vals converged globally, but may not be optimal, doing another iteration", + ) + end + end + return conv_bool_total::Bool +end + + + diff --git a/src/solvers/subspace_expansions/linalg/rsvd_linalg.jl b/src/solvers/subspace_expansions/linalg/rsvd_linalg.jl new file mode 100644 index 00000000..d53d2030 --- /dev/null +++ b/src/solvers/subspace_expansions/linalg/rsvd_linalg.jl @@ -0,0 +1,259 @@ +function rsvd_iterative( + T, + A::ITensors.ITensorNetworkMaps.ITensorNetworkMap, + linds::Vector{<:Index}; + theflux=nothing, + svd_kwargs..., + ) + + #translate from in/out to l/r logic + ininds = ITensors.ITensorNetworkMaps.input_inds(A) + outinds = ITensors.ITensorNetworkMaps.output_inds(A) + @assert linds == ininds ##FIXME: this is specific to the way we wrote the subspace expansion, should be fixed in another iteration + rinds = outinds + + CL = combiner(linds...) + CR = combiner(rinds...) + cL = uniqueind(inds(CL), linds) + cR = uniqueind(inds(CR), rinds) + + l = CL * A.itensors[1] + r = A.itensors[end] * CR + + if length(A.itensors) !== 2 + AC = ITensors.ITensorNetworkMaps.ITensorNetworkMap( + [l, A.itensors[2:(length(A.itensors) - 1)]..., r], + commoninds(l, CL), + commoninds(r, CR), + ) + else + AC = ITensors.ITensorNetworkMaps.ITensorNetworkMap( + [l, r], commoninds(l, CL), commoninds(r, CR) + ) + end + ###this initializer part is still a bit ugly + n_init = 1 + p_rule(n) = 2 * n + ndict2 = init_guess_sizes(cR, n_init, p_rule; theflux=theflux) + ndict = init_guess_sizes(cL, n_init, p_rule; theflux=theflux) + if hasqns(ininds) + ndict = merge(ndict, ndict2) + else + ndict = ndict2 + end + M = build_guess_matrix(T, cR, ndict; theflux=theflux) + fact, Q = rsvd_core(AC, M; svd_kwargs...) + n_inc = 2 + ndict = increment_guess_sizes(ndict, n_inc, p_rule) + new_fact = deepcopy(fact) + while true + M = build_guess_matrix(T, cR, ndict; theflux=theflux) + new_fact, Q = rsvd_core(AC, M; svd_kwargs...) + if is_converged!(ndict, fact, new_fact; n_inc, has_qns=hasqns(ininds), svd_kwargs...) + break + else + fact = new_fact + end + end + vals = diag(array(new_fact.S)) + (length(vals) == 1 && vals[1]^2 ≤ get(svd_kwargs, :cutoff, 0.0)) && + return nothing, nothing, nothing + return dag(CL) * Q * new_fact.U, new_fact.S, new_fact.V * dag(CR) + end + + function rsvd_iterative(A::ITensor, linds::Vector{<:Index}; svd_kwargs...) + rinds = uniqueinds(A, linds) + CL = combiner(linds) + CR = combiner(rinds) + AC = CL * A * CR + cL = combinedind(CL) + cR = combinedind(CR) + if inds(AC) != (cL, cR) + AC = permute(AC, cL, cR) + end + n_init = 1 + p_rule(n) = 2 * n + iszero(norm(AC)) && return nothing, nothing, nothing + #@show flux(AC) + ndict2 = init_guess_sizes(cR, n_init, p_rule; theflux=hasqns(AC) ? flux(AC) : nothing) + ndict = init_guess_sizes(cL, n_init, p_rule; theflux=hasqns(AC) ? flux(AC) : nothing) + ndict = merge(ndict, ndict2) + M = build_guess_matrix(eltype(AC), cR, ndict; theflux=hasqns(AC) ? flux(AC) : nothing) + fact, Q = rsvd_core(AC, M; svd_kwargs...) + n_inc = 1 + ndict = increment_guess_sizes(ndict, n_inc, p_rule) + new_fact = deepcopy(fact) + while true + M = build_guess_matrix(eltype(AC), cR, ndict; theflux=hasqns(AC) ? flux(AC) : nothing) + new_fact, Q = rsvd_core(AC, M; svd_kwargs...) + if is_converged!(ndict, fact, new_fact; n_inc, has_qns=hasqns(AC), svd_kwargs...) + break + else + fact = new_fact + end + end + vals = diag(array(new_fact.S)) + (length(vals) == 1 && vals[1]^2 ≤ get(svd_kwargs, :cutoff, 0.0)) && + return nothing, nothing, nothing + #@show flux(dag(CL)*Q*new_fact.U) + #@show flux(new_fact.S)\ + + @assert flux(new_fact.S) == flux(AC) + return dag(CL) * Q * new_fact.U, new_fact.S, new_fact.V * dag(CR) + #ToDo?: handle non-QN case separately because there it is advisable to start with n_init closer to target maxdim_expand + ##not really an issue anymore since we do *2 increase, so only log number of calls + end + + function rsvd_iterative(A::Vector{ITensor}, linds::Vector{<:Index}; svd_kwargs...) + if length(A)==1 + rinds = uniqueinds(only(A), linds) + else + rinds = uniqueinds(A[end],unioninds(A[1:end-1]...)) + end + CL = combiner(linds) + CR = combiner(rinds) + AC=copy(A) + AC[1] = CL*first(AC) + AC[end] = last(AC)*CR + cL = combinedind(CL) + cR = combinedind(CR) + theflux = any(hasqns.(AC)) ? reduce(+,flux.(AC)) : nothing + #@show theflux + n_init = 1 + p_rule(n) = 2 * n + iszero(norm(AC)) && return nothing, nothing, nothing + #@show flux(AC) + ndict2 = init_guess_sizes(cR, n_init, p_rule; theflux) + #@assert isempty(setdiff(keys(Dict(space(cR))), keys(ndict2))) + ndict = init_guess_sizes(cL, n_init, p_rule; theflux) + #FIXME: merging is not the right thing to do, but is a workaround due to the way is_converged! is implemented + ndict = merge(ndict, ndict2) + #@show keys(ndict), keys(Dict(space(cR))) + #@assert isempty(setdiff(keys(Dict(space(cR))), keys(ndict))) + #@assert isempty(setdiff(keys(ndict),keys(Dict(space(cR))))) + + M = build_guess_matrix(eltype(first(AC)), cR, ndict; theflux) + fact, Q = rsvd_core(AC, M; svd_kwargs...) + n_inc = 1 + ndict = increment_guess_sizes(ndict, n_inc, p_rule) + new_fact = deepcopy(fact) + while true + M = build_guess_matrix(eltype(first(AC)), cR, ndict; theflux) + new_fact, Q = rsvd_core(AC, M; svd_kwargs...) + isnothing(Q) && return nothing,nothing,nothing + if is_converged!(ndict, fact, new_fact; n_inc, has_qns=any(hasqns.(AC)), svd_kwargs...) + break + else + fact = new_fact + end + end + vals = diag(array(new_fact.S)) + (length(vals) == 1 && vals[1]^2 ≤ get(svd_kwargs, :cutoff, 0.0)) && + return nothing, nothing, nothing + #@show flux(dag(CL)*Q*new_fact.U) + #@show flux(new_fact.S), theflux + #@assert flux(new_fact.S) == theflux + #@show inds(new_fact.U) + #@show inds(Q) + #@show inds(dag(CL)) + #@show inds(new_fact.S) + #@show inds(new_fact.V) + #@show inds(dag(CR)) + + return dag(CL) * Q * new_fact.U, new_fact.S, new_fact.V * dag(CR) + #ToDo?: handle non-QN case separately because there it is advisable to start with n_init closer to target maxdim_expand + ##not really an issue anymore since we do *2 increase, so only log number of calls + end + + + + + function rsvd(A::ITensor, linds::Vector{<:Index}, n::Int, p::Int; svd_kwargs...) + rinds = uniqueinds(A, linds) + #ToDo handle empty rinds + #boilerplate matricization of tensor for matrix decomp + CL = combiner(linds) + CR = combiner(rinds) + AC = CL * A * CR + cL = combinedind(CL) + cR = combinedind(CR) + if inds(AC) != (cL, cR) + AC = permute(AC, cL, cR) + end + M = build_guess_matrix(eltype(AC), cR, n, p; theflux=hasqns(AC) ? flux(AC) : nothing) + fact, Q = rsvd_core(AC, M; svd_kwargs...) + return dag(CL) * Q * fact.U, fact.S, fact.V * dag(CR) + end + + function rsvd(A::Vector{<:ITensor}, linds::Vector{<:Index}, n::Int, p::Int; svd_kwargs...) + if length(A)==1 + rinds = uniqueinds(only(A), linds) + else + rinds = uniqueinds(A[end],A[end-1]) + end + #ToDo handle empty rinds + #boilerplate matricization of tensor for matrix decomp + CL = combiner(linds) + CR = combiner(rinds) + @assert !isnothing(commonind(CL,first(A))) + @assert !isnothing(commonind(CR,last(A))) + AC=copy(A) + AC[1] = CL*first(AC) + AC[end] = last(AC)*CR + + cL = combinedind(CL) + cR = combinedind(CR) + theflux = any(hasqns.(AC)) ? reduce(+,flux.(AC)) : nothing + #theflux = mapreduce(flux,+,AC) + M = build_guess_matrix(eltype(first(AC)), cR, n, p; theflux) + fact, Q = rsvd_core(AC, M; svd_kwargs...) + return dag(CL) * Q * fact.U, fact.S, fact.V * dag(CR) + end + + function rsvd_core(AC::ITensor, M; svd_kwargs...) + Q = AC * M + #@show dims(Q) + Q = ITensors.qr(Q, commoninds(AC, Q))[1] + QAC = dag(Q) * AC + fact = svd(QAC, uniqueind(QAC, AC); svd_kwargs...) + return fact, Q + end + + function rsvd_core(AC::ITensors.ITensorNetworkMaps.ITensorNetworkMap, M; svd_kwargs...) + #assumes that we want to do a contraction of M with map over its maps output_inds, i.e. a right-multiply + #thus a transpose is necessary + Q = transpose(AC) * M + Q = ITensors.qr(Q, ITensors.ITensorNetworkMaps.input_inds(AC))[1] + QAC = AC * dag(Q) + @assert typeof(QAC) <: ITensor + #@show inds(QAC) + #@assert !iszero(norm(QAC)) + + fact = svd( + QAC, uniqueind(inds(Q), ITensors.ITensorNetworkMaps.input_inds(AC)); svd_kwargs... + ) + return fact, Q + end + +function rsvd_core(AC::Vector{ITensor}, M; svd_kwargs...) + #assumes that we want to do a contraction of M with map over its maps output_inds, i.e. a right-multiply + #thus a transpose is necessary + @assert !isnothing(commonind(last(AC),M)) + Q = foldr(*,AC;init=M) + Q = ITensors.qr(Q, commoninds(Q,first(AC)))[1] + #@show flux(Q) + #@show nnzblocks(Q) + any(isequal(0),dims(Q)) && return nothing, nothing ,nothing ,nothing + QAC = foldl(*,AC,init=dag(Q)) + #@show inds(QAC) + #@show inds(Q) + #@show inds(first(AC)) + #@show inds(last(AC)) + + @assert typeof(QAC) <: ITensor + + fact = svd( + QAC, commoninds(dag(Q), QAC); svd_kwargs... + ) + return fact, Q + end \ No newline at end of file diff --git a/src/solvers/subspace_expansions/linalg/standard_svd.jl b/src/solvers/subspace_expansions/linalg/standard_svd.jl new file mode 100644 index 00000000..dbabecea --- /dev/null +++ b/src/solvers/subspace_expansions/linalg/standard_svd.jl @@ -0,0 +1,12 @@ +#ToDo: pass use_*_cutoff as kwargs +function _svd_solve_normal(envMap, left_ind; maxdim, cutoff) + M = ITensors.ITensorNetworkMaps.contract(envMap) + # ToDo: infer eltype from envMap + norm(M) ≤ eps(Float64) && return nothing, nothing, nothing + U, S, V = svd( + M, left_ind; maxdim, cutoff=cutoff, use_relative_cutoff=false, use_absolute_cutoff=true + ) + vals = diag(array(S)) + (length(vals) == 1 && vals[1]^2 ≤ cutoff) && return nothing, nothing, nothing + return U, S, V +end diff --git a/src/solvers/subspace_expansions/subspace_expansion.jl b/src/solvers/subspace_expansions/subspace_expansion.jl new file mode 100644 index 00000000..fe73c59f --- /dev/null +++ b/src/solvers/subspace_expansions/subspace_expansion.jl @@ -0,0 +1,398 @@ +function two_site_expansion_updater( + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs, +) + (; maxdim, cutoff, time_step) = region_kwargs + #ToDo: handle timestep==Inf for DMRG case + default_updater_kwargs = (; + svd_func_expand=rsvd_iterative, + maxdim=(maxdim == typemax(Int) ? maxdim : Int(ceil((2.0^(1 / 3))) * maxdim)), + cutoff=isinf(time_step) ? cutoff : cutoff / abs(time_step), # ToDo verify that this is the correct way of scaling the cutoff + use_relative_cutoff=false, + use_absolute_cutoff=true, + ) + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) + + # if on edge return without doing anything + region = first(sweep_plan[which_region_update]) + typeof(region) <: NamedEdge && return init, (;) + region = only(region) + #figure out next region, since we want to expand there + #ToDo account for non-integer indices into which_region_update + next_region = if which_region_update == length(sweep_plan) + nothing + else + first(sweep_plan[which_region_update + 1]) + end + previous_region = + which_region_update == 1 ? nothing : first(sweep_plan[which_region_update - 1]) + isnothing(next_region) && return init, (;) + !(typeof(next_region) <: NamedEdge) && return init, (;) + !(region == src(next_region) || region == dst(next_region)) && return init, (;) + next_vertex = src(next_region) == region ? dst(next_region) : src(next_region) + vp = region => next_vertex + + phi, has_changed = _two_site_expand_core( + init, region, vp; projected_operator!, state!, expand_dir=1, updater_kwargs... + ) + !has_changed && return init, (;) + return phi, (;) +end + +""" +Returns a function which applies Nullspace projector to an ITensor with matching indices without +explicitly constructing the projector as an ITensor or Network. +""" +function implicit_nullspace(A, linkind) + # only works when applied in the direction of the environment tensor, not the other way (due to use of ITensorNetworkMap which is not reversible) + outind = uniqueinds(A, linkind) + inind = outind #? + # ToDo: implement without ITensorNetworkMap + P = ITensors.ITensorNetworkMaps.ITensorNetworkMap( + [prime(dag(A), linkind), prime(A, linkind)], inind, outind + ) + return x::ITensor -> x - P(x) +end + +""" +Performs a local expansion using the two-site projected variance. +The input state should have orthogonality center on a single vertex (region), and phi0 is that site_tensors. +Expansion is performed on vertex that would be visited next (if next vertex is a neighbor of region). +""" +function _two_site_expand_core( + phi0, + region, + vertexpair::Pair; + projected_operator!, + state!, + expand_dir=1, + svd_func_expand, + cutoff, + maxdim, + use_relative_cutoff, + use_absolute_cutoff, +) + # preliminaries + theflux = flux(phi0) + svd_func = svd_func_expand + v1 = first(vertexpair) + v2 = last(vertexpair) + verts = [v1, v2] + n1, n2 = 1, 2 + PH = copy(projected_operator![]) + psi = copy(state![]) + # orthogonalize on the edge + next_edge = edgetype(psi)(vertexpair) + psi, phi = extract_local_tensor(psi, next_edge) + psis = map(n -> psi[n], verts) # extract local site tensors + # ToDo: remove following code lines unless we want to truncate before we expand? + ##this block is replaced by extract_local_tensor --- unless we want to truncate here this is the cleanest. + ##otherwise reinsert the explicit block with truncation args (would allow to remove it from higher-level expander) + #left_inds = uniqueinds(psis[n1], psis[n2]) + #U, S, V = svd(psis[n1], left_inds; lefttags=tags(commonind(psis[n1],psis[n2])), righttags=tags(commonind(psis[n1],psis[n2]))) + #psis[n1]= U + #psi[region]=U + #phi = S*V + + # don't expand if we are already at maxdim + old_linkdim = dim(commonind(first(psis), phi)) + linkinds = map(psi -> commonind(psi, phi), psis) + (old_linkdim >= maxdim) && return phi0, false + + # compute nullspace to the left and right + @timeit_debug timer "nullvector" begin + nullVecs = map(zip(psis, linkinds)) do (psi, linkind) + #return nullspace(psi, linkind; atol=atol) + return implicit_nullspace(psi, linkind) + end + end + + # update the projected operator + PH = set_nsite(PH, 2) + PH = position(PH, psi, verts) + + # build environments + g = underlying_graph(PH) + @timeit_debug timer "build environments" begin + envs = map(zip(verts, psis)) do (n, psi) + return noprime( + mapreduce(*, [v => n for v in neighbors(g, n) if !(v ∈ verts)]; init=psi) do e + return PH.environments[NamedEdge(e)] + end * PH.H[n], + ) + end + end + + # apply the projection into nullspace + envs = [last(nullVecs)(last(envs)), first(nullVecs)(first(envs))] + + # assemble ITensorNetworkMap #ToDo: do not rely on ITensorNetworkMap, simplify logic + ininds = uniqueinds(last(psis), phi) + outinds = uniqueinds(first(psis), phi) + cin = combiner(ininds) + cout = combiner(outinds) + envs = [cin * envs[1], cout * envs[2]] + envMap=[last(envs), phi* first(envs)] + #@show inds(envMap[1]) + #@show inds(envMap[2]) + + + + #envMap = ITensors.ITensorNetworkMaps.ITensorNetworkMap( + # [last(envs), phi, first(envs)], + # uniqueinds(inds(cout), outinds), + # uniqueinds(inds(cin), ininds), + #) + #envMapDag = adjoint(envMap) + + # factorize + @timeit_debug timer "svd_func" begin + if svd_func == ITensorNetworks._svd_solve_normal + U, S, V = svd_func( + envMap, uniqueinds(inds(cout), outinds); maxdim=maxdim - old_linkdim, cutoff=cutoff + ) + elseif svd_func == ITensorNetworks.rsvd_iterative + + U, S, V = svd_func( + envMap, + uniqueinds(inds(cout), outinds); + maxdim=maxdim - old_linkdim, + cutoff=cutoff, + use_relative_cutoff=false, + use_absolute_cutoff=true, + ) + + #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, + #use_absolute_cutoff=true) #this one is for debugging in case we want to test the precontracted version + else + U, S, V = svd_func( + eltype(envMap), + envMap, + envMapDag, + uniqueinds(inds(cout), outinds); + flux=theflux, + maxdim=maxdim - old_linkdim, + cutoff=cutoff, + ) + end + end + + # catch cases when we decompose a map of norm==0.0 + (isnothing(U) || iszero(norm(U))) && return phi0, false + #FIXME: somehow the svd funcs sometimes return empty ITensors instead of nothing, that should be caught in the SVD routines instead... + all(isempty.([U, S, V])) && return phi0, false #ToDo: do not rely on isempty here + + # uncombine indices on the non-link-indices + U *= dag(cout) + V *= dag(cin) + #@show inds(U) + #@show inds(V) + #@show inds(S) + + # direct sum the site tensors + @timeit_debug timer "direct sum" begin + new_psis = map(zip(psis, [U, V])) do (psi, exp_basis) + return ITensors.directsum( + psi => commonind(psi, phi), + exp_basis => uniqueind(exp_basis, psi); + tags=tags(commonind(psi, phi)), + ) + end + end + new_inds = [last(x) for x in new_psis] + new_psis = [first(x) for x in new_psis] + + # extract the expanded linkinds from expanded site tensors and create a zero-tensor + phi_indices = replace( + inds(phi), (commonind(phi, psis[n]) => dag(new_inds[n]) for n in 1:2)... + ) + if hasqns(phi) + new_phi = ITensor(eltype(phi), flux(phi), phi_indices...) + fill!(new_phi, 0.0) #ToDo: Check whether this is really necessary. + else + new_phi = ITensor(eltype(phi), phi_indices...) + end + + # set the new bond tensor elements from the old bond tensor + map(eachindex(phi)) do I + v = phi[I] + !iszero(v) && (return new_phi[I] = v) # I think this line errors without the fill! with zeros above + end + + # apply combiners on linkinds #ToDo: figure out why this is strictly necessary + combiners = map( + new_ind -> combiner(new_ind; tags=tags(new_ind), dir=dir(new_ind)), new_inds + ) + for (v, new_psi, C) in zip(verts, new_psis, combiners) + psi[v] = noprime(new_psi * C) + end + new_phi = dag(first(combiners)) * new_phi * dag(last(combiners)) + + # apply expanded bond tensor to site tensor and reset projected operator to region + psi[v1] *= new_phi + new_phi = psi[v1] + PH = set_nsite(PH, 1) + PH = position(PH, psi, [region]) + + projected_operator![] = PH + state![] = psi + return new_phi, true + ##body end +end + +#= +function _full_expand_core_vertex( + PH, psi, phi, region, svd_func; direction, expand_dir=-1, expander_cache=Any[], maxdim, cutoff, atol=1e-8, kwargs..., +) ###ToDo: adapt to current interface, broken as of now. + #@show cutoff + #enforce expand_dir in the tested direction + @assert expand_dir==-1 + #println("in full expand") + ### only on edges + #(typeof(region)!=NamedEdge{Int}) && return psi, phi, PH + ### only on verts + (typeof(region)==NamedEdge{Int}) && return psi, phi, PH + ## kind of hacky - only works for mps. More general? + n1 = first(region) + theflux=flux(psi[n1]) + #expand_dir=-1 + if direction == 1 + n2 = expand_dir == 1 ? n1 - 1 : n1+1 + else + n2 = expand_dir == 1 ? n1 + 1 : n1-1 + end + + (n2 < 1 || n2 > length(vertices(psi))) && return psi,phi,PH + neutralflux=flux(psi[n2]) + + verts = [n1,n2] + + if isempty(expander_cache) + @warn("building environment of H^2 from scratch!") + + g = underlying_graph(PH.H) + H = vertex_data(data_graph(PH.H)) + + H_dag = swapprime.(prime.(dag.(H)), 1,2, tags = "Site") + H2_vd= replaceprime.(map(*, H, H_dag), 2 => 1) + H2_ed = edge_data(data_graph(PH.H)) + + H2 = TTN(ITensorNetwork(DataGraph(g, H2_vd, H2_ed)), PH.H.ortho_center) + PH2 = ProjTTN(H2) + PH2 = set_nsite(PH2, 2) + + push!(expander_cache, PH2) + end + + PH2 = expander_cache[1] + n1 = verts[2] + n2 = verts[1] + + ###do we really need to make the environments two-site? + ###look into how ProjTTN works + PH2 = position(PH2, psi, [n1,n2]) + PH = position(PH, psi, [n1,n2]) + + psi1 = psi[n1] + psi2 = psi[n2] #phi + old_linkdim = dim(commonind(psi1, psi2)) + + # don't expand if we are already at maxdim + ## make more transparent that this is not the normal maxdim arg but maxdim_expand + ## when would this even happen? + #@show old_linkdim, maxdim + + old_linkdim >= maxdim && return psi, phi, PH + #@show "expandin", maxdim, old_linkdim, maxdim-old_linkdim + # compute nullspace to the left and right + linkind_l = commonind(psi1, psi2) + nullVec = implicit_nullspace(psi1, linkind_l) + + # if nullspace is empty (happen's for product states with QNs) + ## ToDo implement norm or equivalent for a generic LinearMap (i guess via sampling a random vector) + #norm(nullVec) == 0.0 && return psi, phi, PH + + ## compute both environments + g = underlying_graph(PH) + + @timeit_debug timer "build environments" begin + env1 = noprime(mapreduce(*, [v => n1 for v in neighbors(g, n1) if v != n2]; init = psi1) do e + return PH.environments[NamedEdge(e)] + end *PH.H[n1] + ) + env2p = noprime(mapreduce(*, [v => n2 for v in neighbors(g, n2) if v != n1]; init = psi2) do e + return PH.environments[NamedEdge(e)] + end *PH.H[n2] + ) + + env2 = mapreduce(*, [v => n2 for v in neighbors(g, n2) if v != n1]; init = psi2) do e + return PH2.environments[NamedEdge(e)] + end * PH2.H[n2]*prime(dag(psi2)) + end + + env1=nullVec(env1) + + outinds = uniqueinds(psi1,psi2) + ininds = dag.(outinds) + cout=combiner(outinds) + env1 *= cout + env2p *= prime(dag(psi[n2]),commonind(dag(psi[n2]),dag(psi[n1]))) + env2p2= replaceprime(env2p * replaceprime(dag(env2p),0=>2),2=>1) + + envMap = ITensors.ITensorNetworkMaps.ITensorNetworkMap([prime(dag(env1)),(env2-env2p2),env1] , prime(dag(uniqueinds(cout,outinds))), uniqueinds(cout,outinds)) + envMapDag=adjoint(envMap) + # svd-decomposition + @timeit_debug timer "svd_func" begin + if svd_func==ITensorNetworks._svd_solve_normal + U,S,_ = svd_func(envMap, uniqueinds(inds(cout),outinds); maxdim=maxdim-old_linkdim, cutoff=cutoff) + elseif svd_func==ITensorNetworks.rsvd_iterative + #@show theflux + envMap=transpose(envMap) + U,S,_ = svd_func(eltype(first(envMap.itensors)),envMap,ITensors.ITensorNetworkMaps.input_inds(envMap);theflux=neutralflux, maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, + use_absolute_cutoff=true) + #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);theflux=theflux, maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, + #use_absolute_cutoff=true) + else + U,S,_= svd_func(eltype(envMap),envMap,envMapDag,uniqueinds(cout,outinds); flux=neutralflux, maxdim=maxdim-old_linkdim, cutoff=cutoff) + end + end + isnothing(U) && return psi,phi,PH + ###FIXME: somehow the svd funcs sometimes return empty ITensors instead of nothing, that should be caught in the SVD routines instead... + all(isempty.([U,S])) && return psi, phi, PH + @assert dim(commonind(U, S)) ≤ maxdim-old_linkdim + nullVec = dag(cout)*U + new_psi1, new_ind1 = ITensors.directsum( + psi1 => uniqueinds(psi1, nullVec), nullVec => uniqueinds(nullVec, psi1); tags=(tags(commonind(psi1,phi)),) + ) + new_ind1 = new_ind1[1] + @assert dim(new_ind1) <= maxdim + + Cl = combiner(new_ind1; tags=tags(new_ind1), dir=dir(new_ind1)) + + phi_indices = replace(inds(phi), commonind(phi,psi1) => dag(new_ind1)) + + if hasqns(phi) + new_phi=ITensor(eltype(phi),flux(phi), phi_indices...) + fill!(new_phi,0.0) + else + new_phi = ITensor(eltype(phi), phi_indices...) + end + + map(eachindex(phi)) do I + v = phi[I] + !iszero(v) && (return new_phi[I]=v) + end + + psi[n1] = noprime(new_psi1*Cl) + new_phi = dag(Cl)*new_phi + + return psi, new_phi, PH +end +=# From 6a97fa822791eb02769e38760eaed91aaf5cd6b9 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Tue, 2 Apr 2024 12:16:58 -0400 Subject: [PATCH 02/19] Add subspace expander, compose_updater, truncating extracter. --- src/solvers/defaults.jl | 2 + src/solvers/extract/extract.jl | 26 +++ src/solvers/local_solvers/expand.jl | 248 ++++++++++++++++++++++++++++ 3 files changed, 276 insertions(+) create mode 100644 src/solvers/local_solvers/expand.jl diff --git a/src/solvers/defaults.jl b/src/solvers/defaults.jl index 9e901af3..52303411 100644 --- a/src/solvers/defaults.jl +++ b/src/solvers/defaults.jl @@ -5,6 +5,8 @@ default_extracter() = default_extracter default_inserter() = default_inserter default_checkdone() = (; kws...) -> false default_transform_operator() = nothing +default_scale_cutoff_by_timestep(cutoff;time_step) = isinf(time_step) ? cutoff : cutoff / abs(time_step) + function default_region_printer(; cutoff, maxdim, diff --git a/src/solvers/extract/extract.jl b/src/solvers/extract/extract.jl index feb57c2f..abce0829 100644 --- a/src/solvers/extract/extract.jl +++ b/src/solvers/extract/extract.jl @@ -24,3 +24,29 @@ function default_extracter(state, projected_operator, region, ortho; internal_kw projected_operator = position(projected_operator, state, region) return state, projected_operator, local_tensor end + +function extract_and_truncate(state,projected_operator, region, ortho; maxdim=nothing, cutoff=nothing,internal_kwargs) + #ToDo: move default into function signature? + #svd_kwargs = NamedTuple() + #svd_kwargs = isnothing(maxdim) ? svd_kwargs : merge(svd_kwargs,(;maxdim)) + #svd_kwargs = isnothing(cutoff) ? svd_kwargs : merge(svd_kwargs,(;cutoff)) + svd_kwargs= (; + (isnothing(maxdim) ? (;) : (;maxdim))..., + (isnothing(cutoff) ? (;) : (;cutoff))... + ) + state = orthogonalize(state, ortho) + if isa(region, AbstractEdge) + other_vertex = only(setdiff(support(region), [ortho])) + left_inds = uniqueinds(state[ortho], state[other_vertex]) + #ToDo: replace with call to factorize + U, S, V = svd( + state[ortho], left_inds; lefttags=tags(state, region), righttags=tags(state, region), svd_kwargs... + ) + state[ortho] = U + local_tensor = S * V + else + local_tensor = prod(state[v] for v in region) + end + projected_operator = position(projected_operator, state, region) + return state, projected_operator, local_tensor +end \ No newline at end of file diff --git a/src/solvers/local_solvers/expand.jl b/src/solvers/local_solvers/expand.jl new file mode 100644 index 00000000..f4867ff5 --- /dev/null +++ b/src/solvers/local_solvers/expand.jl @@ -0,0 +1,248 @@ +#ToDo: is this the correct scaling for infinite timestep? DMRG will not have infinite timestep, unless we set it explicitly +#ToDo: implement this as a closure, to be constructed in sweep_plan (where time_step is known) +#scale_cutoff_by_timestep(cutoff;time_step) = isinf(time_step) ? cutoff : cutoff / abs(time_step) + +function two_site_expansion_updater( + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + internal_kwargs, + svd_func_expand=ITensorNetworks.rsvd_iterative, + maxdim, + maxdim_func= (arg;kwargs...) -> identity(arg), + cutoff, + cutoff_func= (arg;kwargs...) -> identity(arg), + use_relative_cutoff=false, + use_absolute_cutoff=true, + ) + maxdim=maxdim_func(maxdim;internal_kwargs...) + cutoff=cutoff_func(cutoff;internal_kwargs...) + region = first(sweep_plan[which_region_update]) + typeof(region) <: NamedEdge && return init, (;) + region = only(region) + #figure out next region, since we want to expand there + #ToDo account for non-integer indices into which_region_update + next_region = if which_region_update == length(sweep_plan) + nothing + else + first(sweep_plan[which_region_update + 1]) + end + previous_region = + which_region_update == 1 ? nothing : first(sweep_plan[which_region_update - 1]) + isnothing(next_region) && return init, (;) + !(typeof(next_region) <: NamedEdge) && return init, (;) + !(region == src(next_region) || region == dst(next_region)) && return init, (;) + next_vertex = src(next_region) == region ? dst(next_region) : src(next_region) + vp = region => next_vertex + + phi, has_changed = _two_site_expand_core( + init, region, vp; projected_operator!, state!, expand_dir=1, updater_kwargs... + ) + !has_changed && return init, (;) + return phi, (;) + end + + """ +Returns a function which applies Nullspace projector to an ITensor with matching indices without +explicitly constructing the projector as an ITensor or Network. +""" +function implicit_nullspace(A, linkind) + # only works when applied in the direction of the environment tensor, not the other way (due to use of ITensorNetworkMap which is not reversible) + outind = uniqueinds(A, linkind) + inind = outind #? + # ToDo: implement without ITensorNetworkMap + P = ITensors.ITensorNetworkMaps.ITensorNetworkMap( + [prime(dag(A), linkind), prime(A, linkind)], inind, outind + ) + return x::ITensor -> x - P(x) +end + +""" +Performs a local expansion using the two-site projected variance. +The input state should have orthogonality center on a single vertex (region), and phi0 is that site_tensors. +Expansion is performed on vertex that would be visited next (if next vertex is a neighbor of region). +""" +function _two_site_expand_core( + phi0, + region, + vertexpair::Pair; + projected_operator!, + state!, + expand_dir=1, + svd_func_expand, + cutoff, + maxdim, + use_relative_cutoff, + use_absolute_cutoff, +) + # preliminaries + theflux = flux(phi0) + svd_func = svd_func_expand + v1 = first(vertexpair) + v2 = last(vertexpair) + verts = [v1, v2] + n1, n2 = 1, 2 + PH = copy(projected_operator![]) + psi = copy(state![]) + # orthogonalize on the edge + next_edge = edgetype(psi)(vertexpair) + psi, phi = extract_local_tensor(psi, next_edge) + psis = map(n -> psi[n], verts) # extract local site tensors + # ToDo: remove following code lines unless we want to truncate before we expand? + ##this block is replaced by extract_local_tensor --- unless we want to truncate here this is the cleanest. + ##otherwise reinsert the explicit block with truncation args (would allow to remove it from higher-level expander) + #left_inds = uniqueinds(psis[n1], psis[n2]) + #U, S, V = svd(psis[n1], left_inds; lefttags=tags(commonind(psis[n1],psis[n2])), righttags=tags(commonind(psis[n1],psis[n2]))) + #psis[n1]= U + #psi[region]=U + #phi = S*V + + # don't expand if we are already at maxdim + old_linkdim = dim(commonind(first(psis), phi)) + linkinds = map(psi -> commonind(psi, phi), psis) + (old_linkdim >= maxdim) && return phi0, false + + # compute nullspace to the left and right + @timeit_debug timer "nullvector" begin + nullVecs = map(zip(psis, linkinds)) do (psi, linkind) + #return nullspace(psi, linkind; atol=atol) + return implicit_nullspace(psi, linkind) + end + end + + # update the projected operator + PH = set_nsite(PH, 2) + PH = position(PH, psi, verts) + + # build environments + g = underlying_graph(PH) + @timeit_debug timer "build environments" begin + envs = map(zip(verts, psis)) do (n, psi) + return noprime( + mapreduce(*, [v => n for v in neighbors(g, n) if !(v ∈ verts)]; init=psi) do e + return PH.environments[NamedEdge(e)] + end * PH.H[n], + ) + end + end + + # apply the projection into nullspace + envs = [last(nullVecs)(last(envs)), first(nullVecs)(first(envs))] + + # assemble ITensorNetworkMap #ToDo: do not rely on ITensorNetworkMap, simplify logic + ininds = uniqueinds(last(psis), phi) + outinds = uniqueinds(first(psis), phi) + cin = combiner(ininds) + cout = combiner(outinds) + envs = [cin * envs[1], cout * envs[2]] + envMap=[last(envs), phi* first(envs)] + #@show inds(envMap[1]) + #@show inds(envMap[2]) + + + + #envMap = ITensors.ITensorNetworkMaps.ITensorNetworkMap( + # [last(envs), phi, first(envs)], + # uniqueinds(inds(cout), outinds), + # uniqueinds(inds(cin), ininds), + #) + #envMapDag = adjoint(envMap) + + # factorize + @timeit_debug timer "svd_func" begin + if svd_func == ITensorNetworks._svd_solve_normal + U, S, V = svd_func( + envMap, uniqueinds(inds(cout), outinds); maxdim=maxdim - old_linkdim, cutoff=cutoff + ) + elseif svd_func == ITensorNetworks.rsvd_iterative + + U, S, V = svd_func( + envMap, + uniqueinds(inds(cout), outinds); + maxdim=maxdim - old_linkdim, + cutoff=cutoff, + use_relative_cutoff=false, + use_absolute_cutoff=true, + ) + + #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, + #use_absolute_cutoff=true) #this one is for debugging in case we want to test the precontracted version + else + U, S, V = svd_func( + eltype(envMap), + envMap, + envMapDag, + uniqueinds(inds(cout), outinds); + flux=theflux, + maxdim=maxdim - old_linkdim, + cutoff=cutoff, + ) + end + end + + # catch cases when we decompose a map of norm==0.0 + (isnothing(U) || iszero(norm(U))) && return phi0, false + #FIXME: somehow the svd funcs sometimes return empty ITensors instead of nothing, that should be caught in the SVD routines instead... + all(isempty.([U, S, V])) && return phi0, false #ToDo: do not rely on isempty here + + # uncombine indices on the non-link-indices + U *= dag(cout) + V *= dag(cin) + #@show inds(U) + #@show inds(V) + #@show inds(S) + + # direct sum the site tensors + @timeit_debug timer "direct sum" begin + new_psis = map(zip(psis, [U, V])) do (psi, exp_basis) + return ITensors.directsum( + psi => commonind(psi, phi), + exp_basis => uniqueind(exp_basis, psi); + tags=tags(commonind(psi, phi)), + ) + end + end + new_inds = [last(x) for x in new_psis] + new_psis = [first(x) for x in new_psis] + + # extract the expanded linkinds from expanded site tensors and create a zero-tensor + phi_indices = replace( + inds(phi), (commonind(phi, psis[n]) => dag(new_inds[n]) for n in 1:2)... + ) + if hasqns(phi) + new_phi = ITensor(eltype(phi), flux(phi), phi_indices...) + fill!(new_phi, 0.0) #ToDo: Check whether this is really necessary. + else + new_phi = ITensor(eltype(phi), phi_indices...) + end + + # set the new bond tensor elements from the old bond tensor + map(eachindex(phi)) do I + v = phi[I] + !iszero(v) && (return new_phi[I] = v) # I think this line errors without the fill! with zeros above + end + + # apply combiners on linkinds #ToDo: figure out why this is strictly necessary + combiners = map( + new_ind -> combiner(new_ind; tags=tags(new_ind), dir=dir(new_ind)), new_inds + ) + for (v, new_psi, C) in zip(verts, new_psis, combiners) + psi[v] = noprime(new_psi * C) + end + new_phi = dag(first(combiners)) * new_phi * dag(last(combiners)) + + # apply expanded bond tensor to site tensor and reset projected operator to region + psi[v1] *= new_phi + new_phi = psi[v1] + PH = set_nsite(PH, 1) + PH = position(PH, psi, [region]) + + projected_operator![] = PH + state![] = psi + return new_phi, true + ##body end +end \ No newline at end of file From ce1f1709d2b2cbd2b5e8a6d7495f4f18f2afeba7 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Tue, 2 Apr 2024 12:17:24 -0400 Subject: [PATCH 03/19] Actually add compose_updater now. --- src/solvers/solver_utils.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/solvers/solver_utils.jl b/src/solvers/solver_utils.jl index 68911a65..0f213a37 100644 --- a/src/solvers/solver_utils.jl +++ b/src/solvers/solver_utils.jl @@ -86,3 +86,26 @@ function cache_operator_to_disk( end return operator end + +#ToDo: Move? This belongs more into local_solvers +function compose_updaters(; kwargs...#solver=eigsolve_updater, expander=local_subspace_expansion) + funcs=values(kwargs) + kwarg_symbols=Symbol.(keys(kwargs),"_kwargs") + info_symbols=Symbol.(keys(kwargs),"_info") + + function composed_updater(init; kwargs...) + info=(;) + kwargs_for_updaters=map(x -> kwargs[x], kwarg_symbols) + other_kwargs=Base.structdiff(kwargs,NamedTuple(zip(kwarg_symbols,kwargs_for_updaters))) + + for (func,kwargs_for_updater,info_symbol) in zip(funcs,kwargs_for_updaters,info_symbols) + init, new_info=func(init; + other_kwargs..., + kwargs_for_updater... + ) + #figure out a better way to handle composing the info? + info=(;info...,NamedTuple((info_symbol=>new_info,))...) + + return init, info + return composed_updater +end \ No newline at end of file From 99a95a788750d99dae9e8fb6e1f31c797a6954b4 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 3 Apr 2024 18:46:45 -0400 Subject: [PATCH 04/19] update includes --- src/ITensorNetworks.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index d5e61da6..7e4e5a6e 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -56,6 +56,10 @@ include("solvers/solver_utils.jl") include("solvers/defaults.jl") include("solvers/insert/insert.jl") include("solvers/extract/extract.jl") +include("solvers/extract/subspace_expansions/subspace_expansion.jl") +include("solvers/extract/subspace_expansions/linalg/standard_svd.jl") +include("solvers/extract/subspace_expansions/linalg/rsvd_aux.jl") +include("solvers/extract/subspace_expansions/linalg/rsvd_linalg.jl") include("solvers/alternating_update/alternating_update.jl") include("solvers/alternating_update/region_update.jl") include("solvers/tdvp.jl") From 4393544e8e0d9d49a4705287d08b16d6b6894d16 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 3 Apr 2024 18:47:58 -0400 Subject: [PATCH 05/19] Slurp kwargs in scaling by cutoff, and cosmetic change to compose_updater. --- src/solvers/defaults.jl | 2 +- src/solvers/solver_utils.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/defaults.jl b/src/solvers/defaults.jl index ad4b57e3..13dbafe5 100644 --- a/src/solvers/defaults.jl +++ b/src/solvers/defaults.jl @@ -7,7 +7,7 @@ default_extracter() = default_extracter default_inserter() = default_inserter default_checkdone() = (; kws...) -> false default_transform_operator() = nothing -default_scale_cutoff_by_timestep(cutoff;time_step) = isinf(time_step) ? cutoff : cutoff / abs(time_step) +default_scale_cutoff_by_timestep(cutoff;time_step, kwargs...) = isinf(time_step) ? cutoff : cutoff / abs(time_step) function default_region_printer(; cutoff, diff --git a/src/solvers/solver_utils.jl b/src/solvers/solver_utils.jl index 918a70c8..74507d92 100644 --- a/src/solvers/solver_utils.jl +++ b/src/solvers/solver_utils.jl @@ -97,7 +97,7 @@ function compose_updaters(; kwargs...#solver=eigsolve_updater, expander=local_su function composed_updater(init; kwargs...) info=(;) kwargs_for_updaters=map(x -> kwargs[x], kwarg_symbols) - other_kwargs=Base.structdiff(kwargs,NamedTuple(zip(kwarg_symbols,kwargs_for_updaters))) + other_kwargs=Base.structdiff(kwargs,NamedTuple(kwarg_symbols .=> kwargs_for_updaters)) for (func,kwargs_for_updater,info_symbol) in zip(funcs,kwargs_for_updaters,info_symbols) init, new_info=func(init; From 39c6529762152bc11d11cf53f724b570a3a8a34e Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 10:11:18 -0400 Subject: [PATCH 06/19] Bump compat entry for ITensors. --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index aa488a06..c78198a6 100644 --- a/Project.toml +++ b/Project.toml @@ -54,12 +54,12 @@ DocStringExtensions = "0.8, 0.9" EinExprs = "0.6.4" Graphs = "1.8" GraphsFlows = "0.1.1" -ITensors = "0.3.58" +ITensors = "0.4" IsApprox = "0.1" IterTools = "1.4.0" KrylovKit = "0.6, 0.7" -NamedGraphs = "0.1.23" NDTensors = "0.2, 0.3" +NamedGraphs = "0.1.23" Observers = "0.2" PackageExtensionCompat = "1" Requires = "1.3" From 22565dfff5fe9f4a19dc4811a631182df1b466e8 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 10:23:19 -0400 Subject: [PATCH 07/19] Fix typos in compose_updater. --- src/solvers/solver_utils.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solvers/solver_utils.jl b/src/solvers/solver_utils.jl index 74507d92..d77e2894 100644 --- a/src/solvers/solver_utils.jl +++ b/src/solvers/solver_utils.jl @@ -89,11 +89,10 @@ function cache_operator_to_disk( end #ToDo: Move? This belongs more into local_solvers -function compose_updaters(; kwargs...#solver=eigsolve_updater, expander=local_subspace_expansion) +function compose_updaters(;kwargs...) funcs=values(kwargs) kwarg_symbols=Symbol.(keys(kwargs),"_kwargs") info_symbols=Symbol.(keys(kwargs),"_info") - function composed_updater(init; kwargs...) info=(;) kwargs_for_updaters=map(x -> kwargs[x], kwarg_symbols) @@ -106,7 +105,8 @@ function compose_updaters(; kwargs...#solver=eigsolve_updater, expander=local_su ) #figure out a better way to handle composing the info? info=(;info...,NamedTuple((info_symbol=>new_info,))...) - + end return init, info + end return composed_updater end \ No newline at end of file From b5428817a1a5e5207623b1616960a3b0b3e7e120 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 10:23:44 -0400 Subject: [PATCH 08/19] Comment out instrumentation in subspace expansion. --- .../subspace_expansions/subspace_expansion.jl | 100 +++++++++--------- 1 file changed, 50 insertions(+), 50 deletions(-) diff --git a/src/solvers/subspace_expansions/subspace_expansion.jl b/src/solvers/subspace_expansions/subspace_expansion.jl index fe73c59f..55fbce5e 100644 --- a/src/solvers/subspace_expansions/subspace_expansion.jl +++ b/src/solvers/subspace_expansions/subspace_expansion.jl @@ -107,12 +107,12 @@ function _two_site_expand_core( (old_linkdim >= maxdim) && return phi0, false # compute nullspace to the left and right - @timeit_debug timer "nullvector" begin - nullVecs = map(zip(psis, linkinds)) do (psi, linkind) - #return nullspace(psi, linkind; atol=atol) - return implicit_nullspace(psi, linkind) - end + #@timeit_debug timer "nullvector" begin + nullVecs = map(zip(psis, linkinds)) do (psi, linkind) + #return nullspace(psi, linkind; atol=atol) + return implicit_nullspace(psi, linkind) end + #end # update the projected operator PH = set_nsite(PH, 2) @@ -120,15 +120,15 @@ function _two_site_expand_core( # build environments g = underlying_graph(PH) - @timeit_debug timer "build environments" begin - envs = map(zip(verts, psis)) do (n, psi) - return noprime( - mapreduce(*, [v => n for v in neighbors(g, n) if !(v ∈ verts)]; init=psi) do e - return PH.environments[NamedEdge(e)] - end * PH.H[n], - ) - end + #@timeit_debug timer "build environments" begin + envs = map(zip(verts, psis)) do (n, psi) + return noprime( + mapreduce(*, [v => n for v in neighbors(g, n) if !(v ∈ verts)]; init=psi) do e + return PH.environments[NamedEdge(e)] + end * PH.H[n], + ) end + #end # apply the projection into nullspace envs = [last(nullVecs)(last(envs)), first(nullVecs)(first(envs))] @@ -153,36 +153,36 @@ function _two_site_expand_core( #envMapDag = adjoint(envMap) # factorize - @timeit_debug timer "svd_func" begin - if svd_func == ITensorNetworks._svd_solve_normal - U, S, V = svd_func( - envMap, uniqueinds(inds(cout), outinds); maxdim=maxdim - old_linkdim, cutoff=cutoff - ) - elseif svd_func == ITensorNetworks.rsvd_iterative - - U, S, V = svd_func( - envMap, - uniqueinds(inds(cout), outinds); - maxdim=maxdim - old_linkdim, - cutoff=cutoff, - use_relative_cutoff=false, - use_absolute_cutoff=true, - ) - - #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, - #use_absolute_cutoff=true) #this one is for debugging in case we want to test the precontracted version - else - U, S, V = svd_func( - eltype(envMap), - envMap, - envMapDag, - uniqueinds(inds(cout), outinds); - flux=theflux, - maxdim=maxdim - old_linkdim, - cutoff=cutoff, - ) - end + #@timeit_debug timer "svd_func" begin + if svd_func == ITensorNetworks._svd_solve_normal + U, S, V = svd_func( + envMap, uniqueinds(inds(cout), outinds); maxdim=maxdim - old_linkdim, cutoff=cutoff + ) + elseif svd_func == ITensorNetworks.rsvd_iterative + + U, S, V = svd_func( + envMap, + uniqueinds(inds(cout), outinds); + maxdim=maxdim - old_linkdim, + cutoff=cutoff, + use_relative_cutoff=false, + use_absolute_cutoff=true, + ) + + #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, + #use_absolute_cutoff=true) #this one is for debugging in case we want to test the precontracted version + else + U, S, V = svd_func( + eltype(envMap), + envMap, + envMapDag, + uniqueinds(inds(cout), outinds); + flux=theflux, + maxdim=maxdim - old_linkdim, + cutoff=cutoff, + ) end + #end # catch cases when we decompose a map of norm==0.0 (isnothing(U) || iszero(norm(U))) && return phi0, false @@ -197,15 +197,15 @@ function _two_site_expand_core( #@show inds(S) # direct sum the site tensors - @timeit_debug timer "direct sum" begin - new_psis = map(zip(psis, [U, V])) do (psi, exp_basis) - return ITensors.directsum( - psi => commonind(psi, phi), - exp_basis => uniqueind(exp_basis, psi); - tags=tags(commonind(psi, phi)), - ) - end + #@timeit_debug timer "direct sum" begin + new_psis = map(zip(psis, [U, V])) do (psi, exp_basis) + return ITensors.directsum( + psi => commonind(psi, phi), + exp_basis => uniqueind(exp_basis, psi); + tags=tags(commonind(psi, phi)), + ) end + #end new_inds = [last(x) for x in new_psis] new_psis = [first(x) for x in new_psis] From 94efab9c1e94bea42d2a20f40645821779928dd0 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 10:24:02 -0400 Subject: [PATCH 09/19] Fix wrong paths in includes. --- src/ITensorNetworks.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 74ee2a1e..5beb14f5 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -53,10 +53,10 @@ include("solvers/solver_utils.jl") include("solvers/defaults.jl") include("solvers/insert/insert.jl") include("solvers/extract/extract.jl") -include("solvers/extract/subspace_expansions/subspace_expansion.jl") -include("solvers/extract/subspace_expansions/linalg/standard_svd.jl") -include("solvers/extract/subspace_expansions/linalg/rsvd_aux.jl") -include("solvers/extract/subspace_expansions/linalg/rsvd_linalg.jl") +include("solvers/subspace_expansions/subspace_expansion.jl") +include("solvers/subspace_expansions/linalg/standard_svd.jl") +include("solvers/subspace_expansions/linalg/rsvd_aux.jl") +include("solvers/subspace_expansions/linalg/rsvd_linalg.jl") include("solvers/alternating_update/alternating_update.jl") include("solvers/alternating_update/region_update.jl") include("solvers/tdvp.jl") From 679da47665faaf90d8ffa91af72c72552c754156 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 17:03:16 -0400 Subject: [PATCH 10/19] Fix lacking conversion to NamedTuple in compose_updaters. --- src/solvers/solver_utils.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/solver_utils.jl b/src/solvers/solver_utils.jl index d77e2894..ad42fea1 100644 --- a/src/solvers/solver_utils.jl +++ b/src/solvers/solver_utils.jl @@ -93,10 +93,11 @@ function compose_updaters(;kwargs...) funcs=values(kwargs) kwarg_symbols=Symbol.(keys(kwargs),"_kwargs") info_symbols=Symbol.(keys(kwargs),"_info") + function composed_updater(init; kwargs...) info=(;) kwargs_for_updaters=map(x -> kwargs[x], kwarg_symbols) - other_kwargs=Base.structdiff(kwargs,NamedTuple(kwarg_symbols .=> kwargs_for_updaters)) + other_kwargs=Base.structdiff((;kwargs...),NamedTuple(kwarg_symbols .=> kwargs_for_updaters)) for (func,kwargs_for_updater,info_symbol) in zip(funcs,kwargs_for_updaters,info_symbols) init, new_info=func(init; From 4dac808e1674e9ea7a705a865076238246fb7eb6 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 17:06:08 -0400 Subject: [PATCH 11/19] Fix bugs in expand.jl and two-site expansion --- src/solvers/local_solvers/expand.jl | 99 ++++++++----------- .../local_solvers/expand_and_exponentiate.jl | 94 ------------------ 2 files changed, 42 insertions(+), 151 deletions(-) delete mode 100644 src/solvers/local_solvers/expand_and_exponentiate.jl diff --git a/src/solvers/local_solvers/expand.jl b/src/solvers/local_solvers/expand.jl index f4867ff5..5b8a0f89 100644 --- a/src/solvers/local_solvers/expand.jl +++ b/src/solvers/local_solvers/expand.jl @@ -19,6 +19,7 @@ function two_site_expansion_updater( use_relative_cutoff=false, use_absolute_cutoff=true, ) + #return init, (;) maxdim=maxdim_func(maxdim;internal_kwargs...) cutoff=cutoff_func(cutoff;internal_kwargs...) region = first(sweep_plan[which_region_update]) @@ -40,7 +41,12 @@ function two_site_expansion_updater( vp = region => next_vertex phi, has_changed = _two_site_expand_core( - init, region, vp; projected_operator!, state!, expand_dir=1, updater_kwargs... + init, region, vp; projected_operator!, state!, expand_dir=1, + svd_func_expand, + maxdim, + cutoff, + use_relative_cutoff, + use_absolute_cutoff, ) !has_changed && return init, (;) return phi, (;) @@ -90,89 +96,73 @@ function _two_site_expand_core( psi = copy(state![]) # orthogonalize on the edge next_edge = edgetype(psi)(vertexpair) - psi, phi = extract_local_tensor(psi, next_edge) + #update of environment not strictly necessary here + psi, PH, phi = default_extracter(psi, PH, next_edge, v1;internal_kwargs=(;))#default_extracter(psi,next_edge)#extract_local_tensor(psi, next_edge) + psis = map(n -> psi[n], verts) # extract local site tensors - # ToDo: remove following code lines unless we want to truncate before we expand? - ##this block is replaced by extract_local_tensor --- unless we want to truncate here this is the cleanest. - ##otherwise reinsert the explicit block with truncation args (would allow to remove it from higher-level expander) - #left_inds = uniqueinds(psis[n1], psis[n2]) - #U, S, V = svd(psis[n1], left_inds; lefttags=tags(commonind(psis[n1],psis[n2])), righttags=tags(commonind(psis[n1],psis[n2]))) - #psis[n1]= U - #psi[region]=U - #phi = S*V - - # don't expand if we are already at maxdim + + #if we are already past the maximum allowed bond dimension, return the input + #ToDo: move up so we don't do any orthogonalization old_linkdim = dim(commonind(first(psis), phi)) linkinds = map(psi -> commonind(psi, phi), psis) (old_linkdim >= maxdim) && return phi0, false # compute nullspace to the left and right - @timeit_debug timer "nullvector" begin + #@timeit_debug timer "nullvector" begin nullVecs = map(zip(psis, linkinds)) do (psi, linkind) #return nullspace(psi, linkind; atol=atol) - return implicit_nullspace(psi, linkind) + return ITensorNetworks.implicit_nullspace(psi, linkind) end - end + #end - # update the projected operator + # update the projected operator + #ToDo: may not be necessary if we use the extracter anyway PH = set_nsite(PH, 2) PH = position(PH, psi, verts) # build environments g = underlying_graph(PH) - @timeit_debug timer "build environments" begin + #@timeit_debug timer "build environments" begin envs = map(zip(verts, psis)) do (n, psi) return noprime( mapreduce(*, [v => n for v in neighbors(g, n) if !(v ∈ verts)]; init=psi) do e return PH.environments[NamedEdge(e)] - end * PH.H[n], + end * PH.operator[n], ) end - end + #end - # apply the projection into nullspace - envs = [last(nullVecs)(last(envs)), first(nullVecs)(first(envs))] - - # assemble ITensorNetworkMap #ToDo: do not rely on ITensorNetworkMap, simplify logic + # apply the projection into nullspace, + # FIXME: think through again whether we really want to evaluate null space projector already + envs=[ first(nullVecs)(first(envs)),last(nullVecs)(last(envs))] ininds = uniqueinds(last(psis), phi) outinds = uniqueinds(first(psis), phi) cin = combiner(ininds) cout = combiner(outinds) - envs = [cin * envs[1], cout * envs[2]] - envMap=[last(envs), phi* first(envs)] - #@show inds(envMap[1]) - #@show inds(envMap[2]) - - - - #envMap = ITensors.ITensorNetworkMaps.ITensorNetworkMap( - # [last(envs), phi, first(envs)], - # uniqueinds(inds(cout), outinds), - # uniqueinds(inds(cin), ininds), - #) - #envMapDag = adjoint(envMap) + envMap=[cout * envs[1], phi* (cin * envs[2])] # factorize - @timeit_debug timer "svd_func" begin + #FIXME: remove specialization on svd_func, define interface for svd_funcs + #@timeit_debug timer "svd_func" begin if svd_func == ITensorNetworks._svd_solve_normal - U, S, V = svd_func( + U, S, V = svd_func( envMap, uniqueinds(inds(cout), outinds); maxdim=maxdim - old_linkdim, cutoff=cutoff - ) + ) elseif svd_func == ITensorNetworks.rsvd_iterative - - U, S, V = svd_func( + + U, S, V = svd_func( envMap, uniqueinds(inds(cout), outinds); maxdim=maxdim - old_linkdim, cutoff=cutoff, use_relative_cutoff=false, use_absolute_cutoff=true, - ) - - #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, - #use_absolute_cutoff=true) #this one is for debugging in case we want to test the precontracted version + ) + + #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, + #use_absolute_cutoff=true) #this one is for debugging in case we want to test the precontracted version else - U, S, V = svd_func( + U, S, V = svd_func( eltype(envMap), envMap, envMapDag, @@ -180,9 +170,9 @@ function _two_site_expand_core( flux=theflux, maxdim=maxdim - old_linkdim, cutoff=cutoff, - ) - end - end + ) +end + #end # catch cases when we decompose a map of norm==0.0 (isnothing(U) || iszero(norm(U))) && return phi0, false @@ -192,12 +182,9 @@ function _two_site_expand_core( # uncombine indices on the non-link-indices U *= dag(cout) V *= dag(cin) - #@show inds(U) - #@show inds(V) - #@show inds(S) - + # direct sum the site tensors - @timeit_debug timer "direct sum" begin + #@timeit_debug timer "direct sum" begin new_psis = map(zip(psis, [U, V])) do (psi, exp_basis) return ITensors.directsum( psi => commonind(psi, phi), @@ -205,7 +192,7 @@ function _two_site_expand_core( tags=tags(commonind(psi, phi)), ) end - end + #end new_inds = [last(x) for x in new_psis] new_psis = [first(x) for x in new_psis] @@ -240,9 +227,7 @@ function _two_site_expand_core( new_phi = psi[v1] PH = set_nsite(PH, 1) PH = position(PH, psi, [region]) - projected_operator![] = PH state![] = psi return new_phi, true - ##body end end \ No newline at end of file diff --git a/src/solvers/local_solvers/expand_and_exponentiate.jl b/src/solvers/local_solvers/expand_and_exponentiate.jl deleted file mode 100644 index 6df90efe..00000000 --- a/src/solvers/local_solvers/expand_and_exponentiate.jl +++ /dev/null @@ -1,94 +0,0 @@ -function local_expand_and_exponentiate_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - region_kwargs, - updater_kwargs, -) - # expansion - expand_kwargs = updater_kwargs.expand_kwargs # defaults are set in the expansion_updater - expanded_init, _ = two_site_expansion_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - region_kwargs, - updater_kwargs=expand_kwargs, - ) - - # exponentiate - # ToDo: also call exponentiate_updater, instead of reimplementing it here - # - default_exponentiate_kwargs = (; - krylovdim=30, - maxiter=100, - verbosity=0, - tol=1E-8, - ishermitian=true, - issymmetric=true, - eager=true, - ) - exponentiate_kwargs = updater_kwargs.exponentiate_kwargs - exponentiate_kwargs = merge(default_exponentiate_kwargs, exponentiate_kwargs) #last collection has precedence - result, exp_info = exponentiate( - projected_operator![], region_kwargs.time_step, expanded_init; exponentiate_kwargs... - ) - - # - #ToDo: return truncation error and append to info - # truncate - # ToDo: refactor - # ToDo: either remove do_truncate, or make accessible as kwarg - # ToDo: think about more elaborate ways to trigger truncation, i.e. truncate only once maxdim is hit? - # ToDo: add truncation info - do_truncate = true - # simply return if we are not truncating - !do_truncate && return result, (; info=exp_info) - # check for cases where we don't want to try truncating - region = first(sweep_plan[which_region_update]) - typeof(region) <: NamedEdge && return result, (; info=exp_info) - (; maxdim, cutoff) = region_kwargs - region = only(region) - next_region = if which_region_update == length(sweep_plan) - nothing - else - first(sweep_plan[which_region_update + 1]) - end - previous_region = - which_region_update == 1 ? nothing : first(sweep_plan[which_region_update - 1]) - isnothing(next_region) && return result, (; info=exp_info) - !(typeof(next_region) <: NamedEdge) && return result, (; info=exp_info) - !(region == src(next_region) || region == dst(next_region)) && - return result, (; info=exp_info) - # actually truncate - left_inds = uniqueinds(state![], next_region) - U, S, V = svd( - result, - left_inds; - lefttags=tags(state![], next_region), - righttags=tags(state![], next_region), - maxdim, - cutoff, - ) - next_vertex = src(next_region) == region ? dst(next_region) : src(next_region) - state = copy(state![]) - #@show inds(V) - state[next_vertex] = state[next_vertex] * V - _nsites = (region isa AbstractEdge) ? 0 : length(region) #should be 1 - #@show _nsites - PH = copy(projected_operator![]) - PH = set_nsite(PH, 2) - PH = position(PH, state, [region, next_vertex]) - PH = set_nsite(PH, _nsites) - PH = position(PH, state, first(sweep_plan[which_region_update])) - state![] = state - projected_operator![] = PH - return U * S, (; info=exp_info) -end From cea57a070df120acdd64e4c93a45f177b5b635d5 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 17:06:35 -0400 Subject: [PATCH 12/19] Remove duplicate implementation of two-site expansion, leaving only commented out full expansion. --- .../subspace_expansions/subspace_expansion.jl | 249 ------------------ 1 file changed, 249 deletions(-) diff --git a/src/solvers/subspace_expansions/subspace_expansion.jl b/src/solvers/subspace_expansions/subspace_expansion.jl index 55fbce5e..2e85b40d 100644 --- a/src/solvers/subspace_expansions/subspace_expansion.jl +++ b/src/solvers/subspace_expansions/subspace_expansion.jl @@ -1,252 +1,3 @@ -function two_site_expansion_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - region_kwargs, - updater_kwargs, -) - (; maxdim, cutoff, time_step) = region_kwargs - #ToDo: handle timestep==Inf for DMRG case - default_updater_kwargs = (; - svd_func_expand=rsvd_iterative, - maxdim=(maxdim == typemax(Int) ? maxdim : Int(ceil((2.0^(1 / 3))) * maxdim)), - cutoff=isinf(time_step) ? cutoff : cutoff / abs(time_step), # ToDo verify that this is the correct way of scaling the cutoff - use_relative_cutoff=false, - use_absolute_cutoff=true, - ) - updater_kwargs = merge(default_updater_kwargs, updater_kwargs) - - # if on edge return without doing anything - region = first(sweep_plan[which_region_update]) - typeof(region) <: NamedEdge && return init, (;) - region = only(region) - #figure out next region, since we want to expand there - #ToDo account for non-integer indices into which_region_update - next_region = if which_region_update == length(sweep_plan) - nothing - else - first(sweep_plan[which_region_update + 1]) - end - previous_region = - which_region_update == 1 ? nothing : first(sweep_plan[which_region_update - 1]) - isnothing(next_region) && return init, (;) - !(typeof(next_region) <: NamedEdge) && return init, (;) - !(region == src(next_region) || region == dst(next_region)) && return init, (;) - next_vertex = src(next_region) == region ? dst(next_region) : src(next_region) - vp = region => next_vertex - - phi, has_changed = _two_site_expand_core( - init, region, vp; projected_operator!, state!, expand_dir=1, updater_kwargs... - ) - !has_changed && return init, (;) - return phi, (;) -end - -""" -Returns a function which applies Nullspace projector to an ITensor with matching indices without -explicitly constructing the projector as an ITensor or Network. -""" -function implicit_nullspace(A, linkind) - # only works when applied in the direction of the environment tensor, not the other way (due to use of ITensorNetworkMap which is not reversible) - outind = uniqueinds(A, linkind) - inind = outind #? - # ToDo: implement without ITensorNetworkMap - P = ITensors.ITensorNetworkMaps.ITensorNetworkMap( - [prime(dag(A), linkind), prime(A, linkind)], inind, outind - ) - return x::ITensor -> x - P(x) -end - -""" -Performs a local expansion using the two-site projected variance. -The input state should have orthogonality center on a single vertex (region), and phi0 is that site_tensors. -Expansion is performed on vertex that would be visited next (if next vertex is a neighbor of region). -""" -function _two_site_expand_core( - phi0, - region, - vertexpair::Pair; - projected_operator!, - state!, - expand_dir=1, - svd_func_expand, - cutoff, - maxdim, - use_relative_cutoff, - use_absolute_cutoff, -) - # preliminaries - theflux = flux(phi0) - svd_func = svd_func_expand - v1 = first(vertexpair) - v2 = last(vertexpair) - verts = [v1, v2] - n1, n2 = 1, 2 - PH = copy(projected_operator![]) - psi = copy(state![]) - # orthogonalize on the edge - next_edge = edgetype(psi)(vertexpair) - psi, phi = extract_local_tensor(psi, next_edge) - psis = map(n -> psi[n], verts) # extract local site tensors - # ToDo: remove following code lines unless we want to truncate before we expand? - ##this block is replaced by extract_local_tensor --- unless we want to truncate here this is the cleanest. - ##otherwise reinsert the explicit block with truncation args (would allow to remove it from higher-level expander) - #left_inds = uniqueinds(psis[n1], psis[n2]) - #U, S, V = svd(psis[n1], left_inds; lefttags=tags(commonind(psis[n1],psis[n2])), righttags=tags(commonind(psis[n1],psis[n2]))) - #psis[n1]= U - #psi[region]=U - #phi = S*V - - # don't expand if we are already at maxdim - old_linkdim = dim(commonind(first(psis), phi)) - linkinds = map(psi -> commonind(psi, phi), psis) - (old_linkdim >= maxdim) && return phi0, false - - # compute nullspace to the left and right - #@timeit_debug timer "nullvector" begin - nullVecs = map(zip(psis, linkinds)) do (psi, linkind) - #return nullspace(psi, linkind; atol=atol) - return implicit_nullspace(psi, linkind) - end - #end - - # update the projected operator - PH = set_nsite(PH, 2) - PH = position(PH, psi, verts) - - # build environments - g = underlying_graph(PH) - #@timeit_debug timer "build environments" begin - envs = map(zip(verts, psis)) do (n, psi) - return noprime( - mapreduce(*, [v => n for v in neighbors(g, n) if !(v ∈ verts)]; init=psi) do e - return PH.environments[NamedEdge(e)] - end * PH.H[n], - ) - end - #end - - # apply the projection into nullspace - envs = [last(nullVecs)(last(envs)), first(nullVecs)(first(envs))] - - # assemble ITensorNetworkMap #ToDo: do not rely on ITensorNetworkMap, simplify logic - ininds = uniqueinds(last(psis), phi) - outinds = uniqueinds(first(psis), phi) - cin = combiner(ininds) - cout = combiner(outinds) - envs = [cin * envs[1], cout * envs[2]] - envMap=[last(envs), phi* first(envs)] - #@show inds(envMap[1]) - #@show inds(envMap[2]) - - - - #envMap = ITensors.ITensorNetworkMaps.ITensorNetworkMap( - # [last(envs), phi, first(envs)], - # uniqueinds(inds(cout), outinds), - # uniqueinds(inds(cin), ininds), - #) - #envMapDag = adjoint(envMap) - - # factorize - #@timeit_debug timer "svd_func" begin - if svd_func == ITensorNetworks._svd_solve_normal - U, S, V = svd_func( - envMap, uniqueinds(inds(cout), outinds); maxdim=maxdim - old_linkdim, cutoff=cutoff - ) - elseif svd_func == ITensorNetworks.rsvd_iterative - - U, S, V = svd_func( - envMap, - uniqueinds(inds(cout), outinds); - maxdim=maxdim - old_linkdim, - cutoff=cutoff, - use_relative_cutoff=false, - use_absolute_cutoff=true, - ) - - #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, - #use_absolute_cutoff=true) #this one is for debugging in case we want to test the precontracted version - else - U, S, V = svd_func( - eltype(envMap), - envMap, - envMapDag, - uniqueinds(inds(cout), outinds); - flux=theflux, - maxdim=maxdim - old_linkdim, - cutoff=cutoff, - ) - end - #end - - # catch cases when we decompose a map of norm==0.0 - (isnothing(U) || iszero(norm(U))) && return phi0, false - #FIXME: somehow the svd funcs sometimes return empty ITensors instead of nothing, that should be caught in the SVD routines instead... - all(isempty.([U, S, V])) && return phi0, false #ToDo: do not rely on isempty here - - # uncombine indices on the non-link-indices - U *= dag(cout) - V *= dag(cin) - #@show inds(U) - #@show inds(V) - #@show inds(S) - - # direct sum the site tensors - #@timeit_debug timer "direct sum" begin - new_psis = map(zip(psis, [U, V])) do (psi, exp_basis) - return ITensors.directsum( - psi => commonind(psi, phi), - exp_basis => uniqueind(exp_basis, psi); - tags=tags(commonind(psi, phi)), - ) - end - #end - new_inds = [last(x) for x in new_psis] - new_psis = [first(x) for x in new_psis] - - # extract the expanded linkinds from expanded site tensors and create a zero-tensor - phi_indices = replace( - inds(phi), (commonind(phi, psis[n]) => dag(new_inds[n]) for n in 1:2)... - ) - if hasqns(phi) - new_phi = ITensor(eltype(phi), flux(phi), phi_indices...) - fill!(new_phi, 0.0) #ToDo: Check whether this is really necessary. - else - new_phi = ITensor(eltype(phi), phi_indices...) - end - - # set the new bond tensor elements from the old bond tensor - map(eachindex(phi)) do I - v = phi[I] - !iszero(v) && (return new_phi[I] = v) # I think this line errors without the fill! with zeros above - end - - # apply combiners on linkinds #ToDo: figure out why this is strictly necessary - combiners = map( - new_ind -> combiner(new_ind; tags=tags(new_ind), dir=dir(new_ind)), new_inds - ) - for (v, new_psi, C) in zip(verts, new_psis, combiners) - psi[v] = noprime(new_psi * C) - end - new_phi = dag(first(combiners)) * new_phi * dag(last(combiners)) - - # apply expanded bond tensor to site tensor and reset projected operator to region - psi[v1] *= new_phi - new_phi = psi[v1] - PH = set_nsite(PH, 1) - PH = position(PH, psi, [region]) - - projected_operator![] = PH - state![] = psi - return new_phi, true - ##body end -end - #= function _full_expand_core_vertex( PH, psi, phi, region, svd_func; direction, expand_dir=-1, expander_cache=Any[], maxdim, cutoff, atol=1e-8, kwargs..., From 0400446549bb43b0aad38abb7e0efde72b4fdea7 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 17:08:18 -0400 Subject: [PATCH 13/19] Update includes. --- src/ITensorNetworks.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 5beb14f5..b5491359 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -42,6 +42,7 @@ include("solvers/local_solvers/exponentiate.jl") include("solvers/local_solvers/dmrg_x.jl") include("solvers/local_solvers/contract.jl") include("solvers/local_solvers/linsolve.jl") +include("solvers/local_solvers/expand.jl") include("treetensornetworks/abstracttreetensornetwork.jl") include("treetensornetworks/ttn.jl") include("treetensornetworks/opsum_to_ttn.jl") @@ -53,10 +54,9 @@ include("solvers/solver_utils.jl") include("solvers/defaults.jl") include("solvers/insert/insert.jl") include("solvers/extract/extract.jl") -include("solvers/subspace_expansions/subspace_expansion.jl") include("solvers/subspace_expansions/linalg/standard_svd.jl") -include("solvers/subspace_expansions/linalg/rsvd_aux.jl") -include("solvers/subspace_expansions/linalg/rsvd_linalg.jl") +#include("solvers/subspace_expansions/linalg/rsvd_aux.jl") +#include("solvers/subspace_expansions/linalg/rsvd_linalg.jl") include("solvers/alternating_update/alternating_update.jl") include("solvers/alternating_update/region_update.jl") include("solvers/tdvp.jl") From 0c5bb69af55617fb0d914ed8e35ed68da5641435 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 17:38:38 -0400 Subject: [PATCH 14/19] Cleanup two-site expansion, unify svd wrapper signature. --- src/solvers/local_solvers/expand.jl | 70 ++++++------------- .../linalg/standard_svd.jl | 11 +-- 2 files changed, 29 insertions(+), 52 deletions(-) diff --git a/src/solvers/local_solvers/expand.jl b/src/solvers/local_solvers/expand.jl index 5b8a0f89..575b242b 100644 --- a/src/solvers/local_solvers/expand.jl +++ b/src/solvers/local_solvers/expand.jl @@ -19,14 +19,13 @@ function two_site_expansion_updater( use_relative_cutoff=false, use_absolute_cutoff=true, ) - #return init, (;) maxdim=maxdim_func(maxdim;internal_kwargs...) cutoff=cutoff_func(cutoff;internal_kwargs...) region = first(sweep_plan[which_region_update]) typeof(region) <: NamedEdge && return init, (;) region = only(region) - #figure out next region, since we want to expand there - #ToDo account for non-integer indices into which_region_update + # figure out next region, since we want to expand there + # ToDo: account for non-integer indices into which_region_update next_region = if which_region_update == length(sweep_plan) nothing else @@ -38,10 +37,11 @@ function two_site_expansion_updater( !(typeof(next_region) <: NamedEdge) && return init, (;) !(region == src(next_region) || region == dst(next_region)) && return init, (;) next_vertex = src(next_region) == region ? dst(next_region) : src(next_region) - vp = region => next_vertex phi, has_changed = _two_site_expand_core( - init, region, vp; projected_operator!, state!, expand_dir=1, + init, region, region => next_vertex; + projected_operator!, + state!, svd_func_expand, maxdim, cutoff, @@ -78,7 +78,6 @@ function _two_site_expand_core( vertexpair::Pair; projected_operator!, state!, - expand_dir=1, svd_func_expand, cutoff, maxdim, @@ -87,25 +86,19 @@ function _two_site_expand_core( ) # preliminaries theflux = flux(phi0) - svd_func = svd_func_expand v1 = first(vertexpair) v2 = last(vertexpair) verts = [v1, v2] - n1, n2 = 1, 2 PH = copy(projected_operator![]) - psi = copy(state![]) + state = copy(state![]) + old_linkdim = dim(commonind(state[v1],state[v2])) + (old_linkdim >= maxdim) && return phi0, false # orthogonalize on the edge - next_edge = edgetype(psi)(vertexpair) #update of environment not strictly necessary here - psi, PH, phi = default_extracter(psi, PH, next_edge, v1;internal_kwargs=(;))#default_extracter(psi,next_edge)#extract_local_tensor(psi, next_edge) - - psis = map(n -> psi[n], verts) # extract local site tensors - - #if we are already past the maximum allowed bond dimension, return the input - #ToDo: move up so we don't do any orthogonalization - old_linkdim = dim(commonind(first(psis), phi)) + state, PH, phi = default_extracter(state, PH, edgetype(state)(vertexpair), v1; + internal_kwargs=(;)) + psis = map(n -> state[n], verts) # extract local site tensors linkinds = map(psi -> commonind(psi, phi), psis) - (old_linkdim >= maxdim) && return phi0, false # compute nullspace to the left and right #@timeit_debug timer "nullvector" begin @@ -118,7 +111,7 @@ function _two_site_expand_core( # update the projected operator #ToDo: may not be necessary if we use the extracter anyway PH = set_nsite(PH, 2) - PH = position(PH, psi, verts) + PH = position(PH, state, verts) # build environments g = underlying_graph(PH) @@ -144,34 +137,14 @@ function _two_site_expand_core( # factorize #FIXME: remove specialization on svd_func, define interface for svd_funcs #@timeit_debug timer "svd_func" begin - if svd_func == ITensorNetworks._svd_solve_normal - U, S, V = svd_func( - envMap, uniqueinds(inds(cout), outinds); maxdim=maxdim - old_linkdim, cutoff=cutoff - ) - elseif svd_func == ITensorNetworks.rsvd_iterative - - U, S, V = svd_func( - envMap, - uniqueinds(inds(cout), outinds); - maxdim=maxdim - old_linkdim, - cutoff=cutoff, - use_relative_cutoff=false, - use_absolute_cutoff=true, - ) - - #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, - #use_absolute_cutoff=true) #this one is for debugging in case we want to test the precontracted version - else - U, S, V = svd_func( - eltype(envMap), + U, S, V = svd_func_expand( envMap, - envMapDag, uniqueinds(inds(cout), outinds); - flux=theflux, maxdim=maxdim - old_linkdim, cutoff=cutoff, + use_relative_cutoff, + use_absolute_cutoff, ) -end #end # catch cases when we decompose a map of norm==0.0 @@ -202,7 +175,8 @@ end ) if hasqns(phi) new_phi = ITensor(eltype(phi), flux(phi), phi_indices...) - fill!(new_phi, 0.0) #ToDo: Check whether this is really necessary. + #ToDo: Remove? This used to be necessary, but may have been fixed with bugfixes in ITensor + #fill!(new_phi, 0.0) else new_phi = ITensor(eltype(phi), phi_indices...) end @@ -218,16 +192,16 @@ end new_ind -> combiner(new_ind; tags=tags(new_ind), dir=dir(new_ind)), new_inds ) for (v, new_psi, C) in zip(verts, new_psis, combiners) - psi[v] = noprime(new_psi * C) + state[v] = noprime(new_psi * C) end new_phi = dag(first(combiners)) * new_phi * dag(last(combiners)) # apply expanded bond tensor to site tensor and reset projected operator to region - psi[v1] *= new_phi - new_phi = psi[v1] + state[v1] *= new_phi + new_phi = state[v1] PH = set_nsite(PH, 1) - PH = position(PH, psi, [region]) + PH = position(PH, state, [region]) projected_operator![] = PH - state![] = psi + state![] = state return new_phi, true end \ No newline at end of file diff --git a/src/solvers/subspace_expansions/linalg/standard_svd.jl b/src/solvers/subspace_expansions/linalg/standard_svd.jl index dbabecea..3b8b66c0 100644 --- a/src/solvers/subspace_expansions/linalg/standard_svd.jl +++ b/src/solvers/subspace_expansions/linalg/standard_svd.jl @@ -1,10 +1,13 @@ #ToDo: pass use_*_cutoff as kwargs -function _svd_solve_normal(envMap, left_ind; maxdim, cutoff) +function _svd_solve_normal(envMap, left_ind; + maxdim, + cutoff, + use_relative_cutoff, + use_absolute_cutoff) M = ITensors.ITensorNetworkMaps.contract(envMap) - # ToDo: infer eltype from envMap - norm(M) ≤ eps(Float64) && return nothing, nothing, nothing + norm(M) ≤ eps(real(eltype(M))) && return nothing, nothing, nothing U, S, V = svd( - M, left_ind; maxdim, cutoff=cutoff, use_relative_cutoff=false, use_absolute_cutoff=true + M, left_ind; maxdim, cutoff, use_relative_cutoff, use_absolute_cutoff ) vals = diag(array(S)) (length(vals) == 1 && vals[1]^2 ≤ cutoff) && return nothing, nothing, nothing From b7d08b69e92064ae835cad762c86830aa64045c0 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 17:38:58 -0400 Subject: [PATCH 15/19] Specify reasonable default for scaling maxdim for subspace expansion. --- src/solvers/defaults.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/defaults.jl b/src/solvers/defaults.jl index 13dbafe5..769285f7 100644 --- a/src/solvers/defaults.jl +++ b/src/solvers/defaults.jl @@ -8,7 +8,7 @@ default_inserter() = default_inserter default_checkdone() = (; kws...) -> false default_transform_operator() = nothing default_scale_cutoff_by_timestep(cutoff;time_step, kwargs...) = isinf(time_step) ? cutoff : cutoff / abs(time_step) - +default_scale_maxdim(maxdim;kwargs...) = Int(ceil(2^(1.0/3.0) * maxdim)) function default_region_printer(; cutoff, maxdim, From 2896eb161566ba87a2ea6feea1e358b39da7509c Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Apr 2024 17:39:44 -0400 Subject: [PATCH 16/19] Add example test for subspace expansion, to be improved later. --- .../notatest_tdvp_subspace_minimal.jl | 133 ++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 test/test_treetensornetworks/test_solvers/notatest_tdvp_subspace_minimal.jl diff --git a/test/test_treetensornetworks/test_solvers/notatest_tdvp_subspace_minimal.jl b/test/test_treetensornetworks/test_solvers/notatest_tdvp_subspace_minimal.jl new file mode 100644 index 00000000..a273dedf --- /dev/null +++ b/test/test_treetensornetworks/test_solvers/notatest_tdvp_subspace_minimal.jl @@ -0,0 +1,133 @@ +using ITensors +using ITensorNetworks +using ITensorNetworks: ModelHamiltonians.tight_binding +using ITensorNetworks: exponentiate_updater, _svd_solve_normal, ttn, tdvp +using ITensorNetworks: local_expand_and_exponentiate_updater, vertices +using ITensorNetworks: two_site_expansion_updater, extract_and_truncate, compose_updaters + +#using KrylovKit: exponentiate +using Observers +using Random +using Test +using ITensorGaussianMPS: hopping_operator +using LinearAlgebra: exp, I + +""" +Propagate a Gaussian density matrix `c` with a noninteracting Hamiltonian `h` up to time `tau`. +""" +function propagate_noninteracting(h::AbstractMatrix, c::AbstractMatrix, tau::Number) + U = exp(-im * tau * h) + return U * c * conj(transpose(U)) +end + +@testset "MPS TDVP" begin + @testset "2s vs 1s + local subspace" begin + ITensors.enable_auto_fermion() + Random.seed!(1234) + cutoff = 1e-14 + N = 20 + D = 8 + t = 1.0 + tp = 0.0 + g = ITensorNetworks.named_path_graph(N) + s = siteinds("Fermion", g; conserve_qns=true) + os = tight_binding(g; t, tp) + H = ttn(os, s) + hmat = hopping_operator(os) + # get the exact result + tf = 1.0 + taus = LinRange(0, tf, 2) + #init=I(N)[:,StatsBase.sample(1:N, div(N,2); replace=false)] + init = I(N)[:, 1:(div(N, 2))] #domain wall initial condition + states = i -> i <= div(N, 2) ? "Occ" : "Emp" + init = conj(init) * transpose(init) + #res=[] + res = zeros(N, length(taus)) + for (i, tau) in enumerate(taus) + res[:, i] = real.(diag(propagate_noninteracting(hmat, init, tau))) #densities only + #plot!(1:N, real.(diag(last(res)))) + end + psi = ttn(ITensorNetwork(ComplexF64,states,s;link_space=1)) + psi0 = deepcopy(psi) + function my_sweep_printer(;which_sweep,state,kwargs...) + m = maximum(ITensorNetworks.edge_data(ITensorNetworks.linkdims(state))) + println("Sweep ", which_sweep,", maximum bond dimension: ",m) + end + ###compose_updater + dt = 0.05 + tdvp_kwargs=( ; + updater=compose_updaters(;expander=two_site_expansion_updater,solver=exponentiate_updater), + updater_kwargs=(; + expander_kwargs=(; + cutoff=cutoff, + maxdim=D, + svd_func_expand=ITensorNetworks._svd_solve_normal, + maxdim_func=ITensorNetworks.default_scale_maxdim, + cutoff_func=ITensorNetworks.default_scale_cutoff_by_timestep, + ), + solver_kwargs=(; + tol=1E-8 + ) + ), + extracter=extract_and_truncate, + extracter_kwargs=(;maxdim=D,cutoff=cutoff), + sweep_printer=my_sweep_printer, + time_step=-im * dt, + reverse_step=true, + order=2, + normalize=true, + maxdim=D, + cutoff=cutoff, + outputlevel=1, + ) + + #= + tdvp_kwargs = ( + time_step=-im * dt, + reverse_step=true, + order=2, + normalize=true, + maxdim=D, + cutoff=cutoff, + outputlevel=1, + updater_kwargs=(; + expand_kwargs=(; + cutoff=cutoff, + maxdim=D, + svd_func_expand=ITensorNetworks._svd_solve_normal + ), + exponentiate_kwargs=(;), + ), + )#,exponentiate_kwargs=(;tol=1e-8))) + =# + success = false + psife = nothing + #while !success + # try + psife = tdvp( + H, + -1im * tf, + psi; + nsites=1, + #updater=ITensorNetworks.local_expand_and_exponentiate_updater, + tdvp_kwargs..., + ) + # success=true + # catch + # println("trying again") + # end + # + #end + + @test inner(psi0, psi) ≈ 1 atol = 1e-12 + #@show maxlinkdim(psife) + #mag_2s=expect("N",psif) + mag_exp = expect("N", psife) + # @test real.([mag_2s[v] for v in vertices(psif)]) ≈ res[:,end] atol=1e-3 + @show maximum(abs.(real.([mag_exp[v] for v in vertices(psife)]) .- res[:, end])) + @test all( + i -> i < 5e-3, abs.(real.([mag_exp[v] for v in vertices(psife)]) .- res[:, end]) + ) + end +end +nothing From 732506dae84e7f0962119868f08200d6cdb2a1a3 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Tue, 23 Apr 2024 15:13:22 -0400 Subject: [PATCH 17/19] Remove code for global subspace expansion. --- .../subspace_expansions/subspace_expansion.jl | 149 ------------------ 1 file changed, 149 deletions(-) delete mode 100644 src/solvers/subspace_expansions/subspace_expansion.jl diff --git a/src/solvers/subspace_expansions/subspace_expansion.jl b/src/solvers/subspace_expansions/subspace_expansion.jl deleted file mode 100644 index 2e85b40d..00000000 --- a/src/solvers/subspace_expansions/subspace_expansion.jl +++ /dev/null @@ -1,149 +0,0 @@ -#= -function _full_expand_core_vertex( - PH, psi, phi, region, svd_func; direction, expand_dir=-1, expander_cache=Any[], maxdim, cutoff, atol=1e-8, kwargs..., -) ###ToDo: adapt to current interface, broken as of now. - #@show cutoff - #enforce expand_dir in the tested direction - @assert expand_dir==-1 - #println("in full expand") - ### only on edges - #(typeof(region)!=NamedEdge{Int}) && return psi, phi, PH - ### only on verts - (typeof(region)==NamedEdge{Int}) && return psi, phi, PH - ## kind of hacky - only works for mps. More general? - n1 = first(region) - theflux=flux(psi[n1]) - #expand_dir=-1 - if direction == 1 - n2 = expand_dir == 1 ? n1 - 1 : n1+1 - else - n2 = expand_dir == 1 ? n1 + 1 : n1-1 - end - - (n2 < 1 || n2 > length(vertices(psi))) && return psi,phi,PH - neutralflux=flux(psi[n2]) - - verts = [n1,n2] - - if isempty(expander_cache) - @warn("building environment of H^2 from scratch!") - - g = underlying_graph(PH.H) - H = vertex_data(data_graph(PH.H)) - - H_dag = swapprime.(prime.(dag.(H)), 1,2, tags = "Site") - H2_vd= replaceprime.(map(*, H, H_dag), 2 => 1) - H2_ed = edge_data(data_graph(PH.H)) - - H2 = TTN(ITensorNetwork(DataGraph(g, H2_vd, H2_ed)), PH.H.ortho_center) - PH2 = ProjTTN(H2) - PH2 = set_nsite(PH2, 2) - - push!(expander_cache, PH2) - end - - PH2 = expander_cache[1] - n1 = verts[2] - n2 = verts[1] - - ###do we really need to make the environments two-site? - ###look into how ProjTTN works - PH2 = position(PH2, psi, [n1,n2]) - PH = position(PH, psi, [n1,n2]) - - psi1 = psi[n1] - psi2 = psi[n2] #phi - old_linkdim = dim(commonind(psi1, psi2)) - - # don't expand if we are already at maxdim - ## make more transparent that this is not the normal maxdim arg but maxdim_expand - ## when would this even happen? - #@show old_linkdim, maxdim - - old_linkdim >= maxdim && return psi, phi, PH - #@show "expandin", maxdim, old_linkdim, maxdim-old_linkdim - # compute nullspace to the left and right - linkind_l = commonind(psi1, psi2) - nullVec = implicit_nullspace(psi1, linkind_l) - - # if nullspace is empty (happen's for product states with QNs) - ## ToDo implement norm or equivalent for a generic LinearMap (i guess via sampling a random vector) - #norm(nullVec) == 0.0 && return psi, phi, PH - - ## compute both environments - g = underlying_graph(PH) - - @timeit_debug timer "build environments" begin - env1 = noprime(mapreduce(*, [v => n1 for v in neighbors(g, n1) if v != n2]; init = psi1) do e - return PH.environments[NamedEdge(e)] - end *PH.H[n1] - ) - env2p = noprime(mapreduce(*, [v => n2 for v in neighbors(g, n2) if v != n1]; init = psi2) do e - return PH.environments[NamedEdge(e)] - end *PH.H[n2] - ) - - env2 = mapreduce(*, [v => n2 for v in neighbors(g, n2) if v != n1]; init = psi2) do e - return PH2.environments[NamedEdge(e)] - end * PH2.H[n2]*prime(dag(psi2)) - end - - env1=nullVec(env1) - - outinds = uniqueinds(psi1,psi2) - ininds = dag.(outinds) - cout=combiner(outinds) - env1 *= cout - env2p *= prime(dag(psi[n2]),commonind(dag(psi[n2]),dag(psi[n1]))) - env2p2= replaceprime(env2p * replaceprime(dag(env2p),0=>2),2=>1) - - envMap = ITensors.ITensorNetworkMaps.ITensorNetworkMap([prime(dag(env1)),(env2-env2p2),env1] , prime(dag(uniqueinds(cout,outinds))), uniqueinds(cout,outinds)) - envMapDag=adjoint(envMap) - # svd-decomposition - @timeit_debug timer "svd_func" begin - if svd_func==ITensorNetworks._svd_solve_normal - U,S,_ = svd_func(envMap, uniqueinds(inds(cout),outinds); maxdim=maxdim-old_linkdim, cutoff=cutoff) - elseif svd_func==ITensorNetworks.rsvd_iterative - #@show theflux - envMap=transpose(envMap) - U,S,_ = svd_func(eltype(first(envMap.itensors)),envMap,ITensors.ITensorNetworkMaps.input_inds(envMap);theflux=neutralflux, maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, - use_absolute_cutoff=true) - #U,S,V = svd_func(contract(envMap),uniqueinds(inds(cout),outinds);theflux=theflux, maxdim=maxdim-old_linkdim, cutoff=cutoff, use_relative_cutoff=false, - #use_absolute_cutoff=true) - else - U,S,_= svd_func(eltype(envMap),envMap,envMapDag,uniqueinds(cout,outinds); flux=neutralflux, maxdim=maxdim-old_linkdim, cutoff=cutoff) - end - end - isnothing(U) && return psi,phi,PH - ###FIXME: somehow the svd funcs sometimes return empty ITensors instead of nothing, that should be caught in the SVD routines instead... - all(isempty.([U,S])) && return psi, phi, PH - @assert dim(commonind(U, S)) ≤ maxdim-old_linkdim - nullVec = dag(cout)*U - new_psi1, new_ind1 = ITensors.directsum( - psi1 => uniqueinds(psi1, nullVec), nullVec => uniqueinds(nullVec, psi1); tags=(tags(commonind(psi1,phi)),) - ) - new_ind1 = new_ind1[1] - @assert dim(new_ind1) <= maxdim - - Cl = combiner(new_ind1; tags=tags(new_ind1), dir=dir(new_ind1)) - - phi_indices = replace(inds(phi), commonind(phi,psi1) => dag(new_ind1)) - - if hasqns(phi) - new_phi=ITensor(eltype(phi),flux(phi), phi_indices...) - fill!(new_phi,0.0) - else - new_phi = ITensor(eltype(phi), phi_indices...) - end - - map(eachindex(phi)) do I - v = phi[I] - !iszero(v) && (return new_phi[I]=v) - end - - psi[n1] = noprime(new_psi1*Cl) - new_phi = dag(Cl)*new_phi - - return psi, new_phi, PH -end -=# From c6235cd9565b49c06adacdc32bb9bdfc22056b8b Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Tue, 23 Apr 2024 15:15:43 -0400 Subject: [PATCH 18/19] Remove outdated comments. --- src/solvers/extract/extract.jl | 6 +----- src/solvers/local_solvers/expand.jl | 3 +-- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/src/solvers/extract/extract.jl b/src/solvers/extract/extract.jl index abce0829..98506ec5 100644 --- a/src/solvers/extract/extract.jl +++ b/src/solvers/extract/extract.jl @@ -26,10 +26,6 @@ function default_extracter(state, projected_operator, region, ortho; internal_kw end function extract_and_truncate(state,projected_operator, region, ortho; maxdim=nothing, cutoff=nothing,internal_kwargs) - #ToDo: move default into function signature? - #svd_kwargs = NamedTuple() - #svd_kwargs = isnothing(maxdim) ? svd_kwargs : merge(svd_kwargs,(;maxdim)) - #svd_kwargs = isnothing(cutoff) ? svd_kwargs : merge(svd_kwargs,(;cutoff)) svd_kwargs= (; (isnothing(maxdim) ? (;) : (;maxdim))..., (isnothing(cutoff) ? (;) : (;cutoff))... @@ -49,4 +45,4 @@ function extract_and_truncate(state,projected_operator, region, ortho; maxdim=no end projected_operator = position(projected_operator, state, region) return state, projected_operator, local_tensor -end \ No newline at end of file +end diff --git a/src/solvers/local_solvers/expand.jl b/src/solvers/local_solvers/expand.jl index 575b242b..5228768a 100644 --- a/src/solvers/local_solvers/expand.jl +++ b/src/solvers/local_solvers/expand.jl @@ -135,7 +135,6 @@ function _two_site_expand_core( envMap=[cout * envs[1], phi* (cin * envs[2])] # factorize - #FIXME: remove specialization on svd_func, define interface for svd_funcs #@timeit_debug timer "svd_func" begin U, S, V = svd_func_expand( envMap, @@ -204,4 +203,4 @@ function _two_site_expand_core( projected_operator![] = PH state![] = state return new_phi, true -end \ No newline at end of file +end From 6f4993b813a47d6d89ac3febd14273ba77ab011b Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sun, 28 Jul 2024 10:05:38 -0400 Subject: [PATCH 19/19] Remove rsvd related files, to be included in separate MR. --- src/ITensorNetworks.jl | 2 - .../subspace_expansions/linalg/rsvd_aux.jl | 251 ----------------- .../subspace_expansions/linalg/rsvd_linalg.jl | 259 ------------------ 3 files changed, 512 deletions(-) delete mode 100644 src/solvers/subspace_expansions/linalg/rsvd_aux.jl delete mode 100644 src/solvers/subspace_expansions/linalg/rsvd_linalg.jl diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 2e861939..a864d7c1 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -53,8 +53,6 @@ include("solvers/defaults.jl") include("solvers/insert/insert.jl") include("solvers/extract/extract.jl") include("solvers/subspace_expansions/linalg/standard_svd.jl") -#include("solvers/subspace_expansions/linalg/rsvd_aux.jl") -#include("solvers/subspace_expansions/linalg/rsvd_linalg.jl") include("solvers/alternating_update/alternating_update.jl") include("solvers/alternating_update/region_update.jl") include("solvers/tdvp.jl") diff --git a/src/solvers/subspace_expansions/linalg/rsvd_aux.jl b/src/solvers/subspace_expansions/linalg/rsvd_aux.jl deleted file mode 100644 index cedc11c1..00000000 --- a/src/solvers/subspace_expansions/linalg/rsvd_aux.jl +++ /dev/null @@ -1,251 +0,0 @@ -function get_column_space(A::Vector{<:ITensor}, lc::Index,rc::Index) - #gets non-zero blocks in rc by sticking in lc and contracting through - viable_sectors=Vector{Pair{QN,Int64}} - for s in space(lc) - qn = first(s) - trial=randomITensor(flux(dag(qn)),lc) - adtrial=foldl(A;init=trial) - nnzblocks(adtrial)==0 && continue - thesector=only(space(only(inds(adtrial)))) - push!(viable_sectors, thesector) - end - return viable_sectors -end - - -function build_guess_matrix( - eltype::Type{<:Number}, ind, sectors::Union{Nothing,Vector{Pair{QN,Int64}}}, n::Int, p::Int - ) - if hasqns(ind) - aux_spaces = Pair{QN,Int64}[] - for s in sectors - thedim = last(s) - qn = first(s) - en = min(n + p, thedim) - push!(aux_spaces, Pair(qn, en)) - end - aux_ind = Index(aux_spaces; dir=dir(ind)) - try - M = randomITensor(eltype, dag(ind), aux_ind) #defaults to zero flux - # @show theflux, flux(M) - catch e - @show e - @show aux_ind - @show ind - error("stopping here something is wrong") - end - @assert nnzblocks(M) != 0 - else - thedim = dim(ind) - en = min(n + p, thedim) - #ep = max(0,en-n) - aux_ind = Index(en) - M = randomITensor(eltype, ind, aux_ind) - end - #@show M - return M -end -qns(t::ITensor) = qns(collect(eachnzblock(t)), t) -function qns(bs::Vector{Block{n}}, t::ITensor) where {n} - return [qns(b, t) for b in bs] -end - -function qns(b::Block{n}, t::ITensor) where {n} - theqns = QN[] - for i in 1:order(t) - theb = getindex(b, i) - theind = inds(t)[i] - push!(theqns, ITensors.qn(space(theind), Int(theb))) - end - return theqns -end - -function build_guess_matrix( - eltype::Type{<:Number}, ind, sectors::Union{Nothing,Vector{Pair{QN,Int64}}}, ndict::Dict; theflux=nothing, auxdir=dir(ind) -) - if hasqns(ind) - #translate_qns = get_qn_dict(ind, theflux; auxdir) - aux_spaces = Pair{QN,Int64}[] - #@show first.(space(ind)) - for s in sectors - thedim = last(s) - qn = first(s) - en = min(ndict[qn], thedim) - push!(aux_spaces, Pair(qn, en)) - end - aux_ind = Index(aux_spaces; dir=dir(ind)) - #M=randomITensor(eltype,-theflux,dag(ind),aux_ind) - #@show aux_ind - try - M = randomITensor(eltype, dag(ind), aux_ind) - catch e - @show e - @show aux_ind - @show ind - error("stopping here something is wrong") - end - @assert nnzblocks(M) != 0 - else - thedim = dim(ind) - en = min(ndict[space(ind)], thedim) - aux_ind = Index(en) - M = randomITensor(eltype, ind, aux_ind) - end - return M -end - -function init_guess_sizes(cind, sectors::Union{Nothing,Vector{Pair{QN,Int64}}}, n::Int, rule; theflux=nothing, auxdir=dir(cind)) - if hasqns(cind) - ndict = Dict{QN,Int64}() - for s in sectors - thedim = last(s) - qn = first(s) - ndict[qn] = min(rule(n), thedim) - end - else - @assert sectors==nothing - thedim = dim(cind) - ndict = Dict{Int64,Int64}() - ndict[thedim] = min(thedim, rule(n)) - end - return ndict -end - -function increment_guess_sizes(ndict::Dict{QN,Int64}, n_inc::Int, rule) - for key in keys(ndict) - #thedim=last(key) - ndict[key] = ndict[key] + rule(n_inc) - end - return ndict -end - -function increment_guess_sizes(ndict::Dict{Int64,Int64}, n_inc::Int, rule) - for key in keys(ndict) - #thedim=key - ndict[key] = ndict[key] + rule(n_inc) - end - return ndict -end - -function approx_the_same(o, n; abs_eps=1e-12, rel_eps=1e-6) - absdev = abs.(o .- n) - reldev = abs.((o .- n) ./ n) - abs_conv = absdev .< abs_eps - rel_conv = reldev .< rel_eps - return all(abs_conv .|| rel_conv) -end - -function is_converged_block(o, n; svd_kwargs...) - maxdim = get(svd_kwargs, :maxdim, Inf) - if length(o) != length(n) - return false - else - r = min(maxdim, length(o)) - #ToDo: Pass kwargs? - return approx_the_same(o[1:r], n[1:r]) - end -end - -function is_converged!(ndict, old_fact, new_fact; n_inc=1, has_qns=true, svd_kwargs...) - oS = old_fact.S - nS = new_fact.S - theflux = flux(nS) - oldflux = flux(oS) - if has_qns - if oldflux == nothing || theflux == nothing - if norm(oS) == 0.0 && norm(nS) == 0.0 - return true - else - return false - end - else - try - @assert theflux == flux(oS) - catch e - @show e - @show theflux - @show oldflux - error("Somehow the fluxes are not matching here! Exiting") - end - end - else - ###not entirely sure if this is legal for empty factorization - if norm(oS) == 0.0 - if norm(nS) == 0.0 - return true - else - return false - end - end - end - - maxdim = get(svd_kwargs, :maxdim, Inf) - os = sort(storage(oS); rev=true) - ns = sort(storage(nS); rev=true) - if length(os) >= maxdim && length(ns) >= maxdim - conv_global = approx_the_same(os[1:maxdim], ns[1:maxdim]) - elseif length(os) != length(ns) - conv_global = false - else - r = length(ns) - conv_global = approx_the_same(os[1:r], ns[1:r]) - end - if !hasqns(oS) - if conv_global == false - #ndict has only one key - ndict[first(keys(ndict))] *= 2 - end - return conv_global - end - conv_bool_total = true - ##a lot of this would be more convenient with ITensor internal_inds_space - ##ToDo: refactor, but not entirely trivial because functionality is implemented on the level of QNIndex, not a set of QNIndices - ##e.g. it is cumbersome to query the collection of QNs associated with a Block{n} of an ITensor with n>1 - soS = space(inds(oS)[1]) - snS = space(inds(nS)[1]) - qns = union(ITensors.qn.(soS), ITensors.qn.(snS)) - - oblocks = eachnzblock(oS) - oblockdict = Int.(getindex.(oblocks, 1)) - oqnindtoblock = Dict(collect(values(oblockdict)) .=> collect(keys(oblockdict))) - - nblocks = eachnzblock(nS) - nblockdict = Int.(getindex.(nblocks, 1)) - nqnindtoblock = Dict(collect(values(nblockdict)) .=> collect(keys(nblockdict))) - - for qn in qns - if qn in ITensors.qn.(snS) && qn in ITensors.qn.(soS) - oqnind = findfirst((first.(soS)) .== [qn]) - nqnind = findfirst((first.(snS)) .== (qn,)) - oblock = oqnindtoblock[oqnind] - nblock = nqnindtoblock[nqnind] - #oblock=ITensors.block(first,inds(oS)[1],qn) - #nblock=ITensors.block(first,inds(nS)[1],qn) - - #make sure blocks are the same QN when we compare them - #@assert first(soS[oqnind])==first(snS[nqnind])# - ovals = diag(oS[oblock]) - nvals = diag(nS[nblock]) - conv_bool = is_converged_block(collect(ovals), collect(nvals); svd_kwargs...) - else - conv_bool = false - end - if conv_bool == false - ndict[qn] *= 2 - end - conv_bool_total *= conv_bool - end - if conv_bool_total == true - @assert conv_global == true - else - if conv_global == true - println( - "Subspace expansion, rand. svd: singular vals converged globally, but may not be optimal, doing another iteration", - ) - end - end - return conv_bool_total::Bool -end - - - diff --git a/src/solvers/subspace_expansions/linalg/rsvd_linalg.jl b/src/solvers/subspace_expansions/linalg/rsvd_linalg.jl deleted file mode 100644 index d53d2030..00000000 --- a/src/solvers/subspace_expansions/linalg/rsvd_linalg.jl +++ /dev/null @@ -1,259 +0,0 @@ -function rsvd_iterative( - T, - A::ITensors.ITensorNetworkMaps.ITensorNetworkMap, - linds::Vector{<:Index}; - theflux=nothing, - svd_kwargs..., - ) - - #translate from in/out to l/r logic - ininds = ITensors.ITensorNetworkMaps.input_inds(A) - outinds = ITensors.ITensorNetworkMaps.output_inds(A) - @assert linds == ininds ##FIXME: this is specific to the way we wrote the subspace expansion, should be fixed in another iteration - rinds = outinds - - CL = combiner(linds...) - CR = combiner(rinds...) - cL = uniqueind(inds(CL), linds) - cR = uniqueind(inds(CR), rinds) - - l = CL * A.itensors[1] - r = A.itensors[end] * CR - - if length(A.itensors) !== 2 - AC = ITensors.ITensorNetworkMaps.ITensorNetworkMap( - [l, A.itensors[2:(length(A.itensors) - 1)]..., r], - commoninds(l, CL), - commoninds(r, CR), - ) - else - AC = ITensors.ITensorNetworkMaps.ITensorNetworkMap( - [l, r], commoninds(l, CL), commoninds(r, CR) - ) - end - ###this initializer part is still a bit ugly - n_init = 1 - p_rule(n) = 2 * n - ndict2 = init_guess_sizes(cR, n_init, p_rule; theflux=theflux) - ndict = init_guess_sizes(cL, n_init, p_rule; theflux=theflux) - if hasqns(ininds) - ndict = merge(ndict, ndict2) - else - ndict = ndict2 - end - M = build_guess_matrix(T, cR, ndict; theflux=theflux) - fact, Q = rsvd_core(AC, M; svd_kwargs...) - n_inc = 2 - ndict = increment_guess_sizes(ndict, n_inc, p_rule) - new_fact = deepcopy(fact) - while true - M = build_guess_matrix(T, cR, ndict; theflux=theflux) - new_fact, Q = rsvd_core(AC, M; svd_kwargs...) - if is_converged!(ndict, fact, new_fact; n_inc, has_qns=hasqns(ininds), svd_kwargs...) - break - else - fact = new_fact - end - end - vals = diag(array(new_fact.S)) - (length(vals) == 1 && vals[1]^2 ≤ get(svd_kwargs, :cutoff, 0.0)) && - return nothing, nothing, nothing - return dag(CL) * Q * new_fact.U, new_fact.S, new_fact.V * dag(CR) - end - - function rsvd_iterative(A::ITensor, linds::Vector{<:Index}; svd_kwargs...) - rinds = uniqueinds(A, linds) - CL = combiner(linds) - CR = combiner(rinds) - AC = CL * A * CR - cL = combinedind(CL) - cR = combinedind(CR) - if inds(AC) != (cL, cR) - AC = permute(AC, cL, cR) - end - n_init = 1 - p_rule(n) = 2 * n - iszero(norm(AC)) && return nothing, nothing, nothing - #@show flux(AC) - ndict2 = init_guess_sizes(cR, n_init, p_rule; theflux=hasqns(AC) ? flux(AC) : nothing) - ndict = init_guess_sizes(cL, n_init, p_rule; theflux=hasqns(AC) ? flux(AC) : nothing) - ndict = merge(ndict, ndict2) - M = build_guess_matrix(eltype(AC), cR, ndict; theflux=hasqns(AC) ? flux(AC) : nothing) - fact, Q = rsvd_core(AC, M; svd_kwargs...) - n_inc = 1 - ndict = increment_guess_sizes(ndict, n_inc, p_rule) - new_fact = deepcopy(fact) - while true - M = build_guess_matrix(eltype(AC), cR, ndict; theflux=hasqns(AC) ? flux(AC) : nothing) - new_fact, Q = rsvd_core(AC, M; svd_kwargs...) - if is_converged!(ndict, fact, new_fact; n_inc, has_qns=hasqns(AC), svd_kwargs...) - break - else - fact = new_fact - end - end - vals = diag(array(new_fact.S)) - (length(vals) == 1 && vals[1]^2 ≤ get(svd_kwargs, :cutoff, 0.0)) && - return nothing, nothing, nothing - #@show flux(dag(CL)*Q*new_fact.U) - #@show flux(new_fact.S)\ - - @assert flux(new_fact.S) == flux(AC) - return dag(CL) * Q * new_fact.U, new_fact.S, new_fact.V * dag(CR) - #ToDo?: handle non-QN case separately because there it is advisable to start with n_init closer to target maxdim_expand - ##not really an issue anymore since we do *2 increase, so only log number of calls - end - - function rsvd_iterative(A::Vector{ITensor}, linds::Vector{<:Index}; svd_kwargs...) - if length(A)==1 - rinds = uniqueinds(only(A), linds) - else - rinds = uniqueinds(A[end],unioninds(A[1:end-1]...)) - end - CL = combiner(linds) - CR = combiner(rinds) - AC=copy(A) - AC[1] = CL*first(AC) - AC[end] = last(AC)*CR - cL = combinedind(CL) - cR = combinedind(CR) - theflux = any(hasqns.(AC)) ? reduce(+,flux.(AC)) : nothing - #@show theflux - n_init = 1 - p_rule(n) = 2 * n - iszero(norm(AC)) && return nothing, nothing, nothing - #@show flux(AC) - ndict2 = init_guess_sizes(cR, n_init, p_rule; theflux) - #@assert isempty(setdiff(keys(Dict(space(cR))), keys(ndict2))) - ndict = init_guess_sizes(cL, n_init, p_rule; theflux) - #FIXME: merging is not the right thing to do, but is a workaround due to the way is_converged! is implemented - ndict = merge(ndict, ndict2) - #@show keys(ndict), keys(Dict(space(cR))) - #@assert isempty(setdiff(keys(Dict(space(cR))), keys(ndict))) - #@assert isempty(setdiff(keys(ndict),keys(Dict(space(cR))))) - - M = build_guess_matrix(eltype(first(AC)), cR, ndict; theflux) - fact, Q = rsvd_core(AC, M; svd_kwargs...) - n_inc = 1 - ndict = increment_guess_sizes(ndict, n_inc, p_rule) - new_fact = deepcopy(fact) - while true - M = build_guess_matrix(eltype(first(AC)), cR, ndict; theflux) - new_fact, Q = rsvd_core(AC, M; svd_kwargs...) - isnothing(Q) && return nothing,nothing,nothing - if is_converged!(ndict, fact, new_fact; n_inc, has_qns=any(hasqns.(AC)), svd_kwargs...) - break - else - fact = new_fact - end - end - vals = diag(array(new_fact.S)) - (length(vals) == 1 && vals[1]^2 ≤ get(svd_kwargs, :cutoff, 0.0)) && - return nothing, nothing, nothing - #@show flux(dag(CL)*Q*new_fact.U) - #@show flux(new_fact.S), theflux - #@assert flux(new_fact.S) == theflux - #@show inds(new_fact.U) - #@show inds(Q) - #@show inds(dag(CL)) - #@show inds(new_fact.S) - #@show inds(new_fact.V) - #@show inds(dag(CR)) - - return dag(CL) * Q * new_fact.U, new_fact.S, new_fact.V * dag(CR) - #ToDo?: handle non-QN case separately because there it is advisable to start with n_init closer to target maxdim_expand - ##not really an issue anymore since we do *2 increase, so only log number of calls - end - - - - - function rsvd(A::ITensor, linds::Vector{<:Index}, n::Int, p::Int; svd_kwargs...) - rinds = uniqueinds(A, linds) - #ToDo handle empty rinds - #boilerplate matricization of tensor for matrix decomp - CL = combiner(linds) - CR = combiner(rinds) - AC = CL * A * CR - cL = combinedind(CL) - cR = combinedind(CR) - if inds(AC) != (cL, cR) - AC = permute(AC, cL, cR) - end - M = build_guess_matrix(eltype(AC), cR, n, p; theflux=hasqns(AC) ? flux(AC) : nothing) - fact, Q = rsvd_core(AC, M; svd_kwargs...) - return dag(CL) * Q * fact.U, fact.S, fact.V * dag(CR) - end - - function rsvd(A::Vector{<:ITensor}, linds::Vector{<:Index}, n::Int, p::Int; svd_kwargs...) - if length(A)==1 - rinds = uniqueinds(only(A), linds) - else - rinds = uniqueinds(A[end],A[end-1]) - end - #ToDo handle empty rinds - #boilerplate matricization of tensor for matrix decomp - CL = combiner(linds) - CR = combiner(rinds) - @assert !isnothing(commonind(CL,first(A))) - @assert !isnothing(commonind(CR,last(A))) - AC=copy(A) - AC[1] = CL*first(AC) - AC[end] = last(AC)*CR - - cL = combinedind(CL) - cR = combinedind(CR) - theflux = any(hasqns.(AC)) ? reduce(+,flux.(AC)) : nothing - #theflux = mapreduce(flux,+,AC) - M = build_guess_matrix(eltype(first(AC)), cR, n, p; theflux) - fact, Q = rsvd_core(AC, M; svd_kwargs...) - return dag(CL) * Q * fact.U, fact.S, fact.V * dag(CR) - end - - function rsvd_core(AC::ITensor, M; svd_kwargs...) - Q = AC * M - #@show dims(Q) - Q = ITensors.qr(Q, commoninds(AC, Q))[1] - QAC = dag(Q) * AC - fact = svd(QAC, uniqueind(QAC, AC); svd_kwargs...) - return fact, Q - end - - function rsvd_core(AC::ITensors.ITensorNetworkMaps.ITensorNetworkMap, M; svd_kwargs...) - #assumes that we want to do a contraction of M with map over its maps output_inds, i.e. a right-multiply - #thus a transpose is necessary - Q = transpose(AC) * M - Q = ITensors.qr(Q, ITensors.ITensorNetworkMaps.input_inds(AC))[1] - QAC = AC * dag(Q) - @assert typeof(QAC) <: ITensor - #@show inds(QAC) - #@assert !iszero(norm(QAC)) - - fact = svd( - QAC, uniqueind(inds(Q), ITensors.ITensorNetworkMaps.input_inds(AC)); svd_kwargs... - ) - return fact, Q - end - -function rsvd_core(AC::Vector{ITensor}, M; svd_kwargs...) - #assumes that we want to do a contraction of M with map over its maps output_inds, i.e. a right-multiply - #thus a transpose is necessary - @assert !isnothing(commonind(last(AC),M)) - Q = foldr(*,AC;init=M) - Q = ITensors.qr(Q, commoninds(Q,first(AC)))[1] - #@show flux(Q) - #@show nnzblocks(Q) - any(isequal(0),dims(Q)) && return nothing, nothing ,nothing ,nothing - QAC = foldl(*,AC,init=dag(Q)) - #@show inds(QAC) - #@show inds(Q) - #@show inds(first(AC)) - #@show inds(last(AC)) - - @assert typeof(QAC) <: ITensor - - fact = svd( - QAC, commoninds(dag(Q), QAC); svd_kwargs... - ) - return fact, Q - end \ No newline at end of file