From ac761235decb43665ed0247bbf132e176fc9a4c3 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 11 Jan 2024 17:16:50 -0500 Subject: [PATCH 01/32] Change required solver function signature. Passes psi and projected Hamiltonian to solver by reference, passes more information about sweep to solver. --- src/treetensornetworks/solvers/tdvp.jl | 25 +++++++++++++++---- src/treetensornetworks/solvers/update_step.jl | 6 ++++- .../test_solvers/test_tdvp.jl | 6 +++-- 3 files changed, 29 insertions(+), 8 deletions(-) diff --git a/src/treetensornetworks/solvers/tdvp.jl b/src/treetensornetworks/solvers/tdvp.jl index 01fe2985..dfc553c6 100644 --- a/src/treetensornetworks/solvers/tdvp.jl +++ b/src/treetensornetworks/solvers/tdvp.jl @@ -1,10 +1,13 @@ function exponentiate_solver() function solver( - H, init; + psi_ref!, + PH_ref!, ishermitian=true, issymmetric=true, region, + sweep_regions, + sweep_step, solver_krylovdim=30, solver_maxiter=100, solver_outputlevel=0, @@ -13,7 +16,15 @@ function exponentiate_solver() normalize, time_step, ) - psi, exp_info = KrylovKit.exponentiate( + #H=copy(PH_ref![]) + H=PH_ref![] ###since we are not changing H we don't need the copy + # let's test whether given region and sweep regions we can find out what the previous and next region were + # this will be needed in subspace expansion + region_ind=sweep_step + next_region=region_ind==length(sweep_regions) ? nothing : first(sweep_regions[region_ind+1]) + previous_region=region_ind==1 ? nothing : first(sweep_regions[region_ind-1]) + + phi, exp_info = KrylovKit.exponentiate( H, time_step, init; @@ -25,16 +36,19 @@ function exponentiate_solver() verbosity=solver_outputlevel, eager=true, ) - return psi, (; info=exp_info) + return phi, (; info=exp_info) end return solver end function applyexp_solver() function solver( - H, init; - tdvp_order, + psi_ref!, + PH_ref!, + region, + sweep_regions, + sweep_step, solver_krylovdim=30, solver_outputlevel=0, solver_tol=1E-8, @@ -42,6 +56,7 @@ function applyexp_solver() time_step, normalize, ) + H=PH_ref![] #applyexp tol is absolute, compute from tol_per_unit_time: tol = abs(time_step) * tol_per_unit_time psi, exp_info = applyexp( diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index 0f69ed52..30e3f29e 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -182,7 +182,11 @@ function local_update( nsites = (region isa AbstractEdge) ? 0 : length(region) PH = set_nsite(PH, nsites) PH = position(PH, psi, region) - phi, info = solver(PH, phi; normalize, region, step_kwargs..., solver_kwargs...) + (psi_ref!) = Ref(psi) # create references, in case solver does (out-of-place) modify PH or psi + (PH_ref!) = Ref(PH) + phi, info = solver(phi;(psi_ref!),(PH_ref!), normalize, region, sweep_regions, sweep_step, step_kwargs..., solver_kwargs...) # args passed by reference are supposed to be modified out of place + psi = psi_ref![] # dereference + PH = PH_ref![] if !(phi isa ITensor && info isa NamedTuple) println("Solver returned the following types: $(typeof(phi)), $(typeof(info))") error("In alternating_update, solver must return an ITensor and a NamedTuple") diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index ec002af1..efa76f02 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -163,10 +163,11 @@ using Test ψ0 = random_mps(s; internal_inds_space=10) - function solver(PH, psi0; time_step, kwargs...) + function solver(psi0; (psi_ref!),(PH_ref!), time_step, kwargs...) solver_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true ) + PH=PH_ref![] psi, exp_info = exponentiate(PH, time_step, psi0; solver_kwargs...) return psi, (; info=exp_info) end @@ -551,10 +552,11 @@ end ψ0 = normalize!(random_ttn(s; link_space=10)) - function solver(PH, psi0; time_step, kwargs...) + function solver(psi0; (psi_ref!),(PH_ref!), time_step, kwargs...) solver_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true ) + PH=PH_ref![] psi, exp_info = exponentiate(PH, time_step, psi0; solver_kwargs...) return psi, (; info=exp_info) end From b23ab06c34d8deb13dc9426c12739760a4bb4d22 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Mon, 15 Jan 2024 20:15:30 -0500 Subject: [PATCH 02/32] Move solver_funcs to separate directory. Adapt dmrg and dmrg-x to new solver interface. --- src/ITensorNetworks.jl | 6 +- src/solvers/applyexp.jl | 27 +++++++ src/solvers/dmrg_x_solver.jl | 19 +++++ src/solvers/eigsolve.jl | 40 ++++++++++ src/solvers/exponentiate.jl | 42 ++++++++++ src/treetensornetworks/solvers/dmrg.jl | 29 +------ src/treetensornetworks/solvers/dmrg_x.jl | 14 +--- src/treetensornetworks/solvers/tdvp.jl | 80 +------------------ src/treetensornetworks/solvers/update_step.jl | 3 +- .../test_solvers/test_tdvp.jl | 13 +-- 10 files changed, 145 insertions(+), 128 deletions(-) create mode 100644 src/solvers/applyexp.jl create mode 100644 src/solvers/dmrg_x_solver.jl create mode 100644 src/solvers/eigsolve.jl create mode 100644 src/solvers/exponentiate.jl diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 5b660bc9..81a2a3d7 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -106,6 +106,11 @@ include("tensornetworkoperators.jl") include(joinpath("ITensorsExt", "itensorutils.jl")) include(joinpath("Graphs", "abstractgraph.jl")) include(joinpath("Graphs", "abstractdatagraph.jl")) +include(joinpath("solvers","eigsolve.jl")) +include(joinpath("solvers","exponentiate.jl")) +include(joinpath("treetensornetworks", "solvers", "applyexp.jl")) #this defines the primitive before the solver function +include(joinpath("solvers","applyexp.jl")) +include(joinpath("solvers","dmrg_x_solver.jl")) include(joinpath("treetensornetworks", "abstracttreetensornetwork.jl")) include(joinpath("treetensornetworks", "ttn.jl")) include(joinpath("treetensornetworks", "opsum_to_ttn.jl")) @@ -114,7 +119,6 @@ include(joinpath("treetensornetworks", "projttns", "projttn.jl")) include(joinpath("treetensornetworks", "projttns", "projttnsum.jl")) include(joinpath("treetensornetworks", "projttns", "projttn_apply.jl")) include(joinpath("treetensornetworks", "solvers", "solver_utils.jl")) -include(joinpath("treetensornetworks", "solvers", "applyexp.jl")) include(joinpath("treetensornetworks", "solvers", "update_step.jl")) include(joinpath("treetensornetworks", "solvers", "alternating_update.jl")) include(joinpath("treetensornetworks", "solvers", "tdvp.jl")) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl new file mode 100644 index 00000000..55a7052d --- /dev/null +++ b/src/solvers/applyexp.jl @@ -0,0 +1,27 @@ +function applyexp_solver() + function solver( + init; + psi_ref!, + PH_ref!, + region, + sweep_regions, + sweep_step, + solver_krylovdim=30, + solver_outputlevel=0, + solver_tol=1E-8, + substep, + time_step, + normalize, + ) + H=PH_ref![] + #applyexp tol is absolute, compute from tol_per_unit_time: + tol = abs(time_step) * tol_per_unit_time + psi, exp_info = applyexp( + H, time_step, init; tol, maxiter=solver_krylovdim, outputlevel=solver_outputlevel + ) + return psi, (; info=exp_info) + end + return solver +end + + diff --git a/src/solvers/dmrg_x_solver.jl b/src/solvers/dmrg_x_solver.jl new file mode 100644 index 00000000..33b1ad2a --- /dev/null +++ b/src/solvers/dmrg_x_solver.jl @@ -0,0 +1,19 @@ +function dmrg_x_solver( + init; + psi_ref!, + PH_ref!, + normalize=nothing, + region, + sweep_regions, + sweep_step, + half_sweep, + step_kwargs... + ) + H = contract(PH_ref![], ITensor(1.0)) + D, U = eigen(H; ishermitian=true) + u = uniqueind(U, H) + max_overlap, max_ind = findmax(abs, array(dag(init) * U)) + U_max = U * dag(onehot(u => max_ind)) + # TODO: improve this to return the energy estimate too + return U_max, NamedTuple() +end \ No newline at end of file diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl new file mode 100644 index 00000000..da3ba5c1 --- /dev/null +++ b/src/solvers/eigsolve.jl @@ -0,0 +1,40 @@ + +function eigsolve_solver(; + solver_which_eigenvalue=:SR, #TODO: settle on pattern to pass solver kwargs + ishermitian=true, + solver_tol=1e-14, + solver_krylovdim=3, + solver_maxiter=1, + solver_verbosity=0, + ) + + function solver( + init; + psi_ref!, + PH_ref!, + normalize, + region, + sweep_regions, + sweep_step, + sweep_kwargs... + # slurp solver_kwargs? #TODO: homogenize how the solver kwargs are passed + ) + howmany = 1 + which = solver_which_eigenvalue + vals, vecs, info = eigsolve( + PH_ref![], + init, + howmany, + which; + ishermitian, + tol=solver_tol, + krylovdim=solver_krylovdim, + maxiter=solver_maxiter, + verbosity=solver_verbosity, + ) + phi = vecs[1] + return phi, (; solver_info=info, energies=vals) + end + return solver + end + diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl new file mode 100644 index 00000000..bb2e98e2 --- /dev/null +++ b/src/solvers/exponentiate.jl @@ -0,0 +1,42 @@ +function exponentiate_solver() + function solver( + init; + psi_ref!, + PH_ref!, + ishermitian=true, + issymmetric=true, + region, + sweep_regions, + sweep_step, + solver_krylovdim=30, + solver_maxiter=100, + solver_outputlevel=0, + solver_tol=1E-12, + substep, + normalize, + time_step, + ) + #H=copy(PH_ref![]) + H=PH_ref![] ###since we are not changing H we don't need the copy + # let's test whether given region and sweep regions we can find out what the previous and next region were + # this will be needed in subspace expansion + region_ind=sweep_step + next_region=region_ind==length(sweep_regions) ? nothing : first(sweep_regions[region_ind+1]) + previous_region=region_ind==1 ? nothing : first(sweep_regions[region_ind-1]) + + phi, exp_info = KrylovKit.exponentiate( + H, + time_step, + init; + ishermitian, + issymmetric, + tol=solver_tol, + krylovdim=solver_krylovdim, + maxiter=solver_maxiter, + verbosity=solver_outputlevel, + eager=true, + ) + return phi, (; info=exp_info) + end + return solver +end diff --git a/src/treetensornetworks/solvers/dmrg.jl b/src/treetensornetworks/solvers/dmrg.jl index c67b3b86..8f24fd63 100644 --- a/src/treetensornetworks/solvers/dmrg.jl +++ b/src/treetensornetworks/solvers/dmrg.jl @@ -1,31 +1,3 @@ -function eigsolve_solver(; - solver_which_eigenvalue=:SR, - ishermitian=true, - solver_tol=1e-14, - solver_krylovdim=3, - solver_maxiter=1, - solver_verbosity=0, -) - function solver(H, init; normalize=nothing, region=nothing, half_sweep=nothing) - howmany = 1 - which = solver_which_eigenvalue - vals, vecs, info = eigsolve( - H, - init, - howmany, - which; - ishermitian, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_verbosity, - ) - psi = vecs[1] - return psi, (; solver_info=info, energies=vals) - end - return solver -end - """ Overload of `ITensors.dmrg`. """ @@ -51,6 +23,7 @@ function dmrg( ), H, init; + reverse_step=false, kwargs..., ) end diff --git a/src/treetensornetworks/solvers/dmrg_x.jl b/src/treetensornetworks/solvers/dmrg_x.jl index 30a97f3a..67ab42c9 100644 --- a/src/treetensornetworks/solvers/dmrg_x.jl +++ b/src/treetensornetworks/solvers/dmrg_x.jl @@ -1,16 +1,4 @@ -function dmrg_x_solver( - PH, init; normalize=nothing, region=nothing, half_sweep=nothing, reverse_step=nothing -) - H = contract(PH, ITensor(1.0)) - D, U = eigen(H; ishermitian=true) - u = uniqueind(U, H) - max_overlap, max_ind = findmax(abs, array(dag(init) * U)) - U_max = U * dag(onehot(u => max_ind)) - # TODO: improve this to return the energy estimate too - return U_max, NamedTuple() -end - function dmrg_x(PH, init::AbstractTTN; kwargs...) - psi = alternating_update(dmrg_x_solver, PH, init; kwargs...) + psi = alternating_update(ITensorNetworks.dmrg_x_solver, PH, init; kwargs...) return psi end diff --git a/src/treetensornetworks/solvers/tdvp.jl b/src/treetensornetworks/solvers/tdvp.jl index dfc553c6..863a36df 100644 --- a/src/treetensornetworks/solvers/tdvp.jl +++ b/src/treetensornetworks/solvers/tdvp.jl @@ -1,72 +1,3 @@ -function exponentiate_solver() - function solver( - init; - psi_ref!, - PH_ref!, - ishermitian=true, - issymmetric=true, - region, - sweep_regions, - sweep_step, - solver_krylovdim=30, - solver_maxiter=100, - solver_outputlevel=0, - solver_tol=1E-12, - substep, - normalize, - time_step, - ) - #H=copy(PH_ref![]) - H=PH_ref![] ###since we are not changing H we don't need the copy - # let's test whether given region and sweep regions we can find out what the previous and next region were - # this will be needed in subspace expansion - region_ind=sweep_step - next_region=region_ind==length(sweep_regions) ? nothing : first(sweep_regions[region_ind+1]) - previous_region=region_ind==1 ? nothing : first(sweep_regions[region_ind-1]) - - phi, exp_info = KrylovKit.exponentiate( - H, - time_step, - init; - ishermitian, - issymmetric, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_outputlevel, - eager=true, - ) - return phi, (; info=exp_info) - end - return solver -end - -function applyexp_solver() - function solver( - init; - psi_ref!, - PH_ref!, - region, - sweep_regions, - sweep_step, - solver_krylovdim=30, - solver_outputlevel=0, - solver_tol=1E-8, - substep, - time_step, - normalize, - ) - H=PH_ref![] - #applyexp tol is absolute, compute from tol_per_unit_time: - tol = abs(time_step) * tol_per_unit_time - psi, exp_info = applyexp( - H, time_step, init; tol, maxiter=solver_krylovdim, outputlevel=solver_outputlevel - ) - return psi, (; info=exp_info) - end - return solver -end - function _compute_nsweeps(nsteps, t, time_step, order) nsweeps_per_step = order / 2 nsweeps = 1 @@ -183,15 +114,6 @@ Optional keyword arguments: * `observer` - object implementing the Observer interface which can perform measurements and stop early * `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations """ -function tdvp(H, t::Number, init::AbstractTTN; solver_backend="exponentiate", kwargs...) - if solver_backend == "exponentiate" - solver = exponentiate_solver - elseif solver_backend == "applyexp" - solver = applyexp_solver - else - error( - "solver_backend=$solver_backend not recognized (options are \"applyexp\" or \"exponentiate\")", - ) - end +function tdvp(H, t::Number, init::AbstractTTN; solver=exponentiate_solver, kwargs...) return tdvp(solver(), H, t, init; kwargs...) end diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index 30e3f29e..a62f7b70 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -44,7 +44,8 @@ function update_step( step_printer=step_printer, (step_observer!)=observer(), sweep::Int=1, - sweep_regions=default_sweep_regions(nsite, psi), + reverse_step::Bool=false, + sweep_regions=default_sweep_regions(nsite, psi;reverse_step), kwargs..., ) PH = copy(PH) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index efa76f02..1d50fba1 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -4,6 +4,7 @@ using KrylovKit: exponentiate using Observers using Random using Test +#exponentiate_solver=ITensorNetworks.exponentiate_solver #ToDo: how to best handle importing etc. @testset "MPS TDVP" begin @testset "Basic TDVP" begin @@ -31,7 +32,7 @@ using Test # #Different backend solvers, default solver_backend = "applyexp" ψ1_exponentiate_backend = tdvp( - H, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, solver_backend="exponentiate" + H, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, solver=exponentiate_solver ) @test ψ1 ≈ ψ1_exponentiate_backend rtol = 1e-7 @@ -48,7 +49,7 @@ using Test #Different backend solvers, default solver_backend = "applyexp" ψ2_exponentiate_backend = tdvp( - H, +0.1im, ψ1; nsteps=1, cutoff, solver_backend="exponentiate" + H, +0.1im, ψ1; nsteps=1, cutoff, solver=exponentiate_solver ) @test ψ2 ≈ ψ2_exponentiate_backend rtol = 1e-7 @@ -84,7 +85,7 @@ using Test #Different backend solvers, default solver_backend = "applyexp" ψ1_exponentiate_backend = tdvp( - Hs, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, solver_backend="exponentiate" + Hs, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, solver=exponentiate_solver ) @test ψ1 ≈ ψ1_exponentiate_backend rtol = 1e-7 @@ -101,7 +102,7 @@ using Test #Different backend solvers, default solver_backend = "applyexp" ψ2_exponentiate_backend = tdvp( - Hs, +0.1im, ψ1; nsteps=1, cutoff, solver_backend="exponentiate" + Hs, +0.1im, ψ1; nsteps=1, cutoff, solver=exponentiate_solver ) @test ψ2 ≈ ψ2_exponentiate_backend rtol = 1e-7 @@ -249,7 +250,7 @@ using Test solver_tol=1e-12, solver_maxiter=500, solver_krylovdim=25, - solver_backend="exponentiate", + solver=exponentiate_solver, ) # TODO: What should `expect` output? Right now # it outputs a dictionary. @@ -385,7 +386,7 @@ using Test reverse_step, normalize=true, solver_krylovdim=15, - solver_backend="exponentiate", + solver=ITensorNetworks.exponentiate_solver, ) end From c2c545552566dfd5d935c682ad431ac552a4bc2f Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Tue, 16 Jan 2024 10:45:44 -0500 Subject: [PATCH 03/32] Format. --- src/ITensorNetworks.jl | 8 +- src/solvers/applyexp.jl | 48 ++++++------ src/solvers/dmrg_x_solver.jl | 24 +++--- src/solvers/eigsolve.jl | 70 ++++++++--------- src/solvers/exponentiate.jl | 77 ++++++++++--------- src/treetensornetworks/solvers/update_step.jl | 14 +++- .../test_solvers/test_tdvp.jl | 8 +- 7 files changed, 128 insertions(+), 121 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 81a2a3d7..8ef00052 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -106,11 +106,11 @@ include("tensornetworkoperators.jl") include(joinpath("ITensorsExt", "itensorutils.jl")) include(joinpath("Graphs", "abstractgraph.jl")) include(joinpath("Graphs", "abstractdatagraph.jl")) -include(joinpath("solvers","eigsolve.jl")) -include(joinpath("solvers","exponentiate.jl")) +include(joinpath("solvers", "eigsolve.jl")) +include(joinpath("solvers", "exponentiate.jl")) include(joinpath("treetensornetworks", "solvers", "applyexp.jl")) #this defines the primitive before the solver function -include(joinpath("solvers","applyexp.jl")) -include(joinpath("solvers","dmrg_x_solver.jl")) +include(joinpath("solvers", "applyexp.jl")) +include(joinpath("solvers", "dmrg_x_solver.jl")) include(joinpath("treetensornetworks", "abstracttreetensornetwork.jl")) include(joinpath("treetensornetworks", "ttn.jl")) include(joinpath("treetensornetworks", "opsum_to_ttn.jl")) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl index 55a7052d..222bb9b6 100644 --- a/src/solvers/applyexp.jl +++ b/src/solvers/applyexp.jl @@ -1,27 +1,25 @@ function applyexp_solver() - function solver( - init; - psi_ref!, - PH_ref!, - region, - sweep_regions, - sweep_step, - solver_krylovdim=30, - solver_outputlevel=0, - solver_tol=1E-8, - substep, - time_step, - normalize, - ) - H=PH_ref![] - #applyexp tol is absolute, compute from tol_per_unit_time: - tol = abs(time_step) * tol_per_unit_time - psi, exp_info = applyexp( - H, time_step, init; tol, maxiter=solver_krylovdim, outputlevel=solver_outputlevel - ) - return psi, (; info=exp_info) - end - return solver + function solver( + init; + psi_ref!, + PH_ref!, + region, + sweep_regions, + sweep_step, + solver_krylovdim=30, + solver_outputlevel=0, + solver_tol=1E-8, + substep, + time_step, + normalize, + ) + H = PH_ref![] + #applyexp tol is absolute, compute from tol_per_unit_time: + tol = abs(time_step) * tol_per_unit_time + psi, exp_info = applyexp( + H, time_step, init; tol, maxiter=solver_krylovdim, outputlevel=solver_outputlevel + ) + return psi, (; info=exp_info) + end + return solver end - - diff --git a/src/solvers/dmrg_x_solver.jl b/src/solvers/dmrg_x_solver.jl index 33b1ad2a..6db73e9b 100644 --- a/src/solvers/dmrg_x_solver.jl +++ b/src/solvers/dmrg_x_solver.jl @@ -1,19 +1,19 @@ function dmrg_x_solver( - init; - psi_ref!, - PH_ref!, - normalize=nothing, - region, - sweep_regions, - sweep_step, - half_sweep, - step_kwargs... - ) + init; + psi_ref!, + PH_ref!, + normalize=nothing, + region, + sweep_regions, + sweep_step, + half_sweep, + step_kwargs..., +) H = contract(PH_ref![], ITensor(1.0)) - D, U = eigen(H; ishermitian=true) + D, U = eigen(H; ishermitian=true) u = uniqueind(U, H) max_overlap, max_ind = findmax(abs, array(dag(init) * U)) U_max = U * dag(onehot(u => max_ind)) # TODO: improve this to return the energy estimate too return U_max, NamedTuple() -end \ No newline at end of file +end diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index da3ba5c1..a1dfa39c 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -1,40 +1,38 @@ function eigsolve_solver(; - solver_which_eigenvalue=:SR, #TODO: settle on pattern to pass solver kwargs - ishermitian=true, - solver_tol=1e-14, - solver_krylovdim=3, - solver_maxiter=1, - solver_verbosity=0, + solver_which_eigenvalue=:SR, #TODO: settle on pattern to pass solver kwargs + ishermitian=true, + solver_tol=1e-14, + solver_krylovdim=3, + solver_maxiter=1, + solver_verbosity=0, +) + function solver( + init; + psi_ref!, + PH_ref!, + normalize, + region, + sweep_regions, + sweep_step, + sweep_kwargs..., + # slurp solver_kwargs? #TODO: homogenize how the solver kwargs are passed ) - - function solver( - init; - psi_ref!, - PH_ref!, - normalize, - region, - sweep_regions, - sweep_step, - sweep_kwargs... - # slurp solver_kwargs? #TODO: homogenize how the solver kwargs are passed - ) - howmany = 1 - which = solver_which_eigenvalue - vals, vecs, info = eigsolve( - PH_ref![], - init, - howmany, - which; - ishermitian, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_verbosity, - ) - phi = vecs[1] - return phi, (; solver_info=info, energies=vals) - end - return solver + howmany = 1 + which = solver_which_eigenvalue + vals, vecs, info = eigsolve( + PH_ref![], + init, + howmany, + which; + ishermitian, + tol=solver_tol, + krylovdim=solver_krylovdim, + maxiter=solver_maxiter, + verbosity=solver_verbosity, + ) + phi = vecs[1] + return phi, (; solver_info=info, energies=vals) end - + return solver +end diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index bb2e98e2..b7d1f95a 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -1,42 +1,43 @@ function exponentiate_solver() - function solver( - init; - psi_ref!, - PH_ref!, - ishermitian=true, - issymmetric=true, - region, - sweep_regions, - sweep_step, - solver_krylovdim=30, - solver_maxiter=100, - solver_outputlevel=0, - solver_tol=1E-12, - substep, - normalize, + function solver( + init; + psi_ref!, + PH_ref!, + ishermitian=true, + issymmetric=true, + region, + sweep_regions, + sweep_step, + solver_krylovdim=30, + solver_maxiter=100, + solver_outputlevel=0, + solver_tol=1E-12, + substep, + normalize, + time_step, + ) + #H=copy(PH_ref![]) + H = PH_ref![] ###since we are not changing H we don't need the copy + # let's test whether given region and sweep regions we can find out what the previous and next region were + # this will be needed in subspace expansion + region_ind = sweep_step + next_region = + region_ind == length(sweep_regions) ? nothing : first(sweep_regions[region_ind + 1]) + previous_region = region_ind == 1 ? nothing : first(sweep_regions[region_ind - 1]) + + phi, exp_info = KrylovKit.exponentiate( + H, time_step, + init; + ishermitian, + issymmetric, + tol=solver_tol, + krylovdim=solver_krylovdim, + maxiter=solver_maxiter, + verbosity=solver_outputlevel, + eager=true, ) - #H=copy(PH_ref![]) - H=PH_ref![] ###since we are not changing H we don't need the copy - # let's test whether given region and sweep regions we can find out what the previous and next region were - # this will be needed in subspace expansion - region_ind=sweep_step - next_region=region_ind==length(sweep_regions) ? nothing : first(sweep_regions[region_ind+1]) - previous_region=region_ind==1 ? nothing : first(sweep_regions[region_ind-1]) - - phi, exp_info = KrylovKit.exponentiate( - H, - time_step, - init; - ishermitian, - issymmetric, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_outputlevel, - eager=true, - ) - return phi, (; info=exp_info) - end - return solver + return phi, (; info=exp_info) + end + return solver end diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index a62f7b70..d538fbc9 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -45,7 +45,7 @@ function update_step( (step_observer!)=observer(), sweep::Int=1, reverse_step::Bool=false, - sweep_regions=default_sweep_regions(nsite, psi;reverse_step), + sweep_regions=default_sweep_regions(nsite, psi; reverse_step), kwargs..., ) PH = copy(PH) @@ -185,7 +185,17 @@ function local_update( PH = position(PH, psi, region) (psi_ref!) = Ref(psi) # create references, in case solver does (out-of-place) modify PH or psi (PH_ref!) = Ref(PH) - phi, info = solver(phi;(psi_ref!),(PH_ref!), normalize, region, sweep_regions, sweep_step, step_kwargs..., solver_kwargs...) # args passed by reference are supposed to be modified out of place + phi, info = solver( + phi; + (psi_ref!), + (PH_ref!), + normalize, + region, + sweep_regions, + sweep_step, + step_kwargs..., + solver_kwargs..., + ) # args passed by reference are supposed to be modified out of place psi = psi_ref![] # dereference PH = PH_ref![] if !(phi isa ITensor && info isa NamedTuple) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index 1d50fba1..81ebeb8b 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -164,11 +164,11 @@ using Test ψ0 = random_mps(s; internal_inds_space=10) - function solver(psi0; (psi_ref!),(PH_ref!), time_step, kwargs...) + function solver(psi0; (psi_ref!), (PH_ref!), time_step, kwargs...) solver_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true ) - PH=PH_ref![] + PH = PH_ref![] psi, exp_info = exponentiate(PH, time_step, psi0; solver_kwargs...) return psi, (; info=exp_info) end @@ -553,11 +553,11 @@ end ψ0 = normalize!(random_ttn(s; link_space=10)) - function solver(psi0; (psi_ref!),(PH_ref!), time_step, kwargs...) + function solver(psi0; (psi_ref!), (PH_ref!), time_step, kwargs...) solver_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true ) - PH=PH_ref![] + PH = PH_ref![] psi, exp_info = exponentiate(PH, time_step, psi0; solver_kwargs...) return psi, (; info=exp_info) end From ed0df91728607d8cae52ec898439b634b0287c56 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 17 Jan 2024 17:04:07 -0500 Subject: [PATCH 04/32] Modify alternating update kwarg naming, structure, and also solver_interface. --- src/solvers/exponentiate.jl | 56 +++---- .../solvers/alternating_update.jl | 47 +++--- src/treetensornetworks/solvers/tdvp.jl | 17 ++- src/treetensornetworks/solvers/update_step.jl | 138 +++++++++--------- .../test_solvers/test_tdvp.jl | 101 +++++++------ 5 files changed, 188 insertions(+), 171 deletions(-) diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index b7d1f95a..e70edc80 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -1,43 +1,47 @@ -function exponentiate_solver() - function solver( + function exponentiate_updater( init; psi_ref!, PH_ref!, + outputlevel, + which_sweep, # keep, change name + region_updates, + which_region_update, + region_kwargs, # region_kwargs (timestep for solver) + updater_kwargs, + ) + default_updater_kwargs=(; + krylovdim=30, #from here only solver kwargs + maxiter=100, + outputlevel=0, + tol=1E-12, ishermitian=true, issymmetric=true, - region, - sweep_regions, - sweep_step, - solver_krylovdim=30, - solver_maxiter=100, - solver_outputlevel=0, - solver_tol=1E-12, - substep, - normalize, - time_step, - ) + eager=true + ) + updater_kwargs=merge(default_updater_kwargs,updater_kwargs) #last collection has precedence #H=copy(PH_ref![]) H = PH_ref![] ###since we are not changing H we don't need the copy # let's test whether given region and sweep regions we can find out what the previous and next region were # this will be needed in subspace expansion - region_ind = sweep_step - next_region = - region_ind == length(sweep_regions) ? nothing : first(sweep_regions[region_ind + 1]) - previous_region = region_ind == 1 ? nothing : first(sweep_regions[region_ind - 1]) + #@show step_kwargs + substep=get(region_kwargs,:substep,nothing) + time_step=get(region_kwargs,:time_step,nothing) + @assert !isnothing(time_step) && !isnothing(substep) + region_ind = which_region_update + next_region = region_ind == length(region_updates) ? nothing : first(region_updates[region_ind + 1]) + previous_region = region_ind == 1 ? nothing : first(region_updates[region_ind - 1]) phi, exp_info = KrylovKit.exponentiate( H, time_step, init; - ishermitian, - issymmetric, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_outputlevel, - eager=true, + ishermitian=updater_kwargs[:ishermitian], + issymmetric=updater_kwargs[:issymmetric], + tol=updater_kwargs[:tol], + krylovdim=updater_kwargs[:krylovdim], + maxiter=updater_kwargs[:maxiter], + verbosity=updater_kwargs[:outputlevel], + eager=updater_kwargs[:eager], ) return phi, (; info=exp_info) end - return solver -end diff --git a/src/treetensornetworks/solvers/alternating_update.jl b/src/treetensornetworks/solvers/alternating_update.jl index 17d1f62d..04b5ce37 100644 --- a/src/treetensornetworks/solvers/alternating_update.jl +++ b/src/treetensornetworks/solvers/alternating_update.jl @@ -26,9 +26,9 @@ function process_sweeps( return maxdim, mindim, cutoff, noise, kwargs end -function sweep_printer(; outputlevel, psi, sweep, sw_time) +function sweep_printer(; outputlevel, psi, which_sweep, sw_time) if outputlevel >= 1 - print("After sweep ", sweep, ":") + print("After sweep ", which_sweep, ":") print(" maxlinkdim=", maxlinkdim(psi)) print(" cpu_time=", round(sw_time; digits=3)) println() @@ -37,7 +37,7 @@ function sweep_printer(; outputlevel, psi, sweep, sw_time) end function alternating_update( - solver, + updater, PH, psi0::AbstractTTN; checkdone=(; kws...) -> false, @@ -46,55 +46,58 @@ function alternating_update( (sweep_observer!)=observer(), sweep_printer=sweep_printer, write_when_maxdim_exceeds::Union{Int,Nothing}=nothing, + updater_kwargs, kwargs..., ) maxdim, mindim, cutoff, noise, kwargs = process_sweeps(nsweeps; kwargs...) psi = copy(psi0) - insert_function!(sweep_observer!, "sweep_printer" => sweep_printer) + insert_function!(sweep_observer!, "sweep_printer" => sweep_printer) # FIX THIS - for sweep in 1:nsweeps - if !isnothing(write_when_maxdim_exceeds) && maxdim[sweep] > write_when_maxdim_exceeds + for which_sweep in 1:nsweeps + if !isnothing(write_when_maxdim_exceeds) && maxdim[which_sweep] > write_when_maxdim_exceeds if outputlevel >= 2 println( - "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim[sweep] = $(maxdim[sweep]), writing environment tensors to disk", + "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim[which_sweep] = $(maxdim[which_sweep]), writing environment tensors to disk", ) end PH = disk(PH) end - + sw_time = @elapsed begin - psi, PH = update_step( - solver, + psi, PH = sweep_update( + updater, PH, psi; outputlevel, - sweep, - maxdim=maxdim[sweep], - mindim=mindim[sweep], - cutoff=cutoff[sweep], - noise=noise[sweep], + which_sweep, + sweep_params=(; + maxdim=maxdim[which_sweep], + mindim=mindim[which_sweep], + cutoff=cutoff[which_sweep], + noise=noise[which_sweep]), + updater_kwargs, kwargs..., ) end - update!(sweep_observer!; psi, sweep, sw_time, outputlevel) + update!(sweep_observer!; psi, which_sweep, sw_time, outputlevel) - checkdone(; psi, sweep, outputlevel, kwargs...) && break + checkdone(; psi, which_sweep, outputlevel, kwargs...) && break end - select!(sweep_observer!, Observers.DataFrames.Not("sweep_printer")) # remove sweep_printer + select!(sweep_observer!, Observers.DataFrames.Not("sweep_printer")) return psi end -function alternating_update(solver, H::AbstractTTN, psi0::AbstractTTN; kwargs...) +function alternating_update(updater, H::AbstractTTN, psi0::AbstractTTN; kwargs...) check_hascommoninds(siteinds, H, psi0) check_hascommoninds(siteinds, H, psi0') # Permute the indices to have a better memory layout # and minimize permutations H = ITensors.permute(H, (linkind, siteinds, linkind)) PH = ProjTTN(H) - return alternating_update(solver, PH, psi0; kwargs...) + return alternating_update(updater, PH, psi0; kwargs...) end """ @@ -116,12 +119,12 @@ each step of the algorithm when optimizing the MPS. Returns: * `psi::MPS` - time-evolved MPS """ -function alternating_update(solver, Hs::Vector{<:AbstractTTN}, psi0::AbstractTTN; kwargs...) +function alternating_update(updater, Hs::Vector{<:AbstractTTN}, psi0::AbstractTTN; kwargs...) for H in Hs check_hascommoninds(siteinds, H, psi0) check_hascommoninds(siteinds, H, psi0') end Hs .= ITensors.permute.(Hs, Ref((linkind, siteinds, linkind))) PHs = ProjTTNSum(Hs) - return alternating_update(solver, PHs, psi0; kwargs...) + return alternating_update(updater, PHs, psi0; kwargs...) end diff --git a/src/treetensornetworks/solvers/tdvp.jl b/src/treetensornetworks/solvers/tdvp.jl index 863a36df..badac5b5 100644 --- a/src/treetensornetworks/solvers/tdvp.jl +++ b/src/treetensornetworks/solvers/tdvp.jl @@ -57,7 +57,7 @@ function tdvp_sweep( end function tdvp( - solver, + updater, H, t::Number, init::AbstractTTN; @@ -68,17 +68,18 @@ function tdvp( (sweep_observer!)=observer(), root_vertex=default_root_vertex(init), reverse_step=true, + updater_kwargs=NamedTuple(;), kwargs..., ) nsweeps = _compute_nsweeps(nsteps, t, time_step, order) - sweep_regions = tdvp_sweep(order, nsite, time_step, init; root_vertex, reverse_step) + region_updates = tdvp_sweep(order, nsite, time_step, init; root_vertex, reverse_step) - function sweep_time_printer(; outputlevel, sweep, kwargs...) + function sweep_time_printer(; outputlevel, which_sweep, kwargs...) if outputlevel >= 1 sweeps_per_step = order ÷ 2 if sweep % sweeps_per_step == 0 - current_time = (sweep / sweeps_per_step) * time_step - println("Current time (sweep $sweep) = ", round(current_time; digits=3)) + current_time = (which_sweep / sweeps_per_step) * time_step + println("Current time (sweep $which_sweep) = ", round(current_time; digits=3)) end end return nothing @@ -87,7 +88,7 @@ function tdvp( insert_function!(sweep_observer!, "sweep_time_printer" => sweep_time_printer) psi = alternating_update( - solver, H, init; nsweeps, sweep_observer!, sweep_regions, nsite, kwargs... + updater, H, init; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... ) # remove sweep_time_printer from sweep_observer! @@ -114,6 +115,6 @@ Optional keyword arguments: * `observer` - object implementing the Observer interface which can perform measurements and stop early * `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations """ -function tdvp(H, t::Number, init::AbstractTTN; solver=exponentiate_solver, kwargs...) - return tdvp(solver(), H, t, init; kwargs...) +function tdvp(H, t::Number, init::AbstractTTN; updater=exponentiate_updater, kwargs...) + return tdvp(updater, H, t, init; kwargs...) end diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index d538fbc9..dfc5ce62 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -1,5 +1,5 @@ -function default_sweep_regions(nsite, graph::AbstractGraph; kwargs...) +function default_sweep_regions(nsite, graph::AbstractGraph; kwargs...) ###move this to a different file, algorithmic level idea return vcat( [ half_sweep( @@ -14,11 +14,12 @@ function default_sweep_regions(nsite, graph::AbstractGraph; kwargs...) ) end -function step_printer(; - cutoff, maxdim, mindim, outputlevel::Int=0, psi, region, spec, sweep_step +function region_printer(; + cutoff, maxdim, mindim, outputlevel::Int=0, psi, region_updates, spec, which_region_update, which_sweep,kwargs... ) if outputlevel >= 2 - @printf("Sweep %d, region=%s \n", sweep, region) + region=first(region_updates[which_region_update]) + @printf("Sweep %d, region=%s \n", which_sweep, region) print(" Truncated using") @printf(" cutoff=%.1E", cutoff) @printf(" maxdim=%d", maxdim) @@ -35,51 +36,51 @@ function step_printer(; end end -function update_step( +function sweep_update( solver, PH, psi::AbstractTTN; - normalize::Bool=false, - nsite::Int=2, + normalize::Bool=false, # ToDo: think about where to put the default, probably this default is best defined at algorithmic level + outputlevel, step_printer=step_printer, - (step_observer!)=observer(), - sweep::Int=1, - reverse_step::Bool=false, - sweep_regions=default_sweep_regions(nsite, psi; reverse_step), - kwargs..., + (region_observer!)=observer(), # ToDo: change name to region_observer! ? + which_sweep::Int, + sweep_params::NamedTuple, + region_updates,# =default_sweep_regions(nsite, psi; reverse_step), #move default up to algorithmic level + updater_kwargs, ) - PH = copy(PH) - psi = copy(psi) - - insert_function!(step_observer!, "step_printer" => step_printer) + insert_function!(region_observer!, "region_printer" => region_printer) #ToDo fix this # Append empty namedtuple to each element if not already present - # (Needed to handle user-provided sweep_regions) - sweep_regions = append_missing_namedtuple.(to_tuple.(sweep_regions)) + # (Needed to handle user-provided region_updates) + region_updates = append_missing_namedtuple.(to_tuple.(region_updates)) if nv(psi) == 1 error( "`alternating_update` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.exponentiate`, etc.", ) end - - for (sweep_step, (region, step_kwargs)) in enumerate(sweep_regions) - psi, PH = local_update( + + for which_region_update in eachindex(region_updates) + # merge sweep params in step_kwargs + (region, region_kwargs)=region_updates[which_region_update] + region_kwargs=merge(region_kwargs, sweep_params) #in this case sweep params has precedence over step_kwargs --- correct behaviour? + psi, PH = region_update( solver, PH, - psi, - region; + psi; normalize, - step_kwargs, - step_observer!, - sweep, - sweep_regions, - sweep_step, - kwargs..., + outputlevel, + which_sweep, + region_updates, + which_region_update, + region_kwargs, + region_observer!, + updater_kwargs, ) end - select!(step_observer!, Observers.DataFrames.Not("step_printer")) # remove step_printer + select!(region_observer!, Observers.DataFrames.Not("region_printer")) # remove step_printer #todo fix this # Just to be sure: normalize && normalize!(psi) @@ -155,46 +156,44 @@ function insert_local_tensor(psi::AbstractTTN, phi::ITensor, e::NamedEdge; kwarg end #TODO: clean this up: +# also can we entirely rely on directionality of edges by construction? current_ortho(::Type{<:Vector{<:V}}, st) where {V} = first(st) current_ortho(::Type{NamedEdge{V}}, st) where {V} = src(st) current_ortho(st) = current_ortho(typeof(st), st) -function local_update( - solver, +function region_update( + updater, PH, - psi, - region; + psi; normalize, - noise, - cutoff::AbstractFloat=1E-16, - maxdim::Int=typemax(Int), - mindim::Int=1, - outputlevel::Int=0, - step_kwargs=NamedTuple(), - step_observer!, - sweep, - sweep_regions, - sweep_step, - solver_kwargs..., + outputlevel, + which_sweep, + region_updates, + which_region_update, + region_kwargs, + region_observer!, + #insertion_kwargs, #ToDo: later + #extraction_kwargs, #ToDo: implement later with possibility to pass custom extraction/insertion func (or code into func) + updater_kwargs ) + region=first(region_updates[which_region_update]) psi = orthogonalize(psi, current_ortho(region)) - psi, phi = extract_local_tensor(psi, region) - - nsites = (region isa AbstractEdge) ? 0 : length(region) + psi, phi = extract_local_tensor(psi, region;) + nsites = (region isa AbstractEdge) ? 0 : length(region) #ToDo move into separate funtion PH = set_nsite(PH, nsites) PH = position(PH, psi, region) - (psi_ref!) = Ref(psi) # create references, in case solver does (out-of-place) modify PH or psi - (PH_ref!) = Ref(PH) - phi, info = solver( + psi_ref! = Ref(psi) # create references, in case solver does (out-of-place) modify PH or psi + PH_ref! = Ref(PH) + phi, info = updater( phi; - (psi_ref!), - (PH_ref!), - normalize, - region, - sweep_regions, - sweep_step, - step_kwargs..., - solver_kwargs..., + psi_ref!, + PH_ref!, + outputlevel, + which_sweep, + region_updates, + which_region_update, + region_kwargs, + updater_kwargs ) # args passed by reference are supposed to be modified out of place psi = psi_ref![] # dereference PH = PH_ref![] @@ -205,30 +204,35 @@ function local_update( normalize && (phi /= norm(phi)) drho = nothing - ortho = "left" + ortho = "left" #i guess with respect to ordered vertices that's valid but may be cleaner to use next_region logic #if noise > 0.0 && isforward(direction) # drho = noise * noiseterm(PH, phi, ortho) # TODO: actually implement this for trees... + # so noiseterm is a solver #end psi, spec = insert_local_tensor( - psi, phi, region; eigen_perturbation=drho, ortho, normalize, maxdim, mindim, cutoff + psi, phi, region; eigen_perturbation=drho, ortho, normalize, + maxdim=region_kwargs[:maxdim], + mindim=region_kwargs[:mindim], + cutoff=region_kwargs[:cutoff] ) update!( - step_observer!; + region_observer!; cutoff, maxdim, mindim, - sweep_step, - total_sweep_steps=length(sweep_regions), - end_of_sweep=(sweep_step == length(sweep_regions)), + which_region_update, + region_updates, + total_sweep_steps=length(region_updates), + end_of_sweep=(which_region_update == length(region_updates)), psi, region, - sweep, + which_sweep, spec, outputlevel, info..., - step_kwargs..., + region_kwargs..., ) return psi, PH end diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index 81ebeb8b..f8fadb1b 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -4,7 +4,7 @@ using KrylovKit: exponentiate using Observers using Random using Test -#exponentiate_solver=ITensorNetworks.exponentiate_solver #ToDo: how to best handle importing etc. +exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best handle importing etc. @testset "MPS TDVP" begin @testset "Basic TDVP" begin @@ -30,9 +30,9 @@ using Test # #TODO: exponentiate is now the default, so switch this to applyexp # - #Different backend solvers, default solver_backend = "applyexp" + #Different backend updaters, default updater_backend = "applyexp" ψ1_exponentiate_backend = tdvp( - H, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, solver=exponentiate_solver + H, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, updater=exponentiate_updater ) @test ψ1 ≈ ψ1_exponentiate_backend rtol = 1e-7 @@ -47,9 +47,9 @@ using Test # Time evolve backwards: ψ2 = tdvp(H, +0.1im, ψ1; nsteps=1, cutoff) - #Different backend solvers, default solver_backend = "applyexp" + #Different backend updaters, default updater_backend = "applyexp" ψ2_exponentiate_backend = tdvp( - H, +0.1im, ψ1; nsteps=1, cutoff, solver=exponentiate_solver + H, +0.1im, ψ1; nsteps=1, cutoff, updater=exponentiate_updater ) @test ψ2 ≈ ψ2_exponentiate_backend rtol = 1e-7 @@ -83,9 +83,9 @@ using Test ψ1 = tdvp(Hs, -0.1im, ψ0; nsteps=1, cutoff, nsite=1) - #Different backend solvers, default solver_backend = "applyexp" + #Different backend updaters, default updater_backend = "applyexp" ψ1_exponentiate_backend = tdvp( - Hs, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, solver=exponentiate_solver + Hs, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, updater=exponentiate_updater ) @test ψ1 ≈ ψ1_exponentiate_backend rtol = 1e-7 @@ -100,9 +100,9 @@ using Test # Time evolve backwards: ψ2 = tdvp(Hs, +0.1im, ψ1; nsteps=1, cutoff) - #Different backend solvers, default solver_backend = "applyexp" + #Different backend updaters, default updater_backend = "applyexp" ψ2_exponentiate_backend = tdvp( - Hs, +0.1im, ψ1; nsteps=1, cutoff, solver=exponentiate_solver + Hs, +0.1im, ψ1; nsteps=1, cutoff, updater=exponentiate_updater ) @test ψ2 ≈ ψ2_exponentiate_backend rtol = 1e-7 @@ -146,8 +146,8 @@ using Test # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - - @testset "Custom solver in TDVP" begin + #= + @testset_broken "Custom updater in TDVP" begin N = 10 cutoff = 1e-12 @@ -164,16 +164,18 @@ using Test ψ0 = random_mps(s; internal_inds_space=10) - function solver(psi0; (psi_ref!), (PH_ref!), time_step, kwargs...) - solver_kwargs = (; + function updater(psi0; (psi_ref!), (PH_ref!), kwargs...) + updater_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true ) + time_step=get(step_kwargs,:time_step,nothing) + @assert !isnothing(time_step) PH = PH_ref![] - psi, exp_info = exponentiate(PH, time_step, psi0; solver_kwargs...) + psi, exp_info = exponentiate(PH, time_step, psi0; updater_kwargs...) return psi, (; info=exp_info) end - ψ1 = tdvp(solver, H, -0.1im, ψ0; cutoff, nsite=1) + ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsite=1) @test norm(ψ1) ≈ 1.0 @@ -191,7 +193,7 @@ using Test # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - + =# @testset "Accuracy Test" begin N = 4 tau = 0.1 @@ -233,9 +235,10 @@ using Test psi; cutoff, normalize=false, - solver_tol=1e-12, - solver_maxiter=500, - solver_krylovdim=25, + updater_kwargs=(; + tol=1e-12, + maxiter=500, + krylovdim=25,) ) # TODO: What should `expect` output? Right now # it outputs a dictionary. @@ -247,10 +250,11 @@ using Test psi2; cutoff, normalize=false, - solver_tol=1e-12, - solver_maxiter=500, - solver_krylovdim=25, - solver=exponentiate_solver, + updater_kwargs=(; + tol=1e-12, + maxiter=500, + krylovdim=25,), + updater=exponentiate_updater, ) # TODO: What should `expect` output? Right now # it outputs a dictionary. @@ -314,7 +318,7 @@ using Test nsite = (step <= 3 ? 2 : 1) phi = tdvp( - H, -tau * im, phi; nsteps=1, cutoff, nsite, normalize=true, solver_krylovdim=15 + H, -tau * im, phi; nsteps=1, cutoff, nsite, normalize=true, updater_kwargs=(;krylovdim=15) ) Sz1[step] = real(expect("Sz", psi; vertices=[c])[c]) @@ -374,9 +378,9 @@ using Test for (step, t) in enumerate(trange) nsite = (step <= 10 ? 2 : 1) psi = tdvp( - H, -tau, psi; cutoff, nsite, reverse_step, normalize=true, solver_krylovdim=15 + H, -tau, psi; cutoff, nsite, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) ) - #Different backend solvers, default solver_backend = "applyexp" + #Different backend updaters, default updater_backend = "applyexp" psi2 = tdvp( H, -tau, @@ -385,8 +389,8 @@ using Test nsite, reverse_step, normalize=true, - solver_krylovdim=15, - solver=ITensorNetworks.exponentiate_solver, + updater_kwargs=(;krylovdim=15), + updater=ITensorNetworks.exponentiate_updater, ) end @@ -427,7 +431,7 @@ using Test get_info(; info) = info step_measure_sz(; psi) = expect("Sz", psi; vertices=[c])[c] step_measure_en(; psi) = real(inner(psi', H, psi)) - step_obs = Observer( + region_obs = Observer( "Sz" => step_measure_sz, "En" => step_measure_en, "info" => get_info ) @@ -440,16 +444,16 @@ using Test cutoff, normalize=false, (sweep_observer!)=sweep_obs, - (step_observer!)=step_obs, + (region_observer!)=region_obs, root_vertex=N, # defaults to 1, which breaks observer equality ) Sz2 = sweep_obs.Sz En2 = sweep_obs.En - Sz2_step = step_obs.Sz - En2_step = step_obs.En - infos = step_obs.info + Sz2_step = region_obs.Sz + En2_step = region_obs.En + infos = region_obs.info #@show sweep_obs #@show step_obs @@ -538,8 +542,8 @@ end # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - - @testset "Custom solver in TDVP" begin +#= + @testset "Custom updater in TDVP" begin cutoff = 1e-12 tooth_lengths = fill(2, 3) @@ -553,19 +557,19 @@ end ψ0 = normalize!(random_ttn(s; link_space=10)) - function solver(psi0; (psi_ref!), (PH_ref!), time_step, kwargs...) - solver_kwargs = (; + function updater(psi0; (psi_ref!), (PH_ref!), time_step, kwargs...) + updater_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true ) PH = PH_ref![] - psi, exp_info = exponentiate(PH, time_step, psi0; solver_kwargs...) + psi, exp_info = exponentiate(PH, time_step, psi0; updater_kwargs...) return psi, (; info=exp_info) end - ψ1 = tdvp(solver, H, -0.1im, ψ0; cutoff, nsite=1) + ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsite=1) - #@test ψ1 ≈ tdvp(solver, -0.1im, H, ψ0; cutoff, nsite=1) - #@test ψ1 ≈ tdvp(solver, H, ψ0, -0.1im; cutoff, nsite=1) + #@test ψ1 ≈ tdvp(updater, -0.1im, H, ψ0; cutoff, nsite=1) + #@test ψ1 ≈ tdvp(updater, H, ψ0, -0.1im; cutoff, nsite=1) @test norm(ψ1) ≈ 1.0 @@ -583,7 +587,7 @@ end # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - +=# @testset "Accuracy Test" begin tau = 0.1 ttotal = 1.0 @@ -620,9 +624,10 @@ end psi; cutoff, normalize=false, - solver_tol=1e-12, - solver_maxiter=500, - solver_krylovdim=25, + updater_kwargs=(; + tol=1e-12, + maxiter=500, + krylovdim=25,) ) push!(Sz_tdvp, real(expect("Sz", psi; vertices=[c])[c])) push!(Sz_exact, real(scalar(dag(prime(psix, s[c])) * Szc * psix))) @@ -680,7 +685,7 @@ end nsite = (step <= 3 ? 2 : 1) phi = tdvp( - H, -tau * im, phi; nsteps=1, cutoff, nsite, normalize=true, solver_krylovdim=15 + H, -tau * im, phi; nsteps=1, cutoff, nsite, normalize=true, updater_kwargs=(;krylovdim=15) ) Sz1[step] = real(expect("Sz", psi; vertices=[c])[c]) @@ -705,7 +710,7 @@ end time_step=-im * tau, cutoff, normalize=false, - (step_observer!)=obs, + (region_observer!)=obs, root_vertex=(3, 2), ) @@ -731,7 +736,7 @@ end for (step, t) in enumerate(trange) nsite = (step <= 10 ? 2 : 1) psi = tdvp( - H, -tau, psi; cutoff, nsite, reverse_step, normalize=true, solver_krylovdim=15 + H, -tau, psi; cutoff, nsite, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) ) end From db524f70dfee64a73061e0fe2264a16715b197e5 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 17 Jan 2024 18:02:27 -0500 Subject: [PATCH 05/32] Adapt eigsolve and dmrg to new interfaces, fix dmrg tests (and remove a broken one for now). --- src/solvers/eigsolve.jl | 61 ++++++++--------- src/solvers/exponentiate.jl | 4 +- src/treetensornetworks/solvers/dmrg.jl | 65 +++++++++++++------ .../test_solvers/test_dmrg.jl | 22 ++++--- 4 files changed, 90 insertions(+), 62 deletions(-) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index a1dfa39c..0703adaa 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -1,38 +1,39 @@ -function eigsolve_solver(; - solver_which_eigenvalue=:SR, #TODO: settle on pattern to pass solver kwargs +function eigsolve_updater( + init; + psi_ref!, + PH_ref!, + outputlevel, + which_sweep, + region_updates, + which_region_update, + region_kwargs, + updater_kwargs, + ) + + default_updater_kwargs=(; + solver_which_eigenvalue=:SR, ishermitian=true, - solver_tol=1e-14, - solver_krylovdim=3, - solver_maxiter=1, - solver_verbosity=0, -) - function solver( - init; - psi_ref!, - PH_ref!, - normalize, - region, - sweep_regions, - sweep_step, - sweep_kwargs..., - # slurp solver_kwargs? #TODO: homogenize how the solver kwargs are passed + tol=1e-14, + krylovdim=3, + maxiter=1, + outputlevel=0, + eager=false, ) - howmany = 1 - which = solver_which_eigenvalue - vals, vecs, info = eigsolve( - PH_ref![], + updater_kwargs=merge(default_updater_kwargs,updater_kwargs) #last collection has precedence + howmany=1 + which=updater_kwargs[:solver_which_eigenvalue] + vals, vecs, info = KrylovKit.eigsolve(PH_ref![], init, howmany, which; - ishermitian, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_verbosity, + ishermitian=updater_kwargs[:ishermitian], + tol=updater_kwargs[:tol], + krylovdim=updater_kwargs[:krylovdim], + maxiter=updater_kwargs[:maxiter], + verbosity=updater_kwargs[:outputlevel], + eager=updater_kwargs[:eager], ) - phi = vecs[1] - return phi, (; solver_info=info, energies=vals) - end - return solver + return vecs[1], (; info, energies=vals) end + diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index e70edc80..cf08c4d0 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -3,10 +3,10 @@ psi_ref!, PH_ref!, outputlevel, - which_sweep, # keep, change name + which_sweep, region_updates, which_region_update, - region_kwargs, # region_kwargs (timestep for solver) + region_kwargs, updater_kwargs, ) default_updater_kwargs=(; diff --git a/src/treetensornetworks/solvers/dmrg.jl b/src/treetensornetworks/solvers/dmrg.jl index 8f24fd63..be4c3039 100644 --- a/src/treetensornetworks/solvers/dmrg.jl +++ b/src/treetensornetworks/solvers/dmrg.jl @@ -1,34 +1,59 @@ """ Overload of `ITensors.dmrg`. """ + + +function dmrg_sweep( + nsite::Int, + graph::AbstractGraph; + root_vertex=default_root_vertex(graph), +) + return tdvp_sweep(2,nsite,Inf,graph;root_vertex,reverse_step=false) +end + + +function default_sweep_regions(nsite, graph::AbstractGraph; kwargs...) ###move this to a different file, algorithmic level idea + return vcat( + [ + half_sweep( + direction(half), + graph, + make_region; + nsite, + region_args=(; substep=half), + kwargs..., + ) for half in 1:2 + ]..., + ) +end + function dmrg( + updater, H, init::AbstractTTN; - solver_which_eigenvalue=:SR, - ishermitian=true, - solver_tol=1e-14, - solver_krylovdim=3, - solver_maxiter=1, - solver_verbosity=0, + nsweeps, #it makes sense to require this to be defined + nsite=2, + (sweep_observer!)=observer(), + root_vertex=default_root_vertex(init), + updater_kwargs=NamedTuple(;), kwargs..., -) - return alternating_update( - eigsolve_solver(; - solver_which_eigenvalue, - ishermitian, - solver_tol, - solver_krylovdim, - solver_maxiter, - solver_verbosity, - ), - H, - init; - reverse_step=false, - kwargs..., ) + region_updates = dmrg_sweep(nsite,init; root_vertex) + + psi = alternating_update( + updater, H, init; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... + ) + return psi +end + + +function dmrg(H, init::AbstractTTN; updater=eigsolve_updater, kwargs...) + return dmrg(updater, H, init; kwargs...) end """ Overload of `KrylovKit.eigsolve`. """ eigsolve(H, init::AbstractTTN; kwargs...) = dmrg(H, init; kwargs...) + + diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl index 939b4fee..1e523b28 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl @@ -30,21 +30,23 @@ using Observers psi_mps = MPS([psi[v] for v in 1:nv(psi)]) e2, psi2 = dmrg(H_mpo, psi_mps; nsweeps, maxdim, outputlevel=0) - psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1) + psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1)) @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2) # Alias for `ITensorNetworks.dmrg` psi = eigsolve( - H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1 + H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1) ) @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2) - # Test custom sweep regions + # Test custom sweep regions #BROKEN, ToDo: Make proper custom sweep regions for test + #= orig_E = inner(psi', H, psi) sweep_regions = [[1], [2], [3], [3], [2], [1]] psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_regions) new_E = inner(psi', H, psi) @test new_E ≈ orig_E + =# end @testset "Observers" begin @@ -67,20 +69,20 @@ end # # Make observers # - sweep(; sweep, kw...) = sweep + sweep(; which_sweep, kw...) = which_sweep sweep_observer! = observer(sweep) - region(; region, kw...) = region + region(; which_region_update, region_updates, kw...) = first(region_updates[which_region_update]) energy(; energies, kw...) = energies[1] - step_observer! = observer(region, sweep, energy) + region_observer! = observer(region, sweep, energy) - psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_observer!, step_observer!) + psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_observer!, region_observer!) # # Test out certain values # - @test step_observer![9, :region] == [2, 1] - @test step_observer![30, :energy] < -4.25 + @test region_observer![9, :region] == [2, 1] + @test region_observer![30, :energy] < -4.25 end @testset "Regression test: Arrays of Parameters" begin @@ -126,7 +128,7 @@ end sweeps = Sweeps(nsweeps) # number of sweeps is 5 maxdim!(sweeps, 10, 20, 40, 100) # gradually increase states kept cutoff!(sweeps, cutoff) - psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1) + psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1)) # Compare to `ITensors.MPO` version of `dmrg` linear_order = [4, 1, 2, 5, 3, 6] From d258ccc1307c231d149ffe0e14ec7c7512b7fdb2 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 17 Jan 2024 18:15:44 -0500 Subject: [PATCH 06/32] Remove default_sweep_regions from dmrg. --- src/treetensornetworks/solvers/dmrg.jl | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/src/treetensornetworks/solvers/dmrg.jl b/src/treetensornetworks/solvers/dmrg.jl index be4c3039..bf5c51ed 100644 --- a/src/treetensornetworks/solvers/dmrg.jl +++ b/src/treetensornetworks/solvers/dmrg.jl @@ -11,22 +11,6 @@ function dmrg_sweep( return tdvp_sweep(2,nsite,Inf,graph;root_vertex,reverse_step=false) end - -function default_sweep_regions(nsite, graph::AbstractGraph; kwargs...) ###move this to a different file, algorithmic level idea - return vcat( - [ - half_sweep( - direction(half), - graph, - make_region; - nsite, - region_args=(; substep=half), - kwargs..., - ) for half in 1:2 - ]..., - ) -end - function dmrg( updater, H, From 6ff116ac123a1968cec090216e4118c36e36ced1 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Wed, 17 Jan 2024 18:18:54 -0500 Subject: [PATCH 07/32] Format. --- src/solvers/eigsolve.jl | 49 ++++++----- src/solvers/exponentiate.jl | 81 ++++++++++--------- .../solvers/alternating_update.jl | 20 +++-- src/treetensornetworks/solvers/dmrg.jl | 14 +--- .../test_solvers/test_dmrg.jl | 13 ++- 5 files changed, 90 insertions(+), 87 deletions(-) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 0703adaa..058ca197 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -5,35 +5,34 @@ function eigsolve_updater( PH_ref!, outputlevel, which_sweep, - region_updates, + region_updates, which_region_update, region_kwargs, updater_kwargs, +) + default_updater_kwargs = (; + solver_which_eigenvalue=:SR, + ishermitian=true, + tol=1e-14, + krylovdim=3, + maxiter=1, + outputlevel=0, + eager=false, ) - - default_updater_kwargs=(; - solver_which_eigenvalue=:SR, - ishermitian=true, - tol=1e-14, - krylovdim=3, - maxiter=1, - outputlevel=0, - eager=false, + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence + howmany = 1 + which = updater_kwargs[:solver_which_eigenvalue] + vals, vecs, info = KrylovKit.eigsolve( + PH_ref![], + init, + howmany, + which; + ishermitian=updater_kwargs[:ishermitian], + tol=updater_kwargs[:tol], + krylovdim=updater_kwargs[:krylovdim], + maxiter=updater_kwargs[:maxiter], + verbosity=updater_kwargs[:outputlevel], + eager=updater_kwargs[:eager], ) - updater_kwargs=merge(default_updater_kwargs,updater_kwargs) #last collection has precedence - howmany=1 - which=updater_kwargs[:solver_which_eigenvalue] - vals, vecs, info = KrylovKit.eigsolve(PH_ref![], - init, - howmany, - which; - ishermitian=updater_kwargs[:ishermitian], - tol=updater_kwargs[:tol], - krylovdim=updater_kwargs[:krylovdim], - maxiter=updater_kwargs[:maxiter], - verbosity=updater_kwargs[:outputlevel], - eager=updater_kwargs[:eager], - ) return vecs[1], (; info, energies=vals) end - diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index cf08c4d0..c9682314 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -1,47 +1,48 @@ - function exponentiate_updater( - init; - psi_ref!, - PH_ref!, - outputlevel, - which_sweep, - region_updates, - which_region_update, - region_kwargs, - updater_kwargs, - ) - default_updater_kwargs=(; +function exponentiate_updater( + init; + psi_ref!, + PH_ref!, + outputlevel, + which_sweep, + region_updates, + which_region_update, + region_kwargs, + updater_kwargs, +) + default_updater_kwargs = (; krylovdim=30, #from here only solver kwargs maxiter=100, outputlevel=0, tol=1E-12, ishermitian=true, issymmetric=true, - eager=true - ) - updater_kwargs=merge(default_updater_kwargs,updater_kwargs) #last collection has precedence - #H=copy(PH_ref![]) - H = PH_ref![] ###since we are not changing H we don't need the copy - # let's test whether given region and sweep regions we can find out what the previous and next region were - # this will be needed in subspace expansion - #@show step_kwargs - substep=get(region_kwargs,:substep,nothing) - time_step=get(region_kwargs,:time_step,nothing) - @assert !isnothing(time_step) && !isnothing(substep) - region_ind = which_region_update - next_region = region_ind == length(region_updates) ? nothing : first(region_updates[region_ind + 1]) - previous_region = region_ind == 1 ? nothing : first(region_updates[region_ind - 1]) + eager=true, + ) + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence + #H=copy(PH_ref![]) + H = PH_ref![] ###since we are not changing H we don't need the copy + # let's test whether given region and sweep regions we can find out what the previous and next region were + # this will be needed in subspace expansion + #@show step_kwargs + substep = get(region_kwargs, :substep, nothing) + time_step = get(region_kwargs, :time_step, nothing) + @assert !isnothing(time_step) && !isnothing(substep) + region_ind = which_region_update + next_region = + region_ind == length(region_updates) ? nothing : first(region_updates[region_ind + 1]) + previous_region = region_ind == 1 ? nothing : first(region_updates[region_ind - 1]) - phi, exp_info = KrylovKit.exponentiate( - H, - time_step, - init; - ishermitian=updater_kwargs[:ishermitian], - issymmetric=updater_kwargs[:issymmetric], - tol=updater_kwargs[:tol], - krylovdim=updater_kwargs[:krylovdim], - maxiter=updater_kwargs[:maxiter], - verbosity=updater_kwargs[:outputlevel], - eager=updater_kwargs[:eager], - ) - return phi, (; info=exp_info) - end + phi, exp_info = KrylovKit.exponentiate( + H, + time_step, + init; + ishermitian=updater_kwargs[:ishermitian], + issymmetric=updater_kwargs[:issymmetric], + tol=updater_kwargs[:tol], + krylovdim=updater_kwargs[:krylovdim], + maxiter=updater_kwargs[:maxiter], + verbosity=updater_kwargs[:outputlevel], + eager=updater_kwargs[:eager], + ) + return phi, (; info=exp_info) +end diff --git a/src/treetensornetworks/solvers/alternating_update.jl b/src/treetensornetworks/solvers/alternating_update.jl index 04b5ce37..c7534003 100644 --- a/src/treetensornetworks/solvers/alternating_update.jl +++ b/src/treetensornetworks/solvers/alternating_update.jl @@ -56,7 +56,8 @@ function alternating_update( insert_function!(sweep_observer!, "sweep_printer" => sweep_printer) # FIX THIS for which_sweep in 1:nsweeps - if !isnothing(write_when_maxdim_exceeds) && maxdim[which_sweep] > write_when_maxdim_exceeds + if !isnothing(write_when_maxdim_exceeds) && + maxdim[which_sweep] > write_when_maxdim_exceeds if outputlevel >= 2 println( "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim[which_sweep] = $(maxdim[which_sweep]), writing environment tensors to disk", @@ -64,7 +65,7 @@ function alternating_update( end PH = disk(PH) end - + sw_time = @elapsed begin psi, PH = sweep_update( updater, @@ -73,10 +74,11 @@ function alternating_update( outputlevel, which_sweep, sweep_params=(; - maxdim=maxdim[which_sweep], - mindim=mindim[which_sweep], - cutoff=cutoff[which_sweep], - noise=noise[which_sweep]), + maxdim=maxdim[which_sweep], + mindim=mindim[which_sweep], + cutoff=cutoff[which_sweep], + noise=noise[which_sweep], + ), updater_kwargs, kwargs..., ) @@ -86,7 +88,7 @@ function alternating_update( checkdone(; psi, which_sweep, outputlevel, kwargs...) && break end - select!(sweep_observer!, Observers.DataFrames.Not("sweep_printer")) + select!(sweep_observer!, Observers.DataFrames.Not("sweep_printer")) return psi end @@ -119,7 +121,9 @@ each step of the algorithm when optimizing the MPS. Returns: * `psi::MPS` - time-evolved MPS """ -function alternating_update(updater, Hs::Vector{<:AbstractTTN}, psi0::AbstractTTN; kwargs...) +function alternating_update( + updater, Hs::Vector{<:AbstractTTN}, psi0::AbstractTTN; kwargs... +) for H in Hs check_hascommoninds(siteinds, H, psi0) check_hascommoninds(siteinds, H, psi0') diff --git a/src/treetensornetworks/solvers/dmrg.jl b/src/treetensornetworks/solvers/dmrg.jl index bf5c51ed..1e9ae256 100644 --- a/src/treetensornetworks/solvers/dmrg.jl +++ b/src/treetensornetworks/solvers/dmrg.jl @@ -2,13 +2,10 @@ Overload of `ITensors.dmrg`. """ - function dmrg_sweep( - nsite::Int, - graph::AbstractGraph; - root_vertex=default_root_vertex(graph), + nsite::Int, graph::AbstractGraph; root_vertex=default_root_vertex(graph) ) - return tdvp_sweep(2,nsite,Inf,graph;root_vertex,reverse_step=false) + return tdvp_sweep(2, nsite, Inf, graph; root_vertex, reverse_step=false) end function dmrg( @@ -21,8 +18,8 @@ function dmrg( root_vertex=default_root_vertex(init), updater_kwargs=NamedTuple(;), kwargs..., - ) - region_updates = dmrg_sweep(nsite,init; root_vertex) +) + region_updates = dmrg_sweep(nsite, init; root_vertex) psi = alternating_update( updater, H, init; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... @@ -30,7 +27,6 @@ function dmrg( return psi end - function dmrg(H, init::AbstractTTN; updater=eigsolve_updater, kwargs...) return dmrg(updater, H, init; kwargs...) end @@ -39,5 +35,3 @@ end Overload of `KrylovKit.eigsolve`. """ eigsolve(H, init::AbstractTTN; kwargs...) = dmrg(H, init; kwargs...) - - diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl index 1e523b28..7806e1d6 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl @@ -30,12 +30,14 @@ using Observers psi_mps = MPS([psi[v] for v in 1:nv(psi)]) e2, psi2 = dmrg(H_mpo, psi_mps; nsweeps, maxdim, outputlevel=0) - psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1)) + psi = dmrg( + H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(; krylovdim=3, maxiter=1) + ) @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2) # Alias for `ITensorNetworks.dmrg` psi = eigsolve( - H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1) + H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(; krylovdim=3, maxiter=1) ) @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2) @@ -72,7 +74,8 @@ end sweep(; which_sweep, kw...) = which_sweep sweep_observer! = observer(sweep) - region(; which_region_update, region_updates, kw...) = first(region_updates[which_region_update]) + region(; which_region_update, region_updates, kw...) = + first(region_updates[which_region_update]) energy(; energies, kw...) = energies[1] region_observer! = observer(region, sweep, energy) @@ -128,7 +131,9 @@ end sweeps = Sweeps(nsweeps) # number of sweeps is 5 maxdim!(sweeps, 10, 20, 40, 100) # gradually increase states kept cutoff!(sweeps, cutoff) - psi = dmrg(H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(;krylovdim=3, maxiter=1)) + psi = dmrg( + H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(; krylovdim=3, maxiter=1) + ) # Compare to `ITensors.MPO` version of `dmrg` linear_order = [4, 1, 2, 5, 3, 6] From 0ac2f2d0bf8fe59e9abad43e43c594230a3cc531 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 13:37:10 -0500 Subject: [PATCH 08/32] rename region_[update]_printer and nsite[s] --- src/treetensornetworks/solvers/dmrg.jl | 10 +++-- src/treetensornetworks/solvers/tdvp.jl | 8 ++-- .../solvers/tree_sweeping.jl | 8 ++-- src/treetensornetworks/solvers/update_step.jl | 12 +++--- .../test_solvers/test_dmrg.jl | 14 +++---- .../test_solvers/test_tdvp.jl | 40 +++++++++---------- 6 files changed, 47 insertions(+), 45 deletions(-) diff --git a/src/treetensornetworks/solvers/dmrg.jl b/src/treetensornetworks/solvers/dmrg.jl index 1e9ae256..40f83aad 100644 --- a/src/treetensornetworks/solvers/dmrg.jl +++ b/src/treetensornetworks/solvers/dmrg.jl @@ -3,9 +3,11 @@ Overload of `ITensors.dmrg`. """ function dmrg_sweep( - nsite::Int, graph::AbstractGraph; root_vertex=default_root_vertex(graph) + nsites::Int, graph::AbstractGraph; root_vertex=default_root_vertex(graph) ) - return tdvp_sweep(2, nsite, Inf, graph; root_vertex, reverse_step=false) + order = 2 + time_step = Inf + return tdvp_sweep(order, nsites, time_step, graph; root_vertex, reverse_step=false) end function dmrg( @@ -13,13 +15,13 @@ function dmrg( H, init::AbstractTTN; nsweeps, #it makes sense to require this to be defined - nsite=2, + nsites=2, (sweep_observer!)=observer(), root_vertex=default_root_vertex(init), updater_kwargs=NamedTuple(;), kwargs..., ) - region_updates = dmrg_sweep(nsite, init; root_vertex) + region_updates = dmrg_sweep(nsites, init; root_vertex) psi = alternating_update( updater, H, init; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... diff --git a/src/treetensornetworks/solvers/tdvp.jl b/src/treetensornetworks/solvers/tdvp.jl index badac5b5..1165632a 100644 --- a/src/treetensornetworks/solvers/tdvp.jl +++ b/src/treetensornetworks/solvers/tdvp.jl @@ -32,7 +32,7 @@ end function tdvp_sweep( order::Int, - nsite::Int, + nsites::Int, time_step::Number, graph::AbstractGraph; root_vertex=default_root_vertex(graph), @@ -46,7 +46,7 @@ function tdvp_sweep( graph, make_region; root_vertex, - nsite, + nsites, region_args=(; substep, time_step=sub_time_step), reverse_args=(; substep, time_step=-sub_time_step), reverse_step, @@ -62,7 +62,7 @@ function tdvp( t::Number, init::AbstractTTN; time_step::Number=t, - nsite=2, + nsites=2, nsteps=nothing, order::Integer=2, (sweep_observer!)=observer(), @@ -72,7 +72,7 @@ function tdvp( kwargs..., ) nsweeps = _compute_nsweeps(nsteps, t, time_step, order) - region_updates = tdvp_sweep(order, nsite, time_step, init; root_vertex, reverse_step) + region_updates = tdvp_sweep(order, nsites, time_step, init; root_vertex, reverse_step) function sweep_time_printer(; outputlevel, which_sweep, kwargs...) if outputlevel >= 1 diff --git a/src/treetensornetworks/solvers/tree_sweeping.jl b/src/treetensornetworks/solvers/tree_sweeping.jl index 66171e68..99375ba9 100644 --- a/src/treetensornetworks/solvers/tree_sweeping.jl +++ b/src/treetensornetworks/solvers/tree_sweeping.jl @@ -3,12 +3,12 @@ direction(step_number) = isodd(step_number) ? Base.Forward : Base.Reverse function make_region( edge; last_edge=false, - nsite=1, + nsites=1, region_args=(;), reverse_args=region_args, reverse_step=false, ) - if nsite == 1 + if nsites == 1 site = ([src(edge)], region_args) bond = (edge, reverse_args) region = reverse_step ? (site, bond) : (site,) @@ -17,7 +17,7 @@ function make_region( else return region end - elseif nsite == 2 + elseif nsites == 2 sites_two = ([src(edge), dst(edge)], region_args) sites_one = ([dst(edge)], reverse_args) region = reverse_step ? (sites_two, sites_one) : (sites_two,) @@ -27,7 +27,7 @@ function make_region( return region end else - error("nsite=$nsite not supported in alternating_update / update_step") + error("nsites=$nsites not supported in alternating_update / update_step") end end diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index dfc5ce62..265bace6 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -1,12 +1,12 @@ -function default_sweep_regions(nsite, graph::AbstractGraph; kwargs...) ###move this to a different file, algorithmic level idea +function default_sweep_regions(nsites, graph::AbstractGraph; kwargs...) ###move this to a different file, algorithmic level idea return vcat( [ half_sweep( direction(half), graph, make_region; - nsite, + nsites, region_args=(; half_sweep=half), kwargs..., ) for half in 1:2 @@ -14,7 +14,7 @@ function default_sweep_regions(nsite, graph::AbstractGraph; kwargs...) ###move ) end -function region_printer(; +function region_update_printer(; cutoff, maxdim, mindim, outputlevel::Int=0, psi, region_updates, spec, which_region_update, which_sweep,kwargs... ) if outputlevel >= 2 @@ -42,14 +42,14 @@ function sweep_update( psi::AbstractTTN; normalize::Bool=false, # ToDo: think about where to put the default, probably this default is best defined at algorithmic level outputlevel, - step_printer=step_printer, + region_update_printer=region_update_printer, (region_observer!)=observer(), # ToDo: change name to region_observer! ? which_sweep::Int, sweep_params::NamedTuple, region_updates,# =default_sweep_regions(nsite, psi; reverse_step), #move default up to algorithmic level updater_kwargs, ) - insert_function!(region_observer!, "region_printer" => region_printer) #ToDo fix this + insert_function!(region_observer!, "region_update_printer" => region_update_printer) #ToDo fix this # Append empty namedtuple to each element if not already present # (Needed to handle user-provided region_updates) @@ -80,7 +80,7 @@ function sweep_update( ) end - select!(region_observer!, Observers.DataFrames.Not("region_printer")) # remove step_printer #todo fix this + select!(region_observer!, Observers.DataFrames.Not("region_update_printer")) # remove update_printer # Just to be sure: normalize && normalize!(psi) diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl index 7806e1d6..5d6392bb 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl @@ -5,7 +5,7 @@ using Random using Test using Observers -@testset "MPS DMRG" for nsite in [1, 2] +@testset "MPS DMRG" for nsites in [1, 2] N = 10 cutoff = 1e-12 @@ -31,13 +31,13 @@ using Observers e2, psi2 = dmrg(H_mpo, psi_mps; nsweeps, maxdim, outputlevel=0) psi = dmrg( - H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(; krylovdim=3, maxiter=1) + H, psi; nsweeps, maxdim, cutoff, nsites, updater_kwargs=(; krylovdim=3, maxiter=1) ) @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2) # Alias for `ITensorNetworks.dmrg` psi = eigsolve( - H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(; krylovdim=3, maxiter=1) + H, psi; nsweeps, maxdim, cutoff, nsites, updater_kwargs=(; krylovdim=3, maxiter=1) ) @test inner(psi', H, psi) ≈ inner(psi2', H_mpo, psi2) @@ -113,7 +113,7 @@ end psi = dmrg(H, psi; nsweeps, maxdim, cutoff) end -@testset "Tree DMRG" for nsite in [1, 2] +@testset "Tree DMRG" for nsites in [1, 2] cutoff = 1e-12 tooth_lengths = fill(2, 3) @@ -132,7 +132,7 @@ end maxdim!(sweeps, 10, 20, 40, 100) # gradually increase states kept cutoff!(sweeps, cutoff) psi = dmrg( - H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(; krylovdim=3, maxiter=1) + H, psi; nsweeps, maxdim, cutoff, nsites, updater_kwargs=(; krylovdim=3, maxiter=1) ) # Compare to `ITensors.MPO` version of `dmrg` @@ -148,7 +148,7 @@ end @testset "Regression test: tree truncation" begin maxdim = 4 - nsite = 2 + nsites = 2 nsweeps = 10 c = named_comb_tree((3, 2)) @@ -156,7 +156,7 @@ end os = ITensorNetworks.heisenberg(c) H = TTN(os, s) psi = random_ttn(s; link_space=5) - psi = dmrg(H, psi; nsweeps, maxdim, nsite) + psi = dmrg(H, psi; nsweeps, maxdim, nsites) @test all(edge_data(linkdims(psi)) .<= maxdim) end diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index f8fadb1b..60642a9b 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -25,14 +25,14 @@ exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best han ψ0 = random_mps(s; internal_inds_space=10) # Time evolve forward: - ψ1 = tdvp(H, -0.1im, ψ0; nsteps=1, cutoff, nsite=1) + ψ1 = tdvp(H, -0.1im, ψ0; nsteps=1, cutoff, nsites=1) # #TODO: exponentiate is now the default, so switch this to applyexp # #Different backend updaters, default updater_backend = "applyexp" ψ1_exponentiate_backend = tdvp( - H, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, updater=exponentiate_updater + H, -0.1im, ψ0; nsteps=1, cutoff, nsites=1, updater=exponentiate_updater ) @test ψ1 ≈ ψ1_exponentiate_backend rtol = 1e-7 @@ -81,11 +81,11 @@ exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best han ψ0 = random_mps(s; internal_inds_space=10) - ψ1 = tdvp(Hs, -0.1im, ψ0; nsteps=1, cutoff, nsite=1) + ψ1 = tdvp(Hs, -0.1im, ψ0; nsteps=1, cutoff, nsites=1) #Different backend updaters, default updater_backend = "applyexp" ψ1_exponentiate_backend = tdvp( - Hs, -0.1im, ψ0; nsteps=1, cutoff, nsite=1, updater=exponentiate_updater + Hs, -0.1im, ψ0; nsteps=1, cutoff, nsites=1, updater=exponentiate_updater ) @test ψ1 ≈ ψ1_exponentiate_backend rtol = 1e-7 @@ -131,7 +131,7 @@ exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best han ψ0 = random_mps(s; internal_inds_space=10) # Time evolve forward: - ψ1 = tdvp(H, -0.1im, ψ0; time_step=-0.05im, order, cutoff, nsite=1) + ψ1 = tdvp(H, -0.1im, ψ0; time_step=-0.05im, order, cutoff, nsites=1) @test norm(ψ1) ≈ 1.0 @@ -175,7 +175,7 @@ exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best han return psi, (; info=exp_info) end - ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsite=1) + ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsites=1) @test norm(ψ1) ≈ 1.0 @@ -316,9 +316,9 @@ exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best han psi = apply(gates, psi; cutoff) #normalize!(psi) - nsite = (step <= 3 ? 2 : 1) + nsites = (step <= 3 ? 2 : 1) phi = tdvp( - H, -tau * im, phi; nsteps=1, cutoff, nsite, normalize=true, updater_kwargs=(;krylovdim=15) + H, -tau * im, phi; nsteps=1, cutoff, nsites, normalize=true, updater_kwargs=(;krylovdim=15) ) Sz1[step] = real(expect("Sz", psi; vertices=[c])[c]) @@ -376,9 +376,9 @@ exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best han psi2 = deepcopy(psi) trange = 0.0:tau:ttotal for (step, t) in enumerate(trange) - nsite = (step <= 10 ? 2 : 1) + nsites = (step <= 10 ? 2 : 1) psi = tdvp( - H, -tau, psi; cutoff, nsite, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) + H, -tau, psi; cutoff, nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) ) #Different backend updaters, default updater_backend = "applyexp" psi2 = tdvp( @@ -386,7 +386,7 @@ exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best han -tau, psi2; cutoff, - nsite, + nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15), @@ -482,7 +482,7 @@ end ψ0 = normalize!(random_ttn(s; link_space=10)) # Time evolve forward: - ψ1 = tdvp(H, -0.1im, ψ0; nsteps=1, cutoff, nsite=1) + ψ1 = tdvp(H, -0.1im, ψ0; nsteps=1, cutoff, nsites=1) @test norm(ψ1) ≈ 1.0 @@ -524,7 +524,7 @@ end ψ0 = normalize!(random_ttn(s; link_space=10)) - ψ1 = tdvp(Hs, -0.1im, ψ0; nsteps=1, cutoff, nsite=1) + ψ1 = tdvp(Hs, -0.1im, ψ0; nsteps=1, cutoff, nsites=1) @test norm(ψ1) ≈ 1.0 @@ -566,10 +566,10 @@ end return psi, (; info=exp_info) end - ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsite=1) + ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsites=1) - #@test ψ1 ≈ tdvp(updater, -0.1im, H, ψ0; cutoff, nsite=1) - #@test ψ1 ≈ tdvp(updater, H, ψ0, -0.1im; cutoff, nsite=1) + #@test ψ1 ≈ tdvp(updater, -0.1im, H, ψ0; cutoff, nsites=1) + #@test ψ1 ≈ tdvp(updater, H, ψ0, -0.1im; cutoff, nsites=1) @test norm(ψ1) ≈ 1.0 @@ -683,9 +683,9 @@ end psi = apply(gates, psi; cutoff, maxdim) #normalize!(psi) - nsite = (step <= 3 ? 2 : 1) + nsites = (step <= 3 ? 2 : 1) phi = tdvp( - H, -tau * im, phi; nsteps=1, cutoff, nsite, normalize=true, updater_kwargs=(;krylovdim=15) + H, -tau * im, phi; nsteps=1, cutoff, nsites, normalize=true, updater_kwargs=(;krylovdim=15) ) Sz1[step] = real(expect("Sz", psi; vertices=[c])[c]) @@ -734,9 +734,9 @@ end trange = 0.0:tau:ttotal for (step, t) in enumerate(trange) - nsite = (step <= 10 ? 2 : 1) + nsites = (step <= 10 ? 2 : 1) psi = tdvp( - H, -tau, psi; cutoff, nsite, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) + H, -tau, psi; cutoff, nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) ) end From c02662c743fcfb25f1311167266cba47aa9a19c0 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 13:50:37 -0500 Subject: [PATCH 09/32] Remove applyexp.jl --- src/treetensornetworks/solvers/applyexp.jl | 128 --------------------- 1 file changed, 128 deletions(-) delete mode 100644 src/treetensornetworks/solvers/applyexp.jl diff --git a/src/treetensornetworks/solvers/applyexp.jl b/src/treetensornetworks/solvers/applyexp.jl deleted file mode 100644 index 6e84036e..00000000 --- a/src/treetensornetworks/solvers/applyexp.jl +++ /dev/null @@ -1,128 +0,0 @@ -# -# To Do: -# - implement assembleLanczosVectors -# - check slice ranges - change end value by 1? -# - -function assemble_lanczos_vecs(lanczos_vectors, linear_comb, norm) - #if length(lanczos_vectors) != length(linear_comb) - # @show length(lanczos_vectors) - # @show length(linear_comb) - #end - xt = norm * linear_comb[1] * lanczos_vectors[1] - for i in 2:length(lanczos_vectors) - xt += norm * linear_comb[i] * lanczos_vectors[i] - end - return xt -end - -struct ApplyExpInfo - numops::Int - converged::Int -end - -function applyexp(H, tau::Number, x0; maxiter=30, tol=1e-12, outputlevel=0, normcutoff=1e-7) - # Initialize Lanczos vectors - v1 = copy(x0) - nrm = norm(v1) - v1 /= nrm - lanczos_vectors = [v1] - - ElT = promote_type(typeof(tau), eltype(x0)) - - bigTmat = zeros(ElT, maxiter + 3, maxiter + 3) - - nmatvec = 0 - - v0 = nothing - beta = 0.0 - for iter in 1:maxiter - tmat_size = iter + 1 - - # Matrix-vector multiplication - w = H(v1) - nmatvec += 1 - - avnorm = norm(w) - alpha = dot(w, v1) - - bigTmat[iter, iter] = alpha - - w -= alpha * v1 - if iter > 1 - w -= beta * v0 - end - v0 = copy(v1) - - beta = norm(w) - - # check for Lanczos sequence exhaustion - if abs(beta) < beta_tol - # Assemble the time evolved state - tmat = bigTmat[1:tmat_size, 1:tmat_size] - tmat_exp = exp(tau * tmat) - linear_comb = tmat_exp[:, 1] - xt = assemble_lanczos_vecs(lanczos_vectors, linear_comb, nrm) - return xt, ApplyExpInfo(nmatvec, 1) - end - - # update next lanczos vector - v1 = copy(w) - v1 /= beta - push!(lanczos_vectors, v1) - bigTmat[iter + 1, iter] = beta - bigTmat[iter, iter + 1] = beta - - # Convergence check - if iter > 0 - # Prepare extended T-matrix for exponentiation - tmat_ext_size = tmat_size + 2 - tmat_ext = bigTmat[1:tmat_ext_size, 1:tmat_ext_size] - - tmat_ext[tmat_size - 1, tmat_size] = 0.0 - tmat_ext[tmat_size + 1, tmat_size] = 1.0 - - # Exponentiate extended T-matrix - tmat_ext_exp = exp(tau * tmat_ext) - - ϕ1 = abs(nrm * tmat_ext_exp[tmat_size, 1]) - ϕ2 = abs(nrm * tmat_ext_exp[tmat_size + 1, 1] * avnorm) - - if ϕ1 > 10 * ϕ2 - error = ϕ2 - elseif (ϕ1 > ϕ2) - error = (ϕ1 * ϕ2) / (ϕ1 - ϕ2) - else - error = ϕ1 - end - - if outputlevel >= 3 - @printf(" Iteration: %d, Error: %.2E\n", iter, error) - end - - if ((error < tol) || (iter == maxiter)) - converged = 1 - if (iter == maxiter) - println("warning: applyexp not converged in $maxiter steps") - converged = 0 - end - - # Assemble the time evolved state - linear_comb = tmat_ext_exp[:, 1] - xt = assemble_lanczos_vecs(lanczos_vectors, linear_comb, nrm) - - if outputlevel >= 3 - println(" Number of iterations: $iter") - end - - return xt, ApplyExpInfo(nmatvec, converged) - end - end # end convergence test - end # iter - - if outputlevel >= 0 - println("In applyexp, number of matrix-vector multiplies: ", nmatvec) - end - - return x0 -end From 1a8ff330ede03b4162ea2c7a7e26d7748c8ab35d Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 13:51:36 -0500 Subject: [PATCH 10/32] Fix imports/namespaces. --- src/imports.jl | 2 +- src/solvers/eigsolve.jl | 9 ++++----- src/solvers/exponentiate.jl | 2 +- test/test_treetensornetworks/test_solvers/test_dmrg.jl | 2 +- test/test_treetensornetworks/test_solvers/test_tdvp.jl | 2 +- 5 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/imports.jl b/src/imports.jl index 71406878..28bf0a88 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -34,7 +34,7 @@ import .DataGraphs: import Graphs: SimpleGraph, is_directed, weights -import KrylovKit: eigsolve, linsolve +import KrylovKit: eigsolve, linsolve, exponentiate import LinearAlgebra: factorize, normalize, normalize!, qr, svd diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 058ca197..6af93da0 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -1,4 +1,3 @@ - function eigsolve_updater( init; psi_ref!, @@ -11,7 +10,7 @@ function eigsolve_updater( updater_kwargs, ) default_updater_kwargs = (; - solver_which_eigenvalue=:SR, + which_eigenvalue=:SR, ishermitian=true, tol=1e-14, krylovdim=3, @@ -21,8 +20,8 @@ function eigsolve_updater( ) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence howmany = 1 - which = updater_kwargs[:solver_which_eigenvalue] - vals, vecs, info = KrylovKit.eigsolve( + which = updater_kwargs[:which_eigenvalue] + vals, vecs, info = eigsolve( PH_ref![], init, howmany, @@ -34,5 +33,5 @@ function eigsolve_updater( verbosity=updater_kwargs[:outputlevel], eager=updater_kwargs[:eager], ) - return vecs[1], (; info, energies=vals) + return vecs[1], (; info, eigvals=vals) end diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index c9682314..cccd8b21 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -32,7 +32,7 @@ function exponentiate_updater( region_ind == length(region_updates) ? nothing : first(region_updates[region_ind + 1]) previous_region = region_ind == 1 ? nothing : first(region_updates[region_ind - 1]) - phi, exp_info = KrylovKit.exponentiate( + phi, exp_info = exponentiate( H, time_step, init; diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl index 5d6392bb..92a318d7 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl @@ -76,7 +76,7 @@ end region(; which_region_update, region_updates, kw...) = first(region_updates[which_region_update]) - energy(; energies, kw...) = energies[1] + energy(; eigvals, kw...) = eigvals[1] region_observer! = observer(region, sweep, energy) psi = dmrg(H, psi; nsweeps, maxdim, cutoff, sweep_observer!, region_observer!) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index 60642a9b..3d6ca803 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -1,10 +1,10 @@ using ITensors using ITensorNetworks +using ITensorNetworks: exponentiate_updater using KrylovKit: exponentiate using Observers using Random using Test -exponentiate_updater=ITensorNetworks.exponentiate_updater #ToDo: how to best handle importing etc. @testset "MPS TDVP" begin @testset "Basic TDVP" begin From 80620b85c9518aae84b4b08c1a31439245914535 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 13:52:55 -0500 Subject: [PATCH 11/32] Add ToDo regarding testing applyexp in test_tdvp.jl --- test/test_treetensornetworks/test_solvers/test_tdvp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index 3d6ca803..496d4b21 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -28,7 +28,7 @@ using Test ψ1 = tdvp(H, -0.1im, ψ0; nsteps=1, cutoff, nsites=1) # - #TODO: exponentiate is now the default, so switch this to applyexp + #TODO: applyexp is removed, no need to test a different backend here and below for now. # #Different backend updaters, default updater_backend = "applyexp" ψ1_exponentiate_backend = tdvp( From 4c33d656deed877a3a2b012640318efb011ed79f Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 13:55:18 -0500 Subject: [PATCH 12/32] Define sweep_params outside of function call. --- .../solvers/alternating_update.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/treetensornetworks/solvers/alternating_update.jl b/src/treetensornetworks/solvers/alternating_update.jl index c7534003..37f445f3 100644 --- a/src/treetensornetworks/solvers/alternating_update.jl +++ b/src/treetensornetworks/solvers/alternating_update.jl @@ -65,7 +65,12 @@ function alternating_update( end PH = disk(PH) end - + sweep_params=(; + maxdim=maxdim[which_sweep], + mindim=mindim[which_sweep], + cutoff=cutoff[which_sweep], + noise=noise[which_sweep], + ) sw_time = @elapsed begin psi, PH = sweep_update( updater, @@ -73,12 +78,7 @@ function alternating_update( psi; outputlevel, which_sweep, - sweep_params=(; - maxdim=maxdim[which_sweep], - mindim=mindim[which_sweep], - cutoff=cutoff[which_sweep], - noise=noise[which_sweep], - ), + sweep_params, updater_kwargs, kwargs..., ) From 114ef3eb3a844635e2a53505c04f569e93b13983 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 13:58:32 -0500 Subject: [PATCH 13/32] Change NamedTuple access pattern. --- src/solvers/eigsolve.jl | 14 +++++++------- src/solvers/exponentiate.jl | 14 +++++++------- src/treetensornetworks/solvers/update_step.jl | 6 +++--- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 6af93da0..c0830100 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -20,18 +20,18 @@ function eigsolve_updater( ) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence howmany = 1 - which = updater_kwargs[:which_eigenvalue] + which = updater_kwargs.which_eigenvalue vals, vecs, info = eigsolve( PH_ref![], init, howmany, which; - ishermitian=updater_kwargs[:ishermitian], - tol=updater_kwargs[:tol], - krylovdim=updater_kwargs[:krylovdim], - maxiter=updater_kwargs[:maxiter], - verbosity=updater_kwargs[:outputlevel], - eager=updater_kwargs[:eager], + ishermitian=updater_kwargs.ishermitian, + tol=updater_kwargs.tol, + krylovdim=updater_kwargs.krylovdim, + maxiter=updater_kwargs.maxiter, + verbosity=updater_kwargs.outputlevel, + eager=updater_kwargs.eager, ) return vecs[1], (; info, eigvals=vals) end diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index cccd8b21..23609ae7 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -36,13 +36,13 @@ function exponentiate_updater( H, time_step, init; - ishermitian=updater_kwargs[:ishermitian], - issymmetric=updater_kwargs[:issymmetric], - tol=updater_kwargs[:tol], - krylovdim=updater_kwargs[:krylovdim], - maxiter=updater_kwargs[:maxiter], - verbosity=updater_kwargs[:outputlevel], - eager=updater_kwargs[:eager], + ishermitian=updater_kwargs.ishermitian, + issymmetric=updater_kwargs.issymmetric, + tol=updater_kwargs.tol, + krylovdim=updater_kwargs.krylovdim, + maxiter=updater_kwargs.maxiter, + verbosity=updater_kwargs.outputlevel, + eager=updater_kwargs.eager, ) return phi, (; info=exp_info) end diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index 265bace6..c459e79f 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -212,9 +212,9 @@ function region_update( psi, spec = insert_local_tensor( psi, phi, region; eigen_perturbation=drho, ortho, normalize, - maxdim=region_kwargs[:maxdim], - mindim=region_kwargs[:mindim], - cutoff=region_kwargs[:cutoff] + maxdim=region_kwargs.maxdim, + mindim=region_kwargs.mindim, + cutoff=region_kwargs.cutoff ) update!( From b81758791ba56991be5a98e89602854462ca1106 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 14:01:27 -0500 Subject: [PATCH 14/32] Change NamedTuple(;) to (;). --- src/treetensornetworks/solvers/dmrg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/treetensornetworks/solvers/dmrg.jl b/src/treetensornetworks/solvers/dmrg.jl index 40f83aad..b5525277 100644 --- a/src/treetensornetworks/solvers/dmrg.jl +++ b/src/treetensornetworks/solvers/dmrg.jl @@ -18,7 +18,7 @@ function dmrg( nsites=2, (sweep_observer!)=observer(), root_vertex=default_root_vertex(init), - updater_kwargs=NamedTuple(;), + updater_kwargs=(;), kwargs..., ) region_updates = dmrg_sweep(nsites, init; root_vertex) From 723fae90338d4d15a0e4996bf542b59f78dcc3ef Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 14:02:02 -0500 Subject: [PATCH 15/32] (;) also for tdvp. --- src/treetensornetworks/solvers/tdvp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/treetensornetworks/solvers/tdvp.jl b/src/treetensornetworks/solvers/tdvp.jl index 1165632a..130c92e5 100644 --- a/src/treetensornetworks/solvers/tdvp.jl +++ b/src/treetensornetworks/solvers/tdvp.jl @@ -68,7 +68,7 @@ function tdvp( (sweep_observer!)=observer(), root_vertex=default_root_vertex(init), reverse_step=true, - updater_kwargs=NamedTuple(;), + updater_kwargs=(;), kwargs..., ) nsweeps = _compute_nsweeps(nsteps, t, time_step, order) From 1e40a79ce4caf44afc84355396ff801ed3b83c27 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 14:04:01 -0500 Subject: [PATCH 16/32] Remove second applyexp.jl file. --- src/solvers/applyexp.jl | 25 ------------------------- 1 file changed, 25 deletions(-) delete mode 100644 src/solvers/applyexp.jl diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl deleted file mode 100644 index 222bb9b6..00000000 --- a/src/solvers/applyexp.jl +++ /dev/null @@ -1,25 +0,0 @@ -function applyexp_solver() - function solver( - init; - psi_ref!, - PH_ref!, - region, - sweep_regions, - sweep_step, - solver_krylovdim=30, - solver_outputlevel=0, - solver_tol=1E-8, - substep, - time_step, - normalize, - ) - H = PH_ref![] - #applyexp tol is absolute, compute from tol_per_unit_time: - tol = abs(time_step) * tol_per_unit_time - psi, exp_info = applyexp( - H, time_step, init; tol, maxiter=solver_krylovdim, outputlevel=solver_outputlevel - ) - return psi, (; info=exp_info) - end - return solver -end From 5587391fdf72df8b8e0e0b4c88a05da265a0e9d0 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 14:12:09 -0500 Subject: [PATCH 17/32] Remove applyexp from ITensorNetworks.jl --- src/ITensorNetworks.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 8ef00052..a64a274c 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -108,8 +108,6 @@ include(joinpath("Graphs", "abstractgraph.jl")) include(joinpath("Graphs", "abstractdatagraph.jl")) include(joinpath("solvers", "eigsolve.jl")) include(joinpath("solvers", "exponentiate.jl")) -include(joinpath("treetensornetworks", "solvers", "applyexp.jl")) #this defines the primitive before the solver function -include(joinpath("solvers", "applyexp.jl")) include(joinpath("solvers", "dmrg_x_solver.jl")) include(joinpath("treetensornetworks", "abstracttreetensornetwork.jl")) include(joinpath("treetensornetworks", "ttn.jl")) From 06209435e988c8b4486cf7792180caf37470f15c Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 19 Jan 2024 15:54:53 -0500 Subject: [PATCH 18/32] Start renaming. One tdvp test not passing, observer related. --- src/solvers/dmrg_x_solver.jl | 6 +- src/solvers/eigsolve.jl | 6 +- src/solvers/exponentiate.jl | 10 +- .../solvers/alternating_update.jl | 52 ++++----- src/treetensornetworks/solvers/tdvp.jl | 26 ++--- src/treetensornetworks/solvers/update_step.jl | 96 ++++++++-------- .../test_solvers/test_tdvp.jl | 108 +++++++++--------- 7 files changed, 152 insertions(+), 152 deletions(-) diff --git a/src/solvers/dmrg_x_solver.jl b/src/solvers/dmrg_x_solver.jl index 6db73e9b..41d658f6 100644 --- a/src/solvers/dmrg_x_solver.jl +++ b/src/solvers/dmrg_x_solver.jl @@ -1,7 +1,7 @@ function dmrg_x_solver( init; - psi_ref!, - PH_ref!, + state!, + projected_operator!, normalize=nothing, region, sweep_regions, @@ -9,7 +9,7 @@ function dmrg_x_solver( half_sweep, step_kwargs..., ) - H = contract(PH_ref![], ITensor(1.0)) + H = contract(projected_operator![], ITensor(1.0)) D, U = eigen(H; ishermitian=true) u = uniqueind(U, H) max_overlap, max_ind = findmax(abs, array(dag(init) * U)) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index c0830100..8e16e760 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -1,7 +1,7 @@ function eigsolve_updater( init; - psi_ref!, - PH_ref!, + state!, + projected_operator!, outputlevel, which_sweep, region_updates, @@ -22,7 +22,7 @@ function eigsolve_updater( howmany = 1 which = updater_kwargs.which_eigenvalue vals, vecs, info = eigsolve( - PH_ref![], + projected_operator![], init, howmany, which; diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index 23609ae7..09231d58 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -1,7 +1,7 @@ function exponentiate_updater( init; - psi_ref!, - PH_ref!, + state!, + projected_operator!, outputlevel, which_sweep, region_updates, @@ -19,8 +19,8 @@ function exponentiate_updater( eager=true, ) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence - #H=copy(PH_ref![]) - H = PH_ref![] ###since we are not changing H we don't need the copy + #H=copy(projected_operator![]) + projected_operator = projected_operator![] ###since we are not changing H we don't need the copy # let's test whether given region and sweep regions we can find out what the previous and next region were # this will be needed in subspace expansion #@show step_kwargs @@ -33,7 +33,7 @@ function exponentiate_updater( previous_region = region_ind == 1 ? nothing : first(region_updates[region_ind - 1]) phi, exp_info = exponentiate( - H, + projected_operator, time_step, init; ishermitian=updater_kwargs.ishermitian, diff --git a/src/treetensornetworks/solvers/alternating_update.jl b/src/treetensornetworks/solvers/alternating_update.jl index 37f445f3..a2fe2698 100644 --- a/src/treetensornetworks/solvers/alternating_update.jl +++ b/src/treetensornetworks/solvers/alternating_update.jl @@ -26,10 +26,10 @@ function process_sweeps( return maxdim, mindim, cutoff, noise, kwargs end -function sweep_printer(; outputlevel, psi, which_sweep, sw_time) +function sweep_printer(; outputlevel, state, which_sweep, sw_time) if outputlevel >= 1 print("After sweep ", which_sweep, ":") - print(" maxlinkdim=", maxlinkdim(psi)) + print(" maxlinkdim=", maxlinkdim(state)) print(" cpu_time=", round(sw_time; digits=3)) println() flush(stdout) @@ -38,8 +38,8 @@ end function alternating_update( updater, - PH, - psi0::AbstractTTN; + projected_operator, + init_state::AbstractTTN; checkdone=(; kws...) -> false, outputlevel::Integer=0, nsweeps::Integer=1, @@ -51,7 +51,7 @@ function alternating_update( ) maxdim, mindim, cutoff, noise, kwargs = process_sweeps(nsweeps; kwargs...) - psi = copy(psi0) + state = copy(init_state) insert_function!(sweep_observer!, "sweep_printer" => sweep_printer) # FIX THIS @@ -63,7 +63,7 @@ function alternating_update( "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim[which_sweep] = $(maxdim[which_sweep]), writing environment tensors to disk", ) end - PH = disk(PH) + projected_operator = disk(projected_operator) end sweep_params=(; maxdim=maxdim[which_sweep], @@ -72,10 +72,10 @@ function alternating_update( noise=noise[which_sweep], ) sw_time = @elapsed begin - psi, PH = sweep_update( + state, projected_operator = sweep_update( updater, - PH, - psi; + projected_operator, + state; outputlevel, which_sweep, sweep_params, @@ -84,30 +84,30 @@ function alternating_update( ) end - update!(sweep_observer!; psi, which_sweep, sw_time, outputlevel) + update!(sweep_observer!; state, which_sweep, sw_time, outputlevel) - checkdone(; psi, which_sweep, outputlevel, kwargs...) && break + checkdone(; state, which_sweep, outputlevel, kwargs...) && break end select!(sweep_observer!, Observers.DataFrames.Not("sweep_printer")) - return psi + return state end -function alternating_update(updater, H::AbstractTTN, psi0::AbstractTTN; kwargs...) - check_hascommoninds(siteinds, H, psi0) - check_hascommoninds(siteinds, H, psi0') +function alternating_update(updater, H::AbstractTTN, init_state::AbstractTTN; kwargs...) + check_hascommoninds(siteinds, H, init_state) + check_hascommoninds(siteinds, H, init_state') # Permute the indices to have a better memory layout # and minimize permutations H = ITensors.permute(H, (linkind, siteinds, linkind)) - PH = ProjTTN(H) - return alternating_update(updater, PH, psi0; kwargs...) + projected_operator = ProjTTN(H) + return alternating_update(updater, projected_operator, init_state; kwargs...) end """ - tdvp(Hs::Vector{MPO},psi0::MPS,t::Number; kwargs...) - tdvp(Hs::Vector{MPO},psi0::MPS,t::Number, sweeps::Sweeps; kwargs...) + tdvp(Hs::Vector{MPO},init_state::MPS,t::Number; kwargs...) + tdvp(Hs::Vector{MPO},init_state::MPS,t::Number, sweeps::Sweeps; kwargs...) Use the time dependent variational principle (TDVP) algorithm -to compute `exp(t*H)*psi0` using an efficient algorithm based +to compute `exp(t*H)*init_state` using an efficient algorithm based on alternating optimization of the MPS tensors and local Krylov exponentiation of H. @@ -119,16 +119,16 @@ the set of MPOs [H1,H2,H3,..] is efficiently looped over at each step of the algorithm when optimizing the MPS. Returns: -* `psi::MPS` - time-evolved MPS +* `state::MPS` - time-evolved MPS """ function alternating_update( - updater, Hs::Vector{<:AbstractTTN}, psi0::AbstractTTN; kwargs... + updater, Hs::Vector{<:AbstractTTN}, init_state::AbstractTTN; kwargs... ) for H in Hs - check_hascommoninds(siteinds, H, psi0) - check_hascommoninds(siteinds, H, psi0') + check_hascommoninds(siteinds, H, init_state) + check_hascommoninds(siteinds, H, init_state') end Hs .= ITensors.permute.(Hs, Ref((linkind, siteinds, linkind))) - PHs = ProjTTNSum(Hs) - return alternating_update(updater, PHs, psi0; kwargs...) + projected_operators = ProjTTNSum(Hs) + return alternating_update(updater, projected_operators, init_state; kwargs...) end diff --git a/src/treetensornetworks/solvers/tdvp.jl b/src/treetensornetworks/solvers/tdvp.jl index 130c92e5..88ae3a27 100644 --- a/src/treetensornetworks/solvers/tdvp.jl +++ b/src/treetensornetworks/solvers/tdvp.jl @@ -58,21 +58,21 @@ end function tdvp( updater, - H, + operator, t::Number, - init::AbstractTTN; + init_state::AbstractTTN; time_step::Number=t, nsites=2, nsteps=nothing, order::Integer=2, (sweep_observer!)=observer(), - root_vertex=default_root_vertex(init), + root_vertex=default_root_vertex(init_state), reverse_step=true, updater_kwargs=(;), kwargs..., ) nsweeps = _compute_nsweeps(nsteps, t, time_step, order) - region_updates = tdvp_sweep(order, nsites, time_step, init; root_vertex, reverse_step) + region_updates = tdvp_sweep(order, nsites, time_step, init_state; root_vertex, reverse_step) function sweep_time_printer(; outputlevel, which_sweep, kwargs...) if outputlevel >= 1 @@ -87,26 +87,26 @@ function tdvp( insert_function!(sweep_observer!, "sweep_time_printer" => sweep_time_printer) - psi = alternating_update( - updater, H, init; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... + state = alternating_update( + updater, operator, init_state; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... ) # remove sweep_time_printer from sweep_observer! select!(sweep_observer!, Observers.DataFrames.Not("sweep_time_printer")) - return psi + return state end """ - tdvp(H::TTN, t::Number, psi0::TTN; kwargs...) + tdvp(operator::TTN, t::Number, init_state::TTN; kwargs...) Use the time dependent variational principle (TDVP) algorithm -to approximately compute `exp(H*t)*psi0` using an efficient algorithm based +to approximately compute `exp(operator*t)*init_state` using an efficient algorithm based on alternating optimization of the state tensors and local Krylov -exponentiation of H. The time parameter `t` can be a real or complex number. +exponentiation of operator. The time parameter `t` can be a real or complex number. Returns: -* `psi` - time-evolved state +* `state` - time-evolved state Optional keyword arguments: * `time_step::Number = t` - time step to use when evolving the state. Smaller time steps generally give more accurate results but can make the algorithm take more computational time to run. @@ -115,6 +115,6 @@ Optional keyword arguments: * `observer` - object implementing the Observer interface which can perform measurements and stop early * `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations """ -function tdvp(H, t::Number, init::AbstractTTN; updater=exponentiate_updater, kwargs...) - return tdvp(updater, H, t, init; kwargs...) +function tdvp(operator, t::Number, init_state::AbstractTTN; updater=exponentiate_updater, kwargs...) + return tdvp(updater, operator, t, init_state; kwargs...) end diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index c459e79f..cd48983c 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -15,7 +15,7 @@ function default_sweep_regions(nsites, graph::AbstractGraph; kwargs...) ###move end function region_update_printer(; - cutoff, maxdim, mindim, outputlevel::Int=0, psi, region_updates, spec, which_region_update, which_sweep,kwargs... + cutoff, maxdim, mindim, outputlevel::Int=0, state, region_updates, spec, which_region_update, which_sweep,kwargs... ) if outputlevel >= 2 region=first(region_updates[which_region_update]) @@ -29,7 +29,7 @@ function region_update_printer(; @printf( " Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, - linkdim(psi, edgetype(psi)(region...)) + linkdim(state, edgetype(state)(region...)) ) end flush(stdout) @@ -38,15 +38,15 @@ end function sweep_update( solver, - PH, - psi::AbstractTTN; + projected_operator, + state::AbstractTTN; normalize::Bool=false, # ToDo: think about where to put the default, probably this default is best defined at algorithmic level outputlevel, region_update_printer=region_update_printer, (region_observer!)=observer(), # ToDo: change name to region_observer! ? which_sweep::Int, sweep_params::NamedTuple, - region_updates,# =default_sweep_regions(nsite, psi; reverse_step), #move default up to algorithmic level + region_updates,# =default_sweep_regions(nsite, state; reverse_step), #move default up to algorithmic level updater_kwargs, ) insert_function!(region_observer!, "region_update_printer" => region_update_printer) #ToDo fix this @@ -55,7 +55,7 @@ function sweep_update( # (Needed to handle user-provided region_updates) region_updates = append_missing_namedtuple.(to_tuple.(region_updates)) - if nv(psi) == 1 + if nv(state) == 1 error( "`alternating_update` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.exponentiate`, etc.", ) @@ -65,10 +65,10 @@ function sweep_update( # merge sweep params in step_kwargs (region, region_kwargs)=region_updates[which_region_update] region_kwargs=merge(region_kwargs, sweep_params) #in this case sweep params has precedence over step_kwargs --- correct behaviour? - psi, PH = region_update( + state, projected_operator = region_update( solver, - PH, - psi; + projected_operator, + state; normalize, outputlevel, which_sweep, @@ -82,9 +82,9 @@ function sweep_update( select!(region_observer!, Observers.DataFrames.Not("region_update_printer")) # remove update_printer # Just to be sure: - normalize && normalize!(psi) + normalize && normalize!(state) - return psi, PH + return state, projected_operator end # @@ -98,20 +98,20 @@ end # apart and puts it back into the network. # -function extract_local_tensor(psi::AbstractTTN, pos::Vector) - return psi, prod(psi[v] for v in pos) +function extract_local_tensor(state::AbstractTTN, pos::Vector) + return state, prod(state[v] for v in pos) end -function extract_local_tensor(psi::AbstractTTN, e::NamedEdge) - left_inds = uniqueinds(psi, e) - U, S, V = svd(psi[src(e)], left_inds; lefttags=tags(psi, e), righttags=tags(psi, e)) - psi[src(e)] = U - return psi, S * V +function extract_local_tensor(state::AbstractTTN, e::NamedEdge) + left_inds = uniqueinds(state, e) + U, S, V = svd(state[src(e)], left_inds; lefttags=tags(state, e), righttags=tags(state, e)) + state[src(e)] = U + return state, S * V end # sort of multi-site replacebond!; TODO: use dense TTN constructor instead function insert_local_tensor( - psi::AbstractTTN, + state::AbstractTTN, phi::ITensor, pos::Vector; normalize=false, @@ -125,12 +125,12 @@ function insert_local_tensor( ) spec = nothing for (v, vnext) in IterTools.partition(pos, 2, 1) - e = edgetype(psi)(v, vnext) - indsTe = inds(psi[v]) + e = edgetype(state)(v, vnext) + indsTe = inds(state[v]) L, phi, spec = factorize( phi, indsTe; - tags=tags(psi, e), + tags=tags(state, e), maxdim, mindim, cutoff, @@ -138,21 +138,21 @@ function insert_local_tensor( eigen_perturbation, ortho, ) - psi[v] = L + state[v] = L eigen_perturbation = nothing # TODO: fix this end - psi[last(pos)] = phi - psi = set_ortho_center(psi, [last(pos)]) - @assert isortho(psi) && only(ortho_center(psi)) == last(pos) - normalize && (psi[last(pos)] ./= norm(psi[last(pos)])) + state[last(pos)] = phi + state = set_ortho_center(state, [last(pos)]) + @assert isortho(state) && only(ortho_center(state)) == last(pos) + normalize && (state[last(pos)] ./= norm(state[last(pos)])) # TODO: return maxtruncerr, will not be correct in cases where insertion executes multiple factorizations - return psi, spec + return state, spec end -function insert_local_tensor(psi::AbstractTTN, phi::ITensor, e::NamedEdge; kwargs...) - psi[dst(e)] *= phi - psi = set_ortho_center(psi, [dst(e)]) - return psi, nothing +function insert_local_tensor(state::AbstractTTN, phi::ITensor, e::NamedEdge; kwargs...) + state[dst(e)] *= phi + state = set_ortho_center(state, [dst(e)]) + return state, nothing end #TODO: clean this up: @@ -163,8 +163,8 @@ current_ortho(st) = current_ortho(typeof(st), st) function region_update( updater, - PH, - psi; + projected_operator, + state; normalize, outputlevel, which_sweep, @@ -177,17 +177,17 @@ function region_update( updater_kwargs ) region=first(region_updates[which_region_update]) - psi = orthogonalize(psi, current_ortho(region)) - psi, phi = extract_local_tensor(psi, region;) + state = orthogonalize(state, current_ortho(region)) + state, phi = extract_local_tensor(state, region;) nsites = (region isa AbstractEdge) ? 0 : length(region) #ToDo move into separate funtion - PH = set_nsite(PH, nsites) - PH = position(PH, psi, region) - psi_ref! = Ref(psi) # create references, in case solver does (out-of-place) modify PH or psi - PH_ref! = Ref(PH) + projected_operator = set_nsite(projected_operator, nsites) + projected_operator = position(projected_operator, state, region) + state! = Ref(state) # create references, in case solver does (out-of-place) modify PH or state + projected_operator! = Ref(projected_operator) phi, info = updater( phi; - psi_ref!, - PH_ref!, + state!, + projected_operator!, outputlevel, which_sweep, region_updates, @@ -195,8 +195,8 @@ function region_update( region_kwargs, updater_kwargs ) # args passed by reference are supposed to be modified out of place - psi = psi_ref![] # dereference - PH = PH_ref![] + state = state![] # dereference + projected_operator = projected_operator![] if !(phi isa ITensor && info isa NamedTuple) println("Solver returned the following types: $(typeof(phi)), $(typeof(info))") error("In alternating_update, solver must return an ITensor and a NamedTuple") @@ -210,8 +210,8 @@ function region_update( # so noiseterm is a solver #end - psi, spec = insert_local_tensor( - psi, phi, region; eigen_perturbation=drho, ortho, normalize, + state, spec = insert_local_tensor( + state, phi, region; eigen_perturbation=drho, ortho, normalize, maxdim=region_kwargs.maxdim, mindim=region_kwargs.mindim, cutoff=region_kwargs.cutoff @@ -226,7 +226,7 @@ function region_update( region_updates, total_sweep_steps=length(region_updates), end_of_sweep=(which_region_update == length(region_updates)), - psi, + state, region, which_sweep, spec, @@ -234,5 +234,5 @@ function region_update( info..., region_kwargs..., ) - return psi, PH + return state, projected_operator end diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index 496d4b21..68d9a49b 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -164,15 +164,15 @@ using Test ψ0 = random_mps(s; internal_inds_space=10) - function updater(psi0; (psi_ref!), (PH_ref!), kwargs...) + function updater(psi0; (state!), (projected_operator!), kwargs...) updater_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true ) time_step=get(step_kwargs,:time_step,nothing) @assert !isnothing(time_step) - PH = PH_ref![] - psi, exp_info = exponentiate(PH, time_step, psi0; updater_kwargs...) - return psi, (; info=exp_info) + PH = projected_operator![] + state, exp_info = exponentiate(PH, time_step, psi0; updater_kwargs...) + return state, (; info=exp_info) end ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsites=1) @@ -213,9 +213,9 @@ using Test Ut = exp(-im * tau * HM) - psi = mps(s; states=(n -> isodd(n) ? "Up" : "Dn")) - psi2 = deepcopy(psi) - psix = contract(psi) + state = mps(s; states=(n -> isodd(n) ? "Up" : "Dn")) + psi2 = deepcopy(state) + psix = contract(state) Sz_tdvp = Float64[] Sz_tdvp2 = Float64[] @@ -229,10 +229,10 @@ using Test psix = noprime(Ut * psix) psix /= norm(psix) - psi = tdvp( + state = tdvp( H, -im * tau, - psi; + state; cutoff, normalize=false, updater_kwargs=(; @@ -242,7 +242,7 @@ using Test ) # TODO: What should `expect` output? Right now # it outputs a dictionary. - push!(Sz_tdvp, real(expect("Sz", psi; vertices=[c])[c])) + push!(Sz_tdvp, real(expect("Sz", state; vertices=[c])[c])) psi2 = tdvp( H, @@ -261,7 +261,7 @@ using Test push!(Sz_tdvp2, real(expect("Sz", psi2; vertices=[c])[c])) push!(Sz_exact, real(scalar(dag(prime(psix, s[c])) * Szc * psix))) - F = abs(scalar(dag(psix) * contract(psi))) + F = abs(scalar(dag(psix) * contract(state))) end @test norm(Sz_tdvp - Sz_exact) < 1e-5 @@ -298,8 +298,8 @@ using Test end append!(gates, reverse(gates)) - psi = mps(s; states=(n -> isodd(n) ? "Up" : "Dn")) - phi = deepcopy(psi) + state = mps(s; states=(n -> isodd(n) ? "Up" : "Dn")) + phi = deepcopy(state) c = div(N, 2) # @@ -313,17 +313,17 @@ using Test #En2 = zeros(Nsteps) for step in 1:Nsteps - psi = apply(gates, psi; cutoff) - #normalize!(psi) + state = apply(gates, state; cutoff) + #normalize!(state) nsites = (step <= 3 ? 2 : 1) phi = tdvp( H, -tau * im, phi; nsteps=1, cutoff, nsites, normalize=true, updater_kwargs=(;krylovdim=15) ) - Sz1[step] = real(expect("Sz", psi; vertices=[c])[c]) + Sz1[step] = real(expect("Sz", state; vertices=[c])[c]) #Sz2[step] = real(expect("Sz", phi; vertices=[c])[c]) - En1[step] = real(inner(psi', H, psi)) + En1[step] = real(inner(state', H, state)) #En2[step] = real(inner(phi', H, phi)) end @@ -334,8 +334,8 @@ using Test phi = mps(s; states=(n -> isodd(n) ? "Up" : "Dn")) obs = Observer( - "Sz" => (; psi) -> expect("Sz", psi; vertices=[c])[c], - "En" => (; psi) -> real(inner(psi', H, psi)), + "Sz" => (; state) -> expect("Sz", state; vertices=[c])[c], + "En" => (; state) -> real(inner(state', H, state)), ) phi = tdvp( @@ -372,13 +372,13 @@ using Test H = mpo(os, s) - psi = random_mps(s; internal_inds_space=2) - psi2 = deepcopy(psi) + state = random_mps(s; internal_inds_space=2) + psi2 = deepcopy(state) trange = 0.0:tau:ttotal for (step, t) in enumerate(trange) nsites = (step <= 10 ? 2 : 1) - psi = tdvp( - H, -tau, psi; cutoff, nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) + state = tdvp( + H, -tau, state; cutoff, nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) ) #Different backend updaters, default updater_backend = "applyexp" psi2 = tdvp( @@ -394,9 +394,9 @@ using Test ) end - @test psi ≈ psi2 rtol = 1e-6 + @test state ≈ psi2 rtol = 1e-6 - en1 = inner(psi', H, psi) + en1 = inner(state', H, state) en2 = inner(psi2', H, psi2) @test en1 < -4.25 @test en1 ≈ en2 @@ -424,22 +424,22 @@ using Test # Using Observers.jl # - measure_sz(; psi) = expect("Sz", psi; vertices=[c])[c] - measure_en(; psi) = real(inner(psi', H, psi)) + measure_sz(; state) = expect("Sz", state; vertices=[c])[c] + measure_en(; state) = real(inner(state', H, state)) sweep_obs = Observer("Sz" => measure_sz, "En" => measure_en) get_info(; info) = info - step_measure_sz(; psi) = expect("Sz", psi; vertices=[c])[c] - step_measure_en(; psi) = real(inner(psi', H, psi)) + step_measure_sz(; state) = expect("Sz", state; vertices=[c])[c] + step_measure_en(; state) = real(inner(state', H, state)) region_obs = Observer( "Sz" => step_measure_sz, "En" => step_measure_en, "info" => get_info ) - psi2 = mps(s; states=(n -> isodd(n) ? "Up" : "Dn")) + state2 = mps(s; states=(n -> isodd(n) ? "Up" : "Dn")) tdvp( H, -im * ttotal, - psi2; + state2; time_step=-im * tau, cutoff, normalize=false, @@ -557,11 +557,11 @@ end ψ0 = normalize!(random_ttn(s; link_space=10)) - function updater(psi0; (psi_ref!), (PH_ref!), time_step, kwargs...) + function updater(psi0; (state!), (projected_operator!), time_step, kwargs...) updater_kwargs = (; ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true ) - PH = PH_ref![] + PH = projected_operator![] psi, exp_info = exponentiate(PH, time_step, psi0; updater_kwargs...) return psi, (; info=exp_info) end @@ -604,8 +604,8 @@ end Ut = exp(-im * tau * HM) - psi = TTN(ComplexF64, s, v -> iseven(sum(isodd.(v))) ? "Up" : "Dn") - psix = contract(psi) + state = TTN(ComplexF64, s, v -> iseven(sum(isodd.(v))) ? "Up" : "Dn") + statex = contract(state) Sz_tdvp = Float64[] Sz_exact = Float64[] @@ -615,13 +615,13 @@ end Nsteps = Int(ttotal / tau) for step in 1:Nsteps - psix = noprime(Ut * psix) - psix /= norm(psix) + statex = noprime(Ut * statex) + statex /= norm(statex) - psi = tdvp( + state = tdvp( H, -im * tau, - psi; + state; cutoff, normalize=false, updater_kwargs=(; @@ -629,9 +629,9 @@ end maxiter=500, krylovdim=25,) ) - push!(Sz_tdvp, real(expect("Sz", psi; vertices=[c])[c])) - push!(Sz_exact, real(scalar(dag(prime(psix, s[c])) * Szc * psix))) - F = abs(scalar(dag(psix) * contract(psi))) + push!(Sz_tdvp, real(expect("Sz", state; vertices=[c])[c])) + push!(Sz_exact, real(scalar(dag(prime(statex, s[c])) * Szc * statex))) + F = abs(scalar(dag(statex) * contract(state))) end @test norm(Sz_tdvp - Sz_exact) < 1e-5 @@ -665,8 +665,8 @@ end end append!(gates, reverse(gates)) - psi = TTN(s, v -> iseven(sum(isodd.(v))) ? "Up" : "Dn") - phi = copy(psi) + state = TTN(s, v -> iseven(sum(isodd.(v))) ? "Up" : "Dn") + phi = copy(state) c = (2, 1) # @@ -680,17 +680,17 @@ end En2 = zeros(Nsteps) for step in 1:Nsteps - psi = apply(gates, psi; cutoff, maxdim) - #normalize!(psi) + state = apply(gates, state; cutoff, maxdim) + #normalize!(state) nsites = (step <= 3 ? 2 : 1) phi = tdvp( H, -tau * im, phi; nsteps=1, cutoff, nsites, normalize=true, updater_kwargs=(;krylovdim=15) ) - Sz1[step] = real(expect("Sz", psi; vertices=[c])[c]) + Sz1[step] = real(expect("Sz", state; vertices=[c])[c]) Sz2[step] = real(expect("Sz", phi; vertices=[c])[c]) - En1[step] = real(inner(psi', H, psi)) + En1[step] = real(inner(state', H, state)) En2[step] = real(inner(phi', H, phi)) end @@ -700,8 +700,8 @@ end phi = TTN(s, v -> iseven(sum(isodd.(v))) ? "Up" : "Dn") obs = Observer( - "Sz" => (; psi) -> expect("Sz", psi; vertices=[c])[c], - "En" => (; psi) -> real(inner(psi', H, psi)), + "Sz" => (; state) -> expect("Sz", state; vertices=[c])[c], + "En" => (; state) -> real(inner(state', H, state)), ) phi = tdvp( H, @@ -730,17 +730,17 @@ end os = ITensorNetworks.heisenberg(c) H = TTN(os, s) - psi = normalize!(random_ttn(s; link_space=2)) + state = normalize!(random_ttn(s; link_space=2)) trange = 0.0:tau:ttotal for (step, t) in enumerate(trange) nsites = (step <= 10 ? 2 : 1) - psi = tdvp( - H, -tau, psi; cutoff, nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) + state = tdvp( + H, -tau, state; cutoff, nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) ) end - @test inner(psi', H, psi) < -2.47 + @test inner(state', H, state) < -2.47 end # TODO: verify quantum number suport in ITensorNetworks From 7865f818a751cd828d191bcf5085690aaa2c7c94 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 16:41:51 -0500 Subject: [PATCH 19/32] Cleanup solvers. --- src/solvers/eigsolve.jl | 9 ++------- src/solvers/exponentiate.jl | 29 +++++------------------------ 2 files changed, 7 insertions(+), 31 deletions(-) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 8e16e760..98f9c28f 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -15,7 +15,7 @@ function eigsolve_updater( tol=1e-14, krylovdim=3, maxiter=1, - outputlevel=0, + verbosity=0, eager=false, ) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence @@ -26,12 +26,7 @@ function eigsolve_updater( init, howmany, which; - ishermitian=updater_kwargs.ishermitian, - tol=updater_kwargs.tol, - krylovdim=updater_kwargs.krylovdim, - maxiter=updater_kwargs.maxiter, - verbosity=updater_kwargs.outputlevel, - eager=updater_kwargs.eager, + updater_kwargs... # leaves it to the user to supply only supported kwargs ) return vecs[1], (; info, eigvals=vals) end diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index 09231d58..8cd0e3ff 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -12,37 +12,18 @@ function exponentiate_updater( default_updater_kwargs = (; krylovdim=30, #from here only solver kwargs maxiter=100, - outputlevel=0, + verbosity=0, tol=1E-12, ishermitian=true, issymmetric=true, eager=true, ) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence - #H=copy(projected_operator![]) - projected_operator = projected_operator![] ###since we are not changing H we don't need the copy - # let's test whether given region and sweep regions we can find out what the previous and next region were - # this will be needed in subspace expansion - #@show step_kwargs - substep = get(region_kwargs, :substep, nothing) - time_step = get(region_kwargs, :time_step, nothing) - @assert !isnothing(time_step) && !isnothing(substep) - region_ind = which_region_update - next_region = - region_ind == length(region_updates) ? nothing : first(region_updates[region_ind + 1]) - previous_region = region_ind == 1 ? nothing : first(region_updates[region_ind - 1]) - - phi, exp_info = exponentiate( - projected_operator, + result, exp_info = exponentiate( + projected_operator![], time_step, init; - ishermitian=updater_kwargs.ishermitian, - issymmetric=updater_kwargs.issymmetric, - tol=updater_kwargs.tol, - krylovdim=updater_kwargs.krylovdim, - maxiter=updater_kwargs.maxiter, - verbosity=updater_kwargs.outputlevel, - eager=updater_kwargs.eager, + updater_kwargs... ) - return phi, (; info=exp_info) + return result, (; info=exp_info) end From e08ba96cba8309f9a75b9505fc375a93542759df Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 16:58:07 -0500 Subject: [PATCH 20/32] Rename to sweep_plan. --- src/treetensornetworks/solvers/dmrg.jl | 8 ++++---- src/treetensornetworks/solvers/tdvp.jl | 12 ++++++------ src/treetensornetworks/solvers/update_step.jl | 9 ++++----- 3 files changed, 14 insertions(+), 15 deletions(-) diff --git a/src/treetensornetworks/solvers/dmrg.jl b/src/treetensornetworks/solvers/dmrg.jl index b5525277..653c00c8 100644 --- a/src/treetensornetworks/solvers/dmrg.jl +++ b/src/treetensornetworks/solvers/dmrg.jl @@ -2,12 +2,12 @@ Overload of `ITensors.dmrg`. """ -function dmrg_sweep( +function dmrg_sweep_plan( nsites::Int, graph::AbstractGraph; root_vertex=default_root_vertex(graph) ) order = 2 time_step = Inf - return tdvp_sweep(order, nsites, time_step, graph; root_vertex, reverse_step=false) + return tdvp_sweep_plan(order, nsites, time_step, graph; root_vertex, reverse_step=false) end function dmrg( @@ -21,10 +21,10 @@ function dmrg( updater_kwargs=(;), kwargs..., ) - region_updates = dmrg_sweep(nsites, init; root_vertex) + sweep_plan = dmrg_sweep_plan(nsites, init; root_vertex) psi = alternating_update( - updater, H, init; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... + updater, H, init; nsweeps, sweep_observer!, sweep_plan, updater_kwargs, kwargs... ) return psi end diff --git a/src/treetensornetworks/solvers/tdvp.jl b/src/treetensornetworks/solvers/tdvp.jl index 88ae3a27..98f21849 100644 --- a/src/treetensornetworks/solvers/tdvp.jl +++ b/src/treetensornetworks/solvers/tdvp.jl @@ -30,7 +30,7 @@ function sub_time_steps(order) end end -function tdvp_sweep( +function tdvp_sweep_plan( order::Int, nsites::Int, time_step::Number, @@ -38,7 +38,7 @@ function tdvp_sweep( root_vertex=default_root_vertex(graph), reverse_step=true, ) - sweep = [] + sweep_plan = [] for (substep, fac) in enumerate(sub_time_steps(order)) sub_time_step = time_step * fac half = half_sweep( @@ -51,9 +51,9 @@ function tdvp_sweep( reverse_args=(; substep, time_step=-sub_time_step), reverse_step, ) - append!(sweep, half) + append!(sweep_plan, half) end - return sweep + return sweep_plan end function tdvp( @@ -72,7 +72,7 @@ function tdvp( kwargs..., ) nsweeps = _compute_nsweeps(nsteps, t, time_step, order) - region_updates = tdvp_sweep(order, nsites, time_step, init_state; root_vertex, reverse_step) + sweep_plan = tdvp_sweep_plan(order, nsites, time_step, init_state; root_vertex, reverse_step) function sweep_time_printer(; outputlevel, which_sweep, kwargs...) if outputlevel >= 1 @@ -88,7 +88,7 @@ function tdvp( insert_function!(sweep_observer!, "sweep_time_printer" => sweep_time_printer) state = alternating_update( - updater, operator, init_state; nsweeps, sweep_observer!, region_updates, updater_kwargs, kwargs... + updater, operator, init_state; nsweeps, sweep_observer!, sweep_plan, updater_kwargs, kwargs... ) # remove sweep_time_printer from sweep_observer! diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index cd48983c..90452c43 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -46,7 +46,7 @@ function sweep_update( (region_observer!)=observer(), # ToDo: change name to region_observer! ? which_sweep::Int, sweep_params::NamedTuple, - region_updates,# =default_sweep_regions(nsite, state; reverse_step), #move default up to algorithmic level + sweep_plan, updater_kwargs, ) insert_function!(region_observer!, "region_update_printer" => region_update_printer) #ToDo fix this @@ -61,10 +61,9 @@ function sweep_update( ) end - for which_region_update in eachindex(region_updates) - # merge sweep params in step_kwargs - (region, region_kwargs)=region_updates[which_region_update] - region_kwargs=merge(region_kwargs, sweep_params) #in this case sweep params has precedence over step_kwargs --- correct behaviour? + for which_region_update in eachindex(sweep_plan) + (region, region_kwargs)=sweep_plan[which_region_update] + region_kwargs=merge(region_kwargs, sweep_params) # sweep params has precedence over step_kwargs state, projected_operator = region_update( solver, projected_operator, From 30455a4c09ce4156bd2fc344e4b4d222aa6af206 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 17:32:15 -0500 Subject: [PATCH 21/32] Fix renaming to tdvp_sweep and kwarg handling in updaters. --- src/solvers/eigsolve.jl | 10 +++++++--- src/solvers/exponentiate.jl | 7 ++++--- src/treetensornetworks/solvers/update_step.jl | 20 +++++++++---------- .../test_solvers/test_dmrg.jl | 4 ++-- .../test_solvers/test_tdvp.jl | 3 ++- 5 files changed, 25 insertions(+), 19 deletions(-) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index 98f9c28f..eb9e8aef 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -4,7 +4,7 @@ function eigsolve_updater( projected_operator!, outputlevel, which_sweep, - region_updates, + sweep_plan, which_region_update, region_kwargs, updater_kwargs, @@ -20,13 +20,17 @@ function eigsolve_updater( ) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence howmany = 1 - which = updater_kwargs.which_eigenvalue + which, updater_kwargs = _pop_which_eigenvalue(;updater_kwargs...) vals, vecs, info = eigsolve( projected_operator![], init, howmany, which; - updater_kwargs... # leaves it to the user to supply only supported kwargs + updater_kwargs... #this leaves it ) return vecs[1], (; info, eigvals=vals) end + +function _pop_which_eigenvalue(;which_eigenvalue, kwargs...) + return which_eigenvalue, NamedTuple(kwargs) +end \ No newline at end of file diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index 8cd0e3ff..7d56ad39 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -4,13 +4,13 @@ function exponentiate_updater( projected_operator!, outputlevel, which_sweep, - region_updates, + sweep_plan, which_region_update, region_kwargs, updater_kwargs, ) default_updater_kwargs = (; - krylovdim=30, #from here only solver kwargs + krylovdim=30, maxiter=100, verbosity=0, tol=1E-12, @@ -18,10 +18,11 @@ function exponentiate_updater( issymmetric=true, eager=true, ) + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence result, exp_info = exponentiate( projected_operator![], - time_step, + region_kwargs.time_step, init; updater_kwargs... ) diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index 90452c43..445fb329 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -15,10 +15,10 @@ function default_sweep_regions(nsites, graph::AbstractGraph; kwargs...) ###move end function region_update_printer(; - cutoff, maxdim, mindim, outputlevel::Int=0, state, region_updates, spec, which_region_update, which_sweep,kwargs... + cutoff, maxdim, mindim, outputlevel::Int=0, state, sweep_plan, spec, which_region_update, which_sweep,kwargs... ) if outputlevel >= 2 - region=first(region_updates[which_region_update]) + region=first(sweep_plan[which_region_update]) @printf("Sweep %d, region=%s \n", which_sweep, region) print(" Truncated using") @printf(" cutoff=%.1E", cutoff) @@ -53,7 +53,7 @@ function sweep_update( # Append empty namedtuple to each element if not already present # (Needed to handle user-provided region_updates) - region_updates = append_missing_namedtuple.(to_tuple.(region_updates)) + sweep_plan = append_missing_namedtuple.(to_tuple.(sweep_plan)) if nv(state) == 1 error( @@ -71,7 +71,7 @@ function sweep_update( normalize, outputlevel, which_sweep, - region_updates, + sweep_plan, which_region_update, region_kwargs, region_observer!, @@ -167,7 +167,7 @@ function region_update( normalize, outputlevel, which_sweep, - region_updates, + sweep_plan, which_region_update, region_kwargs, region_observer!, @@ -175,7 +175,7 @@ function region_update( #extraction_kwargs, #ToDo: implement later with possibility to pass custom extraction/insertion func (or code into func) updater_kwargs ) - region=first(region_updates[which_region_update]) + region=first(sweep_plan[which_region_update]) state = orthogonalize(state, current_ortho(region)) state, phi = extract_local_tensor(state, region;) nsites = (region isa AbstractEdge) ? 0 : length(region) #ToDo move into separate funtion @@ -189,7 +189,7 @@ function region_update( projected_operator!, outputlevel, which_sweep, - region_updates, + sweep_plan, which_region_update, region_kwargs, updater_kwargs @@ -222,9 +222,9 @@ function region_update( maxdim, mindim, which_region_update, - region_updates, - total_sweep_steps=length(region_updates), - end_of_sweep=(which_region_update == length(region_updates)), + sweep_plan, + total_sweep_steps=length(sweep_plan), + end_of_sweep=(which_region_update == length(sweep_plan)), state, region, which_sweep, diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl index 92a318d7..9bc3f608 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl @@ -74,8 +74,8 @@ end sweep(; which_sweep, kw...) = which_sweep sweep_observer! = observer(sweep) - region(; which_region_update, region_updates, kw...) = - first(region_updates[which_region_update]) + region(; which_region_update, sweep_plan, kw...) = + first(sweep_plan[which_region_update]) energy(; eigvals, kw...) = eigvals[1] region_observer! = observer(region, sweep, energy) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index 68d9a49b..d53c081c 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -542,7 +542,8 @@ end # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end -#= + # ToDo: Discuss whether the test commented out here is necessary given the new design? + #= @testset "Custom updater in TDVP" begin cutoff = 1e-12 From 86dd746f9c9742fa380f75eb346ba2dcbe08c956 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 17:50:39 -0500 Subject: [PATCH 22/32] Fix dmrg-x. --- src/ITensorNetworks.jl | 2 +- src/solvers/{dmrg_x_solver.jl => dmrg_x.jl} | 20 +++++++++------- src/treetensornetworks/solvers/dmrg_x.jl | 23 +++++++++++++++++-- .../test_solvers/test_dmrg_x.jl | 19 +++++++-------- 4 files changed, 42 insertions(+), 22 deletions(-) rename src/solvers/{dmrg_x_solver.jl => dmrg_x.jl} (51%) diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index a64a274c..af0b4992 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -108,7 +108,7 @@ include(joinpath("Graphs", "abstractgraph.jl")) include(joinpath("Graphs", "abstractdatagraph.jl")) include(joinpath("solvers", "eigsolve.jl")) include(joinpath("solvers", "exponentiate.jl")) -include(joinpath("solvers", "dmrg_x_solver.jl")) +include(joinpath("solvers", "dmrg_x.jl")) include(joinpath("treetensornetworks", "abstracttreetensornetwork.jl")) include(joinpath("treetensornetworks", "ttn.jl")) include(joinpath("treetensornetworks", "opsum_to_ttn.jl")) diff --git a/src/solvers/dmrg_x_solver.jl b/src/solvers/dmrg_x.jl similarity index 51% rename from src/solvers/dmrg_x_solver.jl rename to src/solvers/dmrg_x.jl index 41d658f6..43ba69de 100644 --- a/src/solvers/dmrg_x_solver.jl +++ b/src/solvers/dmrg_x.jl @@ -1,19 +1,23 @@ -function dmrg_x_solver( +function dmrg_x_updater( init; state!, projected_operator!, - normalize=nothing, - region, - sweep_regions, - sweep_step, - half_sweep, - step_kwargs..., + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs, ) + # this updater does not seem to accept any kwargs? + default_updater_kwargs = (; + ) + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) H = contract(projected_operator![], ITensor(1.0)) D, U = eigen(H; ishermitian=true) u = uniqueind(U, H) max_overlap, max_ind = findmax(abs, array(dag(init) * U)) U_max = U * dag(onehot(u => max_ind)) # TODO: improve this to return the energy estimate too - return U_max, NamedTuple() + return U_max, (;) end diff --git a/src/treetensornetworks/solvers/dmrg_x.jl b/src/treetensornetworks/solvers/dmrg_x.jl index 67ab42c9..b5f8405b 100644 --- a/src/treetensornetworks/solvers/dmrg_x.jl +++ b/src/treetensornetworks/solvers/dmrg_x.jl @@ -1,4 +1,23 @@ -function dmrg_x(PH, init::AbstractTTN; kwargs...) - psi = alternating_update(ITensorNetworks.dmrg_x_solver, PH, init; kwargs...) +function dmrg_x( + updater, + operator, + init::AbstractTTN; + nsweeps, #it makes sense to require this to be defined + nsites=2, + (sweep_observer!)=observer(), + root_vertex=default_root_vertex(init), + updater_kwargs=(;), + kwargs..., +) + sweep_plan = dmrg_sweep_plan(nsites, init; root_vertex) + + psi = alternating_update( + updater, operator, init; nsweeps, sweep_observer!, sweep_plan, updater_kwargs, kwargs... + ) return psi end + +function dmrg_x(operator, init::AbstractTTN; updater=dmrg_x_updater, kwargs...) + return dmrg_x(updater, operator, init; kwargs...) +end + diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl index 65871b56..1353b8b8 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl @@ -17,16 +17,16 @@ using Test ψ = mps(s; states=(v -> rand(["↑", "↓"]))) dmrg_x_kwargs = ( - nsweeps=20, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0 + nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0 ) - ϕ = dmrg_x(H, ψ; nsite=2, dmrg_x_kwargs...) + ϕ = dmrg_x(H, ψ; nsites=2, dmrg_x_kwargs...) @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1 @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = 1e-7 @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = 1e-7 - ϕ̃ = dmrg_x(H, ϕ; nsite=1, dmrg_x_kwargs...) + ϕ̃ = dmrg_x(H, ϕ; nsites=1, dmrg_x_kwargs...) @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) rtol = 1e-1 @test inner(H, ϕ̃, H, ϕ̃) ≈ inner(ϕ̃', H, ϕ̃)^2 rtol = 1e-3 @@ -34,10 +34,7 @@ using Test # @test abs(loginner(ϕ̃, ϕ) / n) ≈ 0.0 atol = 1e-6 end -@testset "Tree DMRG-X" for conserve_qns in ( - false, - true, # OpSum → TTN with QNs not working for non-path graphs -) +@testset "Tree DMRG-X" for conserve_qns in (false, true) tooth_lengths = fill(2, 3) root_vertex = (3, 2) c = named_comb_tree(tooth_lengths) @@ -54,18 +51,18 @@ end # TODO: Use `TTN(s; states=v -> rand(["↑", "↓"]))` or # `ttns(s; states=v -> rand(["↑", "↓"]))` ψ = normalize!(TTN(s, v -> rand(["↑", "↓"]))) - + dmrg_x_kwargs = ( - nsweeps=20, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0 + nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0 ) - ϕ = dmrg_x(H, ψ; nsite=2, dmrg_x_kwargs...) + ϕ = dmrg_x(H, ψ; nsites=2, dmrg_x_kwargs...) @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1 @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = 1e-2 @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = 1e-7 - ϕ̃ = dmrg_x(H, ϕ; nsite=1, dmrg_x_kwargs...) + ϕ̃ = dmrg_x(H, ϕ; nsites=1, dmrg_x_kwargs...) @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) rtol = 1e-1 @test inner(H, ϕ̃, H, ϕ̃) ≈ inner(ϕ̃', H, ϕ̃)^2 rtol = 1e-6 From 59f45c17c84f8f093dadd9dd9d90e78e1d5ed13f Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 18:20:49 -0500 Subject: [PATCH 23/32] Adapt contract(_updater) to new interface. --- src/ITensorNetworks.jl | 1 + src/solvers/contract.jl | 19 +++++++++++++++++++ src/treetensornetworks/solvers/contract.jl | 14 ++++---------- .../test_solvers/test_contract.jl | 4 ++-- 4 files changed, 26 insertions(+), 12 deletions(-) create mode 100644 src/solvers/contract.jl diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index af0b4992..2f3b4e03 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -109,6 +109,7 @@ include(joinpath("Graphs", "abstractdatagraph.jl")) include(joinpath("solvers", "eigsolve.jl")) include(joinpath("solvers", "exponentiate.jl")) include(joinpath("solvers", "dmrg_x.jl")) +include(joinpath("solvers", "contract.jl")) include(joinpath("treetensornetworks", "abstracttreetensornetwork.jl")) include(joinpath("treetensornetworks", "ttn.jl")) include(joinpath("treetensornetworks", "opsum_to_ttn.jl")) diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl new file mode 100644 index 00000000..b3e3223d --- /dev/null +++ b/src/solvers/contract.jl @@ -0,0 +1,19 @@ +function contract_updater( + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs, + ) + v = ITensor(1.0) + projected_operator = projected_operator![] + for j in sites(projected_operator) + v *= projected_operator.psi0[j] + end + Hpsi0 = contract(projected_operator, v) + return Hpsi0, (;) + end \ No newline at end of file diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 51f6787f..92e4f744 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -1,18 +1,11 @@ -function contract_solver(PH, psi; normalize, region, half_sweep) - v = ITensor(1.0) - for j in sites(PH) - v *= PH.psi0[j] - end - Hpsi0 = contract(PH, v) - return Hpsi0, NamedTuple() -end - function contract( ::Algorithm"fit", tn1::AbstractTTN, tn2::AbstractTTN; init=random_ttn(flatten_external_indsnetwork(tn1, tn2); link_space=trivial_space(tn1)), nsweeps=1, + nsites=2, # used to be default of call to default_sweep_regions + updater_kwargs=(;), kwargs..., ) n = nv(tn1) @@ -42,7 +35,8 @@ function contract( ## end PH = ProjTTNApply(tn2, tn1) - psi = alternating_update(contract_solver, PH, init; nsweeps, kwargs...) + sweep_plan = default_sweep_regions(nsites, init; kwargs...) + psi = alternating_update(contract_updater, PH, init; nsweeps, sweep_plan, updater_kwargs, kwargs...) return psi end diff --git a/test/test_treetensornetworks/test_solvers/test_contract.jl b/test/test_treetensornetworks/test_solvers/test_contract.jl index 2643b5db..810db6fb 100644 --- a/test/test_treetensornetworks/test_solvers/test_contract.jl +++ b/test/test_treetensornetworks/test_solvers/test_contract.jl @@ -46,7 +46,7 @@ using Test # Test with nsite=1 Hpsi_guess = random_mps(t; internal_inds_space=32) - Hpsi = apply(H, psi; alg="fit", init=Hpsi_guess, nsite=1, nsweeps=4) + Hpsi = apply(H, psi; alg="fit", init=Hpsi_guess, nsites=1, nsweeps=4) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-4 end @@ -84,7 +84,7 @@ end # Test with nsite=1 Hpsi_guess = random_ttn(t; link_space=4) - Hpsi = apply(H, psi; alg="fit", nsite=1, nsweeps=4, init=Hpsi_guess) + Hpsi = apply(H, psi; alg="fit", nsites=1, nsweeps=4, init=Hpsi_guess) @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-4 end From 8bc62fb9dc2e6515c0a1ed2478b6040c2319ddc8 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 18:25:48 -0500 Subject: [PATCH 24/32] Format. --- src/solvers/contract.jl | 36 +++++++++--------- src/solvers/dmrg_x.jl | 5 +-- src/solvers/eigsolve.jl | 8 ++-- src/solvers/exponentiate.jl | 7 +--- .../solvers/alternating_update.jl | 10 ++--- src/treetensornetworks/solvers/contract.jl | 4 +- src/treetensornetworks/solvers/dmrg_x.jl | 1 - src/treetensornetworks/solvers/tdvp.jl | 17 +++++++-- src/treetensornetworks/solvers/update_step.jl | 38 +++++++++++++------ 9 files changed, 74 insertions(+), 52 deletions(-) diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index b3e3223d..fc2b176e 100644 --- a/src/solvers/contract.jl +++ b/src/solvers/contract.jl @@ -1,19 +1,19 @@ function contract_updater( - init; - state!, - projected_operator!, - outputlevel, - which_sweep, - sweep_plan, - which_region_update, - region_kwargs, - updater_kwargs, - ) - v = ITensor(1.0) - projected_operator = projected_operator![] - for j in sites(projected_operator) - v *= projected_operator.psi0[j] - end - Hpsi0 = contract(projected_operator, v) - return Hpsi0, (;) - end \ No newline at end of file + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs, +) + v = ITensor(1.0) + projected_operator = projected_operator![] + for j in sites(projected_operator) + v *= projected_operator.psi0[j] + end + Hpsi0 = contract(projected_operator, v) + return Hpsi0, (;) +end diff --git a/src/solvers/dmrg_x.jl b/src/solvers/dmrg_x.jl index 43ba69de..09618a9f 100644 --- a/src/solvers/dmrg_x.jl +++ b/src/solvers/dmrg_x.jl @@ -10,9 +10,8 @@ function dmrg_x_updater( updater_kwargs, ) # this updater does not seem to accept any kwargs? - default_updater_kwargs = (; - ) - updater_kwargs = merge(default_updater_kwargs, updater_kwargs) + default_updater_kwargs = (;) + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) H = contract(projected_operator![], ITensor(1.0)) D, U = eigen(H; ishermitian=true) u = uniqueind(U, H) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index eb9e8aef..aace5e09 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -20,17 +20,17 @@ function eigsolve_updater( ) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence howmany = 1 - which, updater_kwargs = _pop_which_eigenvalue(;updater_kwargs...) + which, updater_kwargs = _pop_which_eigenvalue(; updater_kwargs...) vals, vecs, info = eigsolve( projected_operator![], init, howmany, which; - updater_kwargs... #this leaves it + updater_kwargs..., #this leaves it ) return vecs[1], (; info, eigvals=vals) end -function _pop_which_eigenvalue(;which_eigenvalue, kwargs...) +function _pop_which_eigenvalue(; which_eigenvalue, kwargs...) return which_eigenvalue, NamedTuple(kwargs) -end \ No newline at end of file +end diff --git a/src/solvers/exponentiate.jl b/src/solvers/exponentiate.jl index 7d56ad39..a4dacebe 100644 --- a/src/solvers/exponentiate.jl +++ b/src/solvers/exponentiate.jl @@ -18,13 +18,10 @@ function exponentiate_updater( issymmetric=true, eager=true, ) - + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence result, exp_info = exponentiate( - projected_operator![], - region_kwargs.time_step, - init; - updater_kwargs... + projected_operator![], region_kwargs.time_step, init; updater_kwargs... ) return result, (; info=exp_info) end diff --git a/src/treetensornetworks/solvers/alternating_update.jl b/src/treetensornetworks/solvers/alternating_update.jl index a2fe2698..9c1ea8b6 100644 --- a/src/treetensornetworks/solvers/alternating_update.jl +++ b/src/treetensornetworks/solvers/alternating_update.jl @@ -65,11 +65,11 @@ function alternating_update( end projected_operator = disk(projected_operator) end - sweep_params=(; - maxdim=maxdim[which_sweep], - mindim=mindim[which_sweep], - cutoff=cutoff[which_sweep], - noise=noise[which_sweep], + sweep_params = (; + maxdim=maxdim[which_sweep], + mindim=mindim[which_sweep], + cutoff=cutoff[which_sweep], + noise=noise[which_sweep], ) sw_time = @elapsed begin state, projected_operator = sweep_update( diff --git a/src/treetensornetworks/solvers/contract.jl b/src/treetensornetworks/solvers/contract.jl index 92e4f744..398138e7 100644 --- a/src/treetensornetworks/solvers/contract.jl +++ b/src/treetensornetworks/solvers/contract.jl @@ -36,7 +36,9 @@ function contract( PH = ProjTTNApply(tn2, tn1) sweep_plan = default_sweep_regions(nsites, init; kwargs...) - psi = alternating_update(contract_updater, PH, init; nsweeps, sweep_plan, updater_kwargs, kwargs...) + psi = alternating_update( + contract_updater, PH, init; nsweeps, sweep_plan, updater_kwargs, kwargs... + ) return psi end diff --git a/src/treetensornetworks/solvers/dmrg_x.jl b/src/treetensornetworks/solvers/dmrg_x.jl index b5f8405b..4e89620e 100644 --- a/src/treetensornetworks/solvers/dmrg_x.jl +++ b/src/treetensornetworks/solvers/dmrg_x.jl @@ -20,4 +20,3 @@ end function dmrg_x(operator, init::AbstractTTN; updater=dmrg_x_updater, kwargs...) return dmrg_x(updater, operator, init; kwargs...) end - diff --git a/src/treetensornetworks/solvers/tdvp.jl b/src/treetensornetworks/solvers/tdvp.jl index 98f21849..f6081f46 100644 --- a/src/treetensornetworks/solvers/tdvp.jl +++ b/src/treetensornetworks/solvers/tdvp.jl @@ -72,7 +72,9 @@ function tdvp( kwargs..., ) nsweeps = _compute_nsweeps(nsteps, t, time_step, order) - sweep_plan = tdvp_sweep_plan(order, nsites, time_step, init_state; root_vertex, reverse_step) + sweep_plan = tdvp_sweep_plan( + order, nsites, time_step, init_state; root_vertex, reverse_step + ) function sweep_time_printer(; outputlevel, which_sweep, kwargs...) if outputlevel >= 1 @@ -88,7 +90,14 @@ function tdvp( insert_function!(sweep_observer!, "sweep_time_printer" => sweep_time_printer) state = alternating_update( - updater, operator, init_state; nsweeps, sweep_observer!, sweep_plan, updater_kwargs, kwargs... + updater, + operator, + init_state; + nsweeps, + sweep_observer!, + sweep_plan, + updater_kwargs, + kwargs..., ) # remove sweep_time_printer from sweep_observer! @@ -115,6 +124,8 @@ Optional keyword arguments: * `observer` - object implementing the Observer interface which can perform measurements and stop early * `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations """ -function tdvp(operator, t::Number, init_state::AbstractTTN; updater=exponentiate_updater, kwargs...) +function tdvp( + operator, t::Number, init_state::AbstractTTN; updater=exponentiate_updater, kwargs... +) return tdvp(updater, operator, t, init_state; kwargs...) end diff --git a/src/treetensornetworks/solvers/update_step.jl b/src/treetensornetworks/solvers/update_step.jl index 445fb329..890a13fa 100644 --- a/src/treetensornetworks/solvers/update_step.jl +++ b/src/treetensornetworks/solvers/update_step.jl @@ -15,10 +15,19 @@ function default_sweep_regions(nsites, graph::AbstractGraph; kwargs...) ###move end function region_update_printer(; - cutoff, maxdim, mindim, outputlevel::Int=0, state, sweep_plan, spec, which_region_update, which_sweep,kwargs... + cutoff, + maxdim, + mindim, + outputlevel::Int=0, + state, + sweep_plan, + spec, + which_region_update, + which_sweep, + kwargs..., ) if outputlevel >= 2 - region=first(sweep_plan[which_region_update]) + region = first(sweep_plan[which_region_update]) @printf("Sweep %d, region=%s \n", which_sweep, region) print(" Truncated using") @printf(" cutoff=%.1E", cutoff) @@ -60,10 +69,10 @@ function sweep_update( "`alternating_update` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.exponentiate`, etc.", ) end - + for which_region_update in eachindex(sweep_plan) - (region, region_kwargs)=sweep_plan[which_region_update] - region_kwargs=merge(region_kwargs, sweep_params) # sweep params has precedence over step_kwargs + (region, region_kwargs) = sweep_plan[which_region_update] + region_kwargs = merge(region_kwargs, sweep_params) # sweep params has precedence over step_kwargs state, projected_operator = region_update( solver, projected_operator, @@ -79,7 +88,7 @@ function sweep_update( ) end - select!(region_observer!, Observers.DataFrames.Not("region_update_printer")) # remove update_printer + select!(region_observer!, Observers.DataFrames.Not("region_update_printer")) # remove update_printer # Just to be sure: normalize && normalize!(state) @@ -173,16 +182,16 @@ function region_update( region_observer!, #insertion_kwargs, #ToDo: later #extraction_kwargs, #ToDo: implement later with possibility to pass custom extraction/insertion func (or code into func) - updater_kwargs + updater_kwargs, ) - region=first(sweep_plan[which_region_update]) + region = first(sweep_plan[which_region_update]) state = orthogonalize(state, current_ortho(region)) state, phi = extract_local_tensor(state, region;) nsites = (region isa AbstractEdge) ? 0 : length(region) #ToDo move into separate funtion projected_operator = set_nsite(projected_operator, nsites) projected_operator = position(projected_operator, state, region) state! = Ref(state) # create references, in case solver does (out-of-place) modify PH or state - projected_operator! = Ref(projected_operator) + projected_operator! = Ref(projected_operator) phi, info = updater( phi; state!, @@ -192,7 +201,7 @@ function region_update( sweep_plan, which_region_update, region_kwargs, - updater_kwargs + updater_kwargs, ) # args passed by reference are supposed to be modified out of place state = state![] # dereference projected_operator = projected_operator![] @@ -210,10 +219,15 @@ function region_update( #end state, spec = insert_local_tensor( - state, phi, region; eigen_perturbation=drho, ortho, normalize, + state, + phi, + region; + eigen_perturbation=drho, + ortho, + normalize, maxdim=region_kwargs.maxdim, mindim=region_kwargs.mindim, - cutoff=region_kwargs.cutoff + cutoff=region_kwargs.cutoff, ) update!( From ac1a923f3b9cabe308da0077599fdba1a719699a Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 20:07:48 -0500 Subject: [PATCH 25/32] Fix tdvp_time_dependent tests. --- .../test_solvers/test_tdvp_time_dependent.jl | 74 +++++++++++++++++-- 1 file changed, 68 insertions(+), 6 deletions(-) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl b/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl index 4abe5613..2448b282 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl @@ -23,6 +23,29 @@ ode_kwargs = (; reltol=1e-8, abstol=1e-8) ω⃗ = [ω₁, ω₂] f⃗ = [t -> cos(ω * t) for ω in ω⃗] +ode_updater_kwargs=(;f=f⃗,solver_alg=ode_alg,ode_kwargs) + +function ode_updater( + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs, +) + time_step=region_kwargs.time_step + f⃗=updater_kwargs.f + ode_kwargs=updater_kwargs.ode_kwargs + solver_alg=updater_kwargs.solver_alg + H⃗₀=projected_operator![] + result, info = ode_solver( + -im * TimeDependentSum(f⃗, H⃗₀), time_step, init; solver_alg, ode_kwargs... + ) + return result, (; info) +end function tdvp_ode_solver(H⃗₀, ψ₀; time_step, kwargs...) psi_t, info = ode_solver( @@ -31,7 +54,9 @@ function tdvp_ode_solver(H⃗₀, ψ₀; time_step, kwargs...) return psi_t, (; info) end + krylov_kwargs = (; tol=1e-8, eager=true) +krylov_updater_kwargs=(;f=f⃗,krylov_kwargs) function krylov_solver(H⃗₀, ψ₀; time_step, ishermitian=false, issymmetric=false, kwargs...) psi_t, info = krylov_solver( @@ -45,6 +70,42 @@ function krylov_solver(H⃗₀, ψ₀; time_step, ishermitian=false, issymmetric return psi_t, (; info) end +function krylov_updater( + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs, +) + default_updater_kwargs = (; + ishermitian=false, + issymmetric=false, + ) + + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedenc + time_step=region_kwargs.time_step + f⃗=updater_kwargs.f + krylov_kwargs=updater_kwargs.krylov_kwargs + ishermitian=updater_kwargs.ishermitian + issymmetric=updater_kwargs.issymmetric + H⃗₀=projected_operator![] + + result, info = krylov_solver( + -im * TimeDependentSum(f⃗, H⃗₀), + time_step, + init; + krylov_kwargs..., + ishermitian, + issymmetric, + ) + return result, (; info) +end + + @testset "MPS: Time dependent Hamiltonian" begin n = 4 J₁ = 1.0 @@ -53,7 +114,7 @@ end time_step = 0.1 time_total = 1.0 - nsite = 2 + nsites = 2 maxdim = 100 cutoff = 1e-8 @@ -65,9 +126,9 @@ end ψ₀ = complex(mps(s; states=(j -> isodd(j) ? "↑" : "↓"))) - ψₜ_ode = tdvp(tdvp_ode_solver, H⃗₀, time_total, ψ₀; time_step, maxdim, cutoff, nsite) + ψₜ_ode = tdvp(ode_updater, H⃗₀, time_total, ψ₀; time_step, maxdim, cutoff, nsites, updater_kwargs=ode_updater_kwargs) - ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_total, ψ₀; time_step, cutoff, nsite) + ψₜ_krylov = tdvp(krylov_updater, H⃗₀, time_total, ψ₀; time_step, cutoff, nsites, updater_kwargs=krylov_updater_kwargs) ψₜ_full, _ = tdvp_ode_solver(contract.(H⃗₀), contract(ψ₀); time_step=time_total) @@ -84,6 +145,7 @@ end @test krylov_err < 1e-2 end + @testset "TTN: Time dependent Hamiltonian" begin tooth_lengths = fill(2, 3) root_vertex = (3, 2) @@ -96,7 +158,7 @@ end time_step = 0.1 time_total = 1.0 - nsite = 2 + nsites = 2 maxdim = 100 cutoff = 1e-8 @@ -108,9 +170,9 @@ end ψ₀ = TTN(ComplexF64, s, v -> iseven(sum(isodd.(v))) ? "↑" : "↓") - ψₜ_ode = tdvp(tdvp_ode_solver, H⃗₀, time_total, ψ₀; time_step, maxdim, cutoff, nsite) + ψₜ_ode = tdvp(ode_updater, H⃗₀, time_total, ψ₀; time_step, maxdim, cutoff, nsites, updater_kwargs=ode_updater_kwargs) - ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_total, ψ₀; time_step, cutoff, nsite) + ψₜ_krylov = tdvp(krylov_updater, H⃗₀, time_total, ψ₀; time_step, cutoff, nsites, updater_kwargs=krylov_updater_kwargs) ψₜ_full, _ = tdvp_ode_solver(contract.(H⃗₀), contract(ψ₀); time_step=time_total) From 2748927d8cee065763d16f0be9f9309a9a47f2bc Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 20:12:38 -0500 Subject: [PATCH 26/32] Format tests. --- .../test_solvers/test_dmrg.jl | 3 +- .../test_solvers/test_dmrg_x.jl | 12 +-- .../test_solvers/test_tdvp.jl | 55 ++++++++---- .../test_solvers/test_tdvp_time_dependent.jl | 86 +++++++++++++------ 4 files changed, 101 insertions(+), 55 deletions(-) diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg.jl b/test/test_treetensornetworks/test_solvers/test_dmrg.jl index 9bc3f608..870012a0 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg.jl @@ -74,8 +74,7 @@ end sweep(; which_sweep, kw...) = which_sweep sweep_observer! = observer(sweep) - region(; which_region_update, sweep_plan, kw...) = - first(sweep_plan[which_region_update]) + region(; which_region_update, sweep_plan, kw...) = first(sweep_plan[which_region_update]) energy(; eigvals, kw...) = eigvals[1] region_observer! = observer(region, sweep, energy) diff --git a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl index 1353b8b8..4bb16268 100644 --- a/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl +++ b/test/test_treetensornetworks/test_solvers/test_dmrg_x.jl @@ -16,9 +16,7 @@ using Test ψ = mps(s; states=(v -> rand(["↑", "↓"]))) - dmrg_x_kwargs = ( - nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0 - ) + dmrg_x_kwargs = (nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0) ϕ = dmrg_x(H, ψ; nsites=2, dmrg_x_kwargs...) @@ -34,7 +32,7 @@ using Test # @test abs(loginner(ϕ̃, ϕ) / n) ≈ 0.0 atol = 1e-6 end -@testset "Tree DMRG-X" for conserve_qns in (false, true) +@testset "Tree DMRG-X" for conserve_qns in (false, true) tooth_lengths = fill(2, 3) root_vertex = (3, 2) c = named_comb_tree(tooth_lengths) @@ -51,10 +49,8 @@ end # TODO: Use `TTN(s; states=v -> rand(["↑", "↓"]))` or # `ttns(s; states=v -> rand(["↑", "↓"]))` ψ = normalize!(TTN(s, v -> rand(["↑", "↓"]))) - - dmrg_x_kwargs = ( - nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0 - ) + + dmrg_x_kwargs = (nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0) ϕ = dmrg_x(H, ψ; nsites=2, dmrg_x_kwargs...) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index d53c081c..e982f81d 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -235,10 +235,7 @@ using Test state; cutoff, normalize=false, - updater_kwargs=(; - tol=1e-12, - maxiter=500, - krylovdim=25,) + updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25), ) # TODO: What should `expect` output? Right now # it outputs a dictionary. @@ -250,10 +247,7 @@ using Test psi2; cutoff, normalize=false, - updater_kwargs=(; - tol=1e-12, - maxiter=500, - krylovdim=25,), + updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25), updater=exponentiate_updater, ) # TODO: What should `expect` output? Right now @@ -318,7 +312,14 @@ using Test nsites = (step <= 3 ? 2 : 1) phi = tdvp( - H, -tau * im, phi; nsteps=1, cutoff, nsites, normalize=true, updater_kwargs=(;krylovdim=15) + H, + -tau * im, + phi; + nsteps=1, + cutoff, + nsites, + normalize=true, + updater_kwargs=(; krylovdim=15), ) Sz1[step] = real(expect("Sz", state; vertices=[c])[c]) @@ -378,7 +379,14 @@ using Test for (step, t) in enumerate(trange) nsites = (step <= 10 ? 2 : 1) state = tdvp( - H, -tau, state; cutoff, nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) + H, + -tau, + state; + cutoff, + nsites, + reverse_step, + normalize=true, + updater_kwargs=(; krylovdim=15), ) #Different backend updaters, default updater_backend = "applyexp" psi2 = tdvp( @@ -389,7 +397,7 @@ using Test nsites, reverse_step, normalize=true, - updater_kwargs=(;krylovdim=15), + updater_kwargs=(; krylovdim=15), updater=ITensorNetworks.exponentiate_updater, ) end @@ -588,7 +596,7 @@ end # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end -=# + =# @testset "Accuracy Test" begin tau = 0.1 ttotal = 1.0 @@ -625,10 +633,7 @@ end state; cutoff, normalize=false, - updater_kwargs=(; - tol=1e-12, - maxiter=500, - krylovdim=25,) + updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25), ) push!(Sz_tdvp, real(expect("Sz", state; vertices=[c])[c])) push!(Sz_exact, real(scalar(dag(prime(statex, s[c])) * Szc * statex))) @@ -686,7 +691,14 @@ end nsites = (step <= 3 ? 2 : 1) phi = tdvp( - H, -tau * im, phi; nsteps=1, cutoff, nsites, normalize=true, updater_kwargs=(;krylovdim=15) + H, + -tau * im, + phi; + nsteps=1, + cutoff, + nsites, + normalize=true, + updater_kwargs=(; krylovdim=15), ) Sz1[step] = real(expect("Sz", state; vertices=[c])[c]) @@ -737,7 +749,14 @@ end for (step, t) in enumerate(trange) nsites = (step <= 10 ? 2 : 1) state = tdvp( - H, -tau, state; cutoff, nsites, reverse_step, normalize=true, updater_kwargs=(;krylovdim=15) + H, + -tau, + state; + cutoff, + nsites, + reverse_step, + normalize=true, + updater_kwargs=(; krylovdim=15), ) end diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl b/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl index 2448b282..763c9df2 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp_time_dependent.jl @@ -23,7 +23,7 @@ ode_kwargs = (; reltol=1e-8, abstol=1e-8) ω⃗ = [ω₁, ω₂] f⃗ = [t -> cos(ω * t) for ω in ω⃗] -ode_updater_kwargs=(;f=f⃗,solver_alg=ode_alg,ode_kwargs) +ode_updater_kwargs = (; f=f⃗, solver_alg=ode_alg, ode_kwargs) function ode_updater( init; @@ -36,12 +36,12 @@ function ode_updater( region_kwargs, updater_kwargs, ) - time_step=region_kwargs.time_step - f⃗=updater_kwargs.f - ode_kwargs=updater_kwargs.ode_kwargs - solver_alg=updater_kwargs.solver_alg - H⃗₀=projected_operator![] - result, info = ode_solver( + time_step = region_kwargs.time_step + f⃗ = updater_kwargs.f + ode_kwargs = updater_kwargs.ode_kwargs + solver_alg = updater_kwargs.solver_alg + H⃗₀ = projected_operator![] + result, info = ode_solver( -im * TimeDependentSum(f⃗, H⃗₀), time_step, init; solver_alg, ode_kwargs... ) return result, (; info) @@ -54,9 +54,8 @@ function tdvp_ode_solver(H⃗₀, ψ₀; time_step, kwargs...) return psi_t, (; info) end - krylov_kwargs = (; tol=1e-8, eager=true) -krylov_updater_kwargs=(;f=f⃗,krylov_kwargs) +krylov_updater_kwargs = (; f=f⃗, krylov_kwargs) function krylov_solver(H⃗₀, ψ₀; time_step, ishermitian=false, issymmetric=false, kwargs...) psi_t, info = krylov_solver( @@ -81,20 +80,17 @@ function krylov_updater( region_kwargs, updater_kwargs, ) - default_updater_kwargs = (; - ishermitian=false, - issymmetric=false, - ) + default_updater_kwargs = (; ishermitian=false, issymmetric=false) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedenc - time_step=region_kwargs.time_step - f⃗=updater_kwargs.f - krylov_kwargs=updater_kwargs.krylov_kwargs - ishermitian=updater_kwargs.ishermitian - issymmetric=updater_kwargs.issymmetric - H⃗₀=projected_operator![] - - result, info = krylov_solver( + time_step = region_kwargs.time_step + f⃗ = updater_kwargs.f + krylov_kwargs = updater_kwargs.krylov_kwargs + ishermitian = updater_kwargs.ishermitian + issymmetric = updater_kwargs.issymmetric + H⃗₀ = projected_operator![] + + result, info = krylov_solver( -im * TimeDependentSum(f⃗, H⃗₀), time_step, init; @@ -105,7 +101,6 @@ function krylov_updater( return result, (; info) end - @testset "MPS: Time dependent Hamiltonian" begin n = 4 J₁ = 1.0 @@ -126,9 +121,28 @@ end ψ₀ = complex(mps(s; states=(j -> isodd(j) ? "↑" : "↓"))) - ψₜ_ode = tdvp(ode_updater, H⃗₀, time_total, ψ₀; time_step, maxdim, cutoff, nsites, updater_kwargs=ode_updater_kwargs) + ψₜ_ode = tdvp( + ode_updater, + H⃗₀, + time_total, + ψ₀; + time_step, + maxdim, + cutoff, + nsites, + updater_kwargs=ode_updater_kwargs, + ) - ψₜ_krylov = tdvp(krylov_updater, H⃗₀, time_total, ψ₀; time_step, cutoff, nsites, updater_kwargs=krylov_updater_kwargs) + ψₜ_krylov = tdvp( + krylov_updater, + H⃗₀, + time_total, + ψ₀; + time_step, + cutoff, + nsites, + updater_kwargs=krylov_updater_kwargs, + ) ψₜ_full, _ = tdvp_ode_solver(contract.(H⃗₀), contract(ψ₀); time_step=time_total) @@ -145,7 +159,6 @@ end @test krylov_err < 1e-2 end - @testset "TTN: Time dependent Hamiltonian" begin tooth_lengths = fill(2, 3) root_vertex = (3, 2) @@ -170,9 +183,28 @@ end ψ₀ = TTN(ComplexF64, s, v -> iseven(sum(isodd.(v))) ? "↑" : "↓") - ψₜ_ode = tdvp(ode_updater, H⃗₀, time_total, ψ₀; time_step, maxdim, cutoff, nsites, updater_kwargs=ode_updater_kwargs) + ψₜ_ode = tdvp( + ode_updater, + H⃗₀, + time_total, + ψ₀; + time_step, + maxdim, + cutoff, + nsites, + updater_kwargs=ode_updater_kwargs, + ) - ψₜ_krylov = tdvp(krylov_updater, H⃗₀, time_total, ψ₀; time_step, cutoff, nsites, updater_kwargs=krylov_updater_kwargs) + ψₜ_krylov = tdvp( + krylov_updater, + H⃗₀, + time_total, + ψ₀; + time_step, + cutoff, + nsites, + updater_kwargs=krylov_updater_kwargs, + ) ψₜ_full, _ = tdvp_ode_solver(contract.(H⃗₀), contract(ψ₀); time_step=time_total) From 1fe92b30cdef567109152f68d4aa1157301aa2ba Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 20:38:25 -0500 Subject: [PATCH 27/32] Remove obsolete tests from test_tdvp (mostly those with an alternative updater backend. --- .../test_solvers/test_tdvp.jl | 141 +----------------- 1 file changed, 2 insertions(+), 139 deletions(-) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index e982f81d..dd69daa0 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -27,15 +27,6 @@ using Test # Time evolve forward: ψ1 = tdvp(H, -0.1im, ψ0; nsteps=1, cutoff, nsites=1) - # - #TODO: applyexp is removed, no need to test a different backend here and below for now. - # - #Different backend updaters, default updater_backend = "applyexp" - ψ1_exponentiate_backend = tdvp( - H, -0.1im, ψ0; nsteps=1, cutoff, nsites=1, updater=exponentiate_updater - ) - @test ψ1 ≈ ψ1_exponentiate_backend rtol = 1e-7 - @test norm(ψ1) ≈ 1.0 ## Should lose fidelity: @@ -47,12 +38,6 @@ using Test # Time evolve backwards: ψ2 = tdvp(H, +0.1im, ψ1; nsteps=1, cutoff) - #Different backend updaters, default updater_backend = "applyexp" - ψ2_exponentiate_backend = tdvp( - H, +0.1im, ψ1; nsteps=1, cutoff, updater=exponentiate_updater - ) - @test ψ2 ≈ ψ2_exponentiate_backend rtol = 1e-7 - @test norm(ψ2) ≈ 1.0 # Should rotate back to original state: @@ -83,12 +68,6 @@ using Test ψ1 = tdvp(Hs, -0.1im, ψ0; nsteps=1, cutoff, nsites=1) - #Different backend updaters, default updater_backend = "applyexp" - ψ1_exponentiate_backend = tdvp( - Hs, -0.1im, ψ0; nsteps=1, cutoff, nsites=1, updater=exponentiate_updater - ) - @test ψ1 ≈ ψ1_exponentiate_backend rtol = 1e-7 - @test norm(ψ1) ≈ 1.0 ## Should lose fidelity: @@ -100,12 +79,6 @@ using Test # Time evolve backwards: ψ2 = tdvp(Hs, +0.1im, ψ1; nsteps=1, cutoff) - #Different backend updaters, default updater_backend = "applyexp" - ψ2_exponentiate_backend = tdvp( - Hs, +0.1im, ψ1; nsteps=1, cutoff, updater=exponentiate_updater - ) - @test ψ2 ≈ ψ2_exponentiate_backend rtol = 1e-7 - @test norm(ψ2) ≈ 1.0 # Should rotate back to original state: @@ -146,54 +119,7 @@ using Test # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - #= - @testset_broken "Custom updater in TDVP" begin - N = 10 - cutoff = 1e-12 - - s = siteinds("S=1/2", N) - - os = OpSum() - for j in 1:(N - 1) - os += 0.5, "S+", j, "S-", j + 1 - os += 0.5, "S-", j, "S+", j + 1 - os += "Sz", j, "Sz", j + 1 - end - - H = mpo(os, s) - - ψ0 = random_mps(s; internal_inds_space=10) - - function updater(psi0; (state!), (projected_operator!), kwargs...) - updater_kwargs = (; - ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true - ) - time_step=get(step_kwargs,:time_step,nothing) - @assert !isnothing(time_step) - PH = projected_operator![] - state, exp_info = exponentiate(PH, time_step, psi0; updater_kwargs...) - return state, (; info=exp_info) - end - - ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsites=1) - - @test norm(ψ1) ≈ 1.0 - - ## Should lose fidelity: - #@test abs(inner(ψ0,ψ1)) < 0.9 - - # Average energy should be conserved: - @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0) - - # Time evolve backwards: - ψ2 = tdvp(H, +0.1im, ψ1; cutoff) - - @test norm(ψ2) ≈ 1.0 - - # Should rotate back to original state: - @test abs(inner(ψ0, ψ2)) > 0.99 - end - =# + @testset "Accuracy Test" begin N = 4 tau = 0.1 @@ -388,22 +314,8 @@ using Test normalize=true, updater_kwargs=(; krylovdim=15), ) - #Different backend updaters, default updater_backend = "applyexp" - psi2 = tdvp( - H, - -tau, - psi2; - cutoff, - nsites, - reverse_step, - normalize=true, - updater_kwargs=(; krylovdim=15), - updater=ITensorNetworks.exponentiate_updater, - ) end - @test state ≈ psi2 rtol = 1e-6 - en1 = inner(state', H, state) en2 = inner(psi2', H, psi2) @test en1 < -4.25 @@ -463,9 +375,6 @@ using Test En2_step = region_obs.En infos = region_obs.info - #@show sweep_obs - #@show step_obs - # # Could use ideas of other things to test here # @@ -550,53 +459,7 @@ end # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - # ToDo: Discuss whether the test commented out here is necessary given the new design? - #= - @testset "Custom updater in TDVP" begin - cutoff = 1e-12 - - tooth_lengths = fill(2, 3) - root_vertex = (3, 2) - c = named_comb_tree(tooth_lengths) - s = siteinds("S=1/2", c) - - os = ITensorNetworks.heisenberg(c) - - H = TTN(os, s) - - ψ0 = normalize!(random_ttn(s; link_space=10)) - - function updater(psi0; (state!), (projected_operator!), time_step, kwargs...) - updater_kwargs = (; - ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true - ) - PH = projected_operator![] - psi, exp_info = exponentiate(PH, time_step, psi0; updater_kwargs...) - return psi, (; info=exp_info) - end - - ψ1 = tdvp(updater, H, -0.1im, ψ0; cutoff, nsites=1) - - #@test ψ1 ≈ tdvp(updater, -0.1im, H, ψ0; cutoff, nsites=1) - #@test ψ1 ≈ tdvp(updater, H, ψ0, -0.1im; cutoff, nsites=1) - - @test norm(ψ1) ≈ 1.0 - - ## Should lose fidelity: - #@test abs(inner(ψ0,ψ1)) < 0.9 - - # Average energy should be conserved: - @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0) - - # Time evolve backwards: - ψ2 = tdvp(H, +0.1im, ψ1; cutoff) - - @test norm(ψ2) ≈ 1.0 - - # Should rotate back to original state: - @test abs(inner(ψ0, ψ2)) > 0.99 - end - =# + @testset "Accuracy Test" begin tau = 0.1 ttotal = 1.0 From 9e8e635e1f816f1c3c695ad9402420e5f606b0b4 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 20:38:44 -0500 Subject: [PATCH 28/32] Remove exponentiate from imports from KrylovKit --- src/imports.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/imports.jl b/src/imports.jl index 28bf0a88..71406878 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -34,7 +34,7 @@ import .DataGraphs: import Graphs: SimpleGraph, is_directed, weights -import KrylovKit: eigsolve, linsolve, exponentiate +import KrylovKit: eigsolve, linsolve import LinearAlgebra: factorize, normalize, normalize!, qr, svd From 23cc545abd7049487a670d5845d89d85d6949784 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 20:39:02 -0500 Subject: [PATCH 29/32] Apply review suggestions. --- src/solvers/contract.jl | 6 +++--- src/solvers/dmrg_x.jl | 2 +- src/solvers/eigsolve.jl | 9 +++++---- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index fc2b176e..588d35ae 100644 --- a/src/solvers/contract.jl +++ b/src/solvers/contract.jl @@ -9,11 +9,11 @@ function contract_updater( region_kwargs, updater_kwargs, ) - v = ITensor(1.0) + v = ITensor(true) projected_operator = projected_operator![] for j in sites(projected_operator) v *= projected_operator.psi0[j] end - Hpsi0 = contract(projected_operator, v) - return Hpsi0, (;) + vp = contract(projected_operator, v) + return vp, (;) end diff --git a/src/solvers/dmrg_x.jl b/src/solvers/dmrg_x.jl index 09618a9f..f1054726 100644 --- a/src/solvers/dmrg_x.jl +++ b/src/solvers/dmrg_x.jl @@ -12,7 +12,7 @@ function dmrg_x_updater( # this updater does not seem to accept any kwargs? default_updater_kwargs = (;) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) - H = contract(projected_operator![], ITensor(1.0)) + H = contract(projected_operator![], ITensor(true)) D, U = eigen(H; ishermitian=true) u = uniqueind(U, H) max_overlap, max_ind = findmax(abs, array(dag(init) * U)) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index aace5e09..abfa3a85 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -10,7 +10,7 @@ function eigsolve_updater( updater_kwargs, ) default_updater_kwargs = (; - which_eigenvalue=:SR, + which_eigval=:SR, ishermitian=true, tol=1e-14, krylovdim=3, @@ -20,13 +20,14 @@ function eigsolve_updater( ) updater_kwargs = merge(default_updater_kwargs, updater_kwargs) #last collection has precedence howmany = 1 - which, updater_kwargs = _pop_which_eigenvalue(; updater_kwargs...) + (; which_eigval) = updater_kwargs + updater_kwargs = Base.structdiff(updater_kwargs, (; which_eigval=nothing)) vals, vecs, info = eigsolve( projected_operator![], init, howmany, - which; - updater_kwargs..., #this leaves it + which_eigval; + updater_kwargs..., ) return vecs[1], (; info, eigvals=vals) end From 3d5ade8ed256f4719661cd0d91b5999c0bba600c Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 20:45:04 -0500 Subject: [PATCH 30/32] Fix test_tdvp.jl --- test/test_treetensornetworks/test_solvers/test_tdvp.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index dd69daa0..59d8cf2f 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -300,7 +300,6 @@ using Test H = mpo(os, s) state = random_mps(s; internal_inds_space=2) - psi2 = deepcopy(state) trange = 0.0:tau:ttotal for (step, t) in enumerate(trange) nsites = (step <= 10 ? 2 : 1) @@ -317,9 +316,7 @@ using Test end en1 = inner(state', H, state) - en2 = inner(psi2', H, psi2) @test en1 < -4.25 - @test en1 ≈ en2 end @testset "Observers" begin From 64ef223c2b95f10810963804b437cabada2d011a Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 20:49:30 -0500 Subject: [PATCH 31/32] Format. --- src/solvers/eigsolve.jl | 6 +----- test/test_treetensornetworks/test_solvers/test_tdvp.jl | 4 ++-- 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/solvers/eigsolve.jl b/src/solvers/eigsolve.jl index abfa3a85..85e99b3f 100644 --- a/src/solvers/eigsolve.jl +++ b/src/solvers/eigsolve.jl @@ -23,11 +23,7 @@ function eigsolve_updater( (; which_eigval) = updater_kwargs updater_kwargs = Base.structdiff(updater_kwargs, (; which_eigval=nothing)) vals, vecs, info = eigsolve( - projected_operator![], - init, - howmany, - which_eigval; - updater_kwargs..., + projected_operator![], init, howmany, which_eigval; updater_kwargs... ) return vecs[1], (; info, eigvals=vals) end diff --git a/test/test_treetensornetworks/test_solvers/test_tdvp.jl b/test/test_treetensornetworks/test_solvers/test_tdvp.jl index 59d8cf2f..d5d1ec49 100644 --- a/test/test_treetensornetworks/test_solvers/test_tdvp.jl +++ b/test/test_treetensornetworks/test_solvers/test_tdvp.jl @@ -119,7 +119,7 @@ using Test # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - + @testset "Accuracy Test" begin N = 4 tau = 0.1 @@ -456,7 +456,7 @@ end # Should rotate back to original state: @test abs(inner(ψ0, ψ2)) > 0.99 end - + @testset "Accuracy Test" begin tau = 0.1 ttotal = 1.0 From 812de92d7fa113108c4f18891f6d89594969dcb2 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Sat, 20 Jan 2024 22:00:09 -0500 Subject: [PATCH 32/32] Adapt linsolve. --- src/ITensorNetworks.jl | 1 + src/solvers/linsolve.jl | 22 ++++++++++ src/treetensornetworks/solvers/linsolve.jl | 40 +++++-------------- .../test_solvers/test_linsolve.jl | 4 +- 4 files changed, 36 insertions(+), 31 deletions(-) create mode 100644 src/solvers/linsolve.jl diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl index 2f3b4e03..9d5f14d2 100644 --- a/src/ITensorNetworks.jl +++ b/src/ITensorNetworks.jl @@ -110,6 +110,7 @@ include(joinpath("solvers", "eigsolve.jl")) include(joinpath("solvers", "exponentiate.jl")) include(joinpath("solvers", "dmrg_x.jl")) include(joinpath("solvers", "contract.jl")) +include(joinpath("solvers", "linsolve.jl")) include(joinpath("treetensornetworks", "abstracttreetensornetwork.jl")) include(joinpath("treetensornetworks", "ttn.jl")) include(joinpath("treetensornetworks", "opsum_to_ttn.jl")) diff --git a/src/solvers/linsolve.jl b/src/solvers/linsolve.jl new file mode 100644 index 00000000..1a595950 --- /dev/null +++ b/src/solvers/linsolve.jl @@ -0,0 +1,22 @@ +function linsolve_updater( + init; + state!, + projected_operator!, + outputlevel, + which_sweep, + sweep_plan, + which_region_update, + region_kwargs, + updater_kwargs, +) + default_updater_kwargs = (; + ishermitian=false, tol=1E-14, krylovdim=30, maxiter=100, verbosity=0, a₀, a₁ + ) + updater_kwargs = merge(default_updater_kwargs, updater_kwargs) + P = projected_operator![] + (; a₀, a₁) = updater_kwargs + updater_kwargs = Base.structdiff(updater_kwargs, (; a₀=nothing, a₁=nothing)) + b = dag(only(proj_mps(P))) + x, info = KrylovKit.linsolve(P, b, init, a₀, a₁; updater_kwargs...) + return x, (;) +end diff --git a/src/treetensornetworks/solvers/linsolve.jl b/src/treetensornetworks/solvers/linsolve.jl index 90ba8572..6f936020 100644 --- a/src/treetensornetworks/solvers/linsolve.jl +++ b/src/treetensornetworks/solvers/linsolve.jl @@ -28,39 +28,21 @@ function linsolve( x₀::AbstractTTN, a₀::Number=0, a₁::Number=1; - normalize, - region, - half_sweep, + updater=linsolve_updater, + nsweeps, #it makes sense to require this to be defined + nsites=2, + (sweep_observer!)=observer(), + root_vertex=default_root_vertex(init), + updater_kwargs=(;), + kwargs..., ) - function linsolve_solver( - P, - x₀; - ishermitian=false, - solver_tol=1E-14, - solver_krylovdim=30, - solver_maxiter=100, - solver_verbosity=0, - ) - b = dag(only(proj_mps(P))) - x, info = KrylovKit.linsolve( - P, - b, - x₀, - a₀, - a₁; - ishermitian, - tol=solver_tol, - krylovdim=solver_krylovdim, - maxiter=solver_maxiter, - verbosity=solver_verbosity, - ) - return x, NamedTuple() - end - + updater_kwargs = (; a₀, a₁, updater_kwargs...) error("`linsolve` for TTN not yet implemented.") + sweep_plan = default_sweep_regions(nsites, x0) # TODO: Define `itensornetwork_cache` # TODO: Define `linsolve_cache` + P = linsolve_cache(itensornetwork_cache(x₀', A, x₀), itensornetwork_cache(x₀', b)) - return alternating_update(linsolve_solver, P, x₀; kwargs...) + return alternating_update(linsolve_updater, P, x₀; sweep_plan, updater_kwargs, kwargs...) end diff --git a/test/test_treetensornetworks/test_solvers/test_linsolve.jl b/test/test_treetensornetworks/test_solvers/test_linsolve.jl index 68e8e953..72291642 100644 --- a/test/test_treetensornetworks/test_solvers/test_linsolve.jl +++ b/test/test_treetensornetworks/test_solvers/test_linsolve.jl @@ -41,11 +41,11 @@ using Random x_c = random_mps(s; states, internal_inds_space=4) + 0.1im * random_mps(s; states, internal_inds_space=2) - b = apply(H, x_c; cutoff) + b = apply(H, x_c; alg="fit", nsweeps=3) #cutoff is unsupported kwarg for apply/contract x0 = random_mps(s; states, internal_inds_space=10) x = @test_broken linsolve( - H, b, x0; cutoff, maxdim, nsweeps, ishermitian=true, solver_tol=1E-6 + H, b, x0; cutoff, maxdim, nsweeps, updater_kwargs=(; tol=1E-6, ishermitian=true) ) # @test norm(x - x_c) < 1E-3