diff --git a/Project.toml b/Project.toml index 15e422a..2fd21d7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorTDVP" uuid = "25707e16-a4db-4a07-99d9-4d67b7af0342" authors = ["Matthew Fishman and contributors"] -version = "0.4.0" +version = "0.4.1" [deps] Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" diff --git a/README.md b/README.md index 3130aa6..78c9263 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,23 @@ However, as noted above we now recommend installing and loading `ITensorMPS` ins ## News +### ITensorTDVP.jl v0.4.1 release notes + +#### Breaking changes + +#### New features + +- A new (experimental) `expand` function has been introduced for performing global Krylov expansion based on [arXiv:2005.06104](https://arxiv.org/abs/2005.06104), which can help with the accuracy of TDVP evolution in certain cases. See the docstrings of `expand` for more details: +```julia +julia> using ITensorTDVP + +julia> ? + +help?> expand +# ... +``` +Users are not given many customization options just yet as we gain more experience on the right balance between efficacy of the expansion and performance, and default values and keyword arguments are subject to change as we learn more about how to best use the method. + ### ITensorTDVP.jl v0.4 release notes #### Breaking changes diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl index 6bced2c..4457cd1 100644 --- a/src/ITensorTDVP.jl +++ b/src/ITensorTDVP.jl @@ -1,5 +1,5 @@ module ITensorTDVP -export TimeDependentSum, dmrg_x, linsolve, tdvp, to_vec +export TimeDependentSum, dmrg_x, expand, linsolve, tdvp, to_vec include("ITensorsExtensions.jl") using .ITensorsExtensions: to_vec include("applyexp.jl") @@ -17,6 +17,7 @@ include("contract.jl") include("reducedconstantterm.jl") include("reducedlinearproblem.jl") include("linsolve.jl") +include("expand.jl") using PackageExtensionCompat: @require_extensions function __init__() @require_extensions diff --git a/src/expand.jl b/src/expand.jl new file mode 100644 index 0000000..0d3c62a --- /dev/null +++ b/src/expand.jl @@ -0,0 +1,159 @@ +using ITensors: + ITensors, + Algorithm, + Index, + ITensor, + @Algorithm_str, + δ, + commonind, + dag, + denseblocks, + directsum, + hasqns, + prime, + scalartype, + uniqueinds +using ITensors.ITensorMPS: MPO, MPS, apply, dim, linkind, maxlinkdim, orthogonalize +using LinearAlgebra: normalize, svd, tr + +# Possible improvements: +# - Allow a maxdim argument to be passed to `expand`. +# - Current behavior is letting bond dimension get too +# big when used in imaginary time evolution. +# - Consider switching the default to variational/fit apply +# when building Krylov vectors. +# - Use (1-tau*operator)|state> to generate Krylov vectors +# instead of operator|state>. Is that needed? + +function expand(state, reference; alg, kwargs...) + return expand(Algorithm(alg), state, reference; kwargs...) +end + +function expand_cutoff_doctring() + return """ + The cutoff is used to control the truncation of the expanded + basis and defaults to half the precision of the scalar type + of the input state, i.e. ~1e-8 for `Float64`. + """ +end + +function expand_warning_doctring() + return """ + !!! warning + Users are not given many customization options just yet as we + gain more experience on the right balance between efficacy of the + expansion and performance in various scenarios, and default values + and keyword arguments are subject to change as we learn more about + how to best use the method. + """ +end + +function expand_citation_docstring() + return """ + [^global_expansion]: Time Dependent Variational Principle with Ancillary Krylov Subspace. Mingru Yang, Steven R. White, [arXiv:2005.06104](https://arxiv.org/abs/2005.06104) + """ +end + +""" + expand(state::MPS, references::Vector{MPS}; alg="orthogonalize", cutoff) + +Given an MPS `state` and a collection of MPS `references`, +returns an MPS which is equal to `state` +(has fidelity 1 with `state`) but whose MPS basis +is expanded to contain a portion of the basis of +the `references` that is orthogonal to the MPS basis of `state`. +See [^global_expansion] for more details. + +$(expand_cutoff_doctring()) + +$(expand_warning_doctring()) + +$(expand_citation_docstring()) +""" +function expand( + ::Algorithm"orthogonalize", + state::MPS, + references::Vector{MPS}; + cutoff=(√(eps(real(scalartype(state))))), +) + n = length(state) + state = orthogonalize(state, n) + references = map(reference -> orthogonalize(reference, n), references) + s = siteinds(state) + for j in reverse(2:n) + # SVD state[j] to compute basisⱼ + linds = [s[j - 1]; linkinds(state, j - 1)] + _, λⱼ, basisⱼ = svd(state[j], linds; righttags="bψ_$j,Link") + rinds = uniqueinds(basisⱼ, λⱼ) + # Make projectorⱼ + idⱼ = prod(r -> denseblocks(δ(scalartype(state), r', dag(r))), rinds) + projectorⱼ = idⱼ - prime(basisⱼ, rinds) * dag(basisⱼ) + # Sum reference density matrices + ρⱼ = sum(reference -> prime(reference[j], rinds) * dag(reference[j]), references) + # TODO: Fix bug that `tr` isn't preserving the element type. + ρⱼ /= scalartype(state)(tr(ρⱼ)) + # Apply projectorⱼ + ρⱼ_projected = apply(apply(projectorⱼ, ρⱼ), projectorⱼ) + expanded_basisⱼ = basisⱼ + if norm(ρⱼ_projected) > 10^3 * eps(real(scalartype(state))) + # Diagonalize projected density matrix ρⱼ_projected + # to compute reference_basisⱼ, which spans part of right basis + # of references which is orthogonal to right basis of state + dⱼ, reference_basisⱼ = eigen( + ρⱼ_projected; cutoff, ishermitian=true, righttags="bϕ_$j,Link" + ) + state_indⱼ = only(commoninds(basisⱼ, λⱼ)) + reference_indⱼ = only(commoninds(reference_basisⱼ, dⱼ)) + expanded_basisⱼ, expanded_indⱼ = directsum( + basisⱼ => state_indⱼ, reference_basisⱼ => reference_indⱼ + ) + end + # Shift ortho center one site left using dag(expanded_basisⱼ) + # and replace tensor at site j with expanded_basisⱼ + state[j - 1] = state[j - 1] * (state[j] * dag(expanded_basisⱼ)) + state[j] = expanded_basisⱼ + for reference in references + reference[j - 1] = reference[j - 1] * (reference[j] * dag(expanded_basisⱼ)) + reference[j] = expanded_basisⱼ + end + end + return state +end + +""" + expand(state::MPS, reference::MPO; alg="global_krylov", krylovdim=2, cutoff) + +Given an MPS `state` and an MPO `reference`, +returns an MPS which is equal to `state` +(has fidelity 1 with state) but whose MPS basis +is expanded to contain a portion of the basis of +the Krylov space built by repeated applications of +`reference` to `state` that is orthogonal +to the MPS basis of `state`. +The `reference` operator is applied to `state` `krylovdim` +number of times, with a default of 2 which should give +a good balance between reliability and performance. +See [^global_expansion] for more details. + +$(expand_cutoff_doctring()) + +$(expand_warning_doctring()) + +$(expand_citation_docstring()) +""" +function expand( + ::Algorithm"global_krylov", + state::MPS, + operator::MPO; + krylovdim=2, + cutoff=(√(eps(real(scalartype(state))))), + apply_kwargs=(; maxdim=maxlinkdim(state) + 1), +) + # TODO: Try replacing this logic with `Base.accumulate`. + references = Vector{MPS}(undef, krylovdim) + for k in 1:krylovdim + previous_reference = get(references, k - 1, state) + references[k] = normalize(apply(operator, previous_reference; apply_kwargs...)) + end + return expand(state, references; alg="orthogonalize", cutoff) +end diff --git a/src/sweep_update.jl b/src/sweep_update.jl index bdaf3d3..ed97685 100644 --- a/src/sweep_update.jl +++ b/src/sweep_update.jl @@ -1,6 +1,14 @@ using ITensors: ITensors, uniqueinds using ITensors.ITensorMPS: - ITensorMPS, MPS, isortho, orthocenter, orthogonalize!, position!, replacebond!, set_nsite! + ITensorMPS, + MPS, + isortho, + noiseterm, + orthocenter, + orthogonalize!, + position!, + replacebond!, + set_nsite! using LinearAlgebra: norm, normalize!, svd using Printf: @printf diff --git a/test/Project.toml b/test/Project.toml index c03e08e..7e322e9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,5 +8,6 @@ Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/test_expand.jl b/test/test_expand.jl new file mode 100644 index 0000000..34d8072 --- /dev/null +++ b/test/test_expand.jl @@ -0,0 +1,94 @@ +@eval module $(gensym()) +using ITensors: scalartype +using ITensors.ITensorMPS: OpSum, MPO, MPS, inner, linkdims, maxlinkdim, randomMPS, siteinds +using ITensorTDVP: dmrg, expand, tdvp +using LinearAlgebra: normalize +using StableRNGs: StableRNG +using Test: @test, @testset +const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) +@testset "expand (eltype=$elt)" for elt in elts + @testset "expand (alg=\"orthogonalize\", conserve_qns=$conserve_qns, eltype=$elt)" for conserve_qns in + ( + false, true + ) + n = 6 + s = siteinds("S=1/2", n; conserve_qns) + rng = StableRNG(1234) + state = randomMPS(rng, elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=4) + reference = randomMPS(rng, elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=2) + state_expanded = expand(state, [reference]; alg="orthogonalize") + @test scalartype(state_expanded) === elt + @test inner(state_expanded, state) ≈ inner(state, state) + @test inner(state_expanded, reference) ≈ inner(state, reference) + end + @testset "expand (alg=\"global_krylov\", conserve_qns=$conserve_qns, eltype=$elt)" for conserve_qns in + ( + false, true + ) + n = 10 + s = siteinds("S=1/2", n; conserve_qns) + opsum = OpSum() + for j in 1:(n - 1) + opsum += 0.5, "S+", j, "S-", j + 1 + opsum += 0.5, "S-", j, "S+", j + 1 + opsum += "Sz", j, "Sz", j + 1 + end + operator = MPO(elt, opsum, s) + state = MPS(elt, s, j -> isodd(j) ? "↑" : "↓") + state_expanded = expand(state, operator; alg="global_krylov") + @test scalartype(state_expanded) === elt + @test maxlinkdim(state_expanded) > 1 + @test inner(state_expanded, state) ≈ inner(state, state) + end + @testset "Decoupled ladder (alg=\"global_krylov\", eltype=$elt)" begin + nx = 10 + ny = 2 + n = nx * ny + s = siteinds("S=1/2", n) + opsum = OpSum() + for j in 1:2:(n - 2) + opsum += 1 / 2, "S+", j, "S-", j + 2 + opsum += 1 / 2, "S-", j, "S+", j + 2 + opsum += "Sz", j, "Sz", j + 2 + end + for j in 2:2:(n - 2) + opsum += 1 / 2, "S+", j, "S-", j + 2 + opsum += 1 / 2, "S-", j, "S+", j + 2 + opsum += "Sz", j, "Sz", j + 2 + end + operator = MPO(elt, opsum, s) + rng = StableRNG(1234) + init = randomMPS(rng, elt, s; linkdims=30) + reference_energy, reference_state = dmrg( + operator, + init; + nsweeps=15, + maxdim=[10, 10, 20, 20, 40, 80, 100], + cutoff=(√(eps(real(elt)))), + noise=(√(eps(real(elt)))), + ) + rng = StableRNG(1234) + state = randomMPS(rng, elt, s) + nexpansions = 10 + tau = elt(0.5) + for step in 1:nexpansions + # TODO: Use `fourthroot`/`∜` in Julia 1.10 and above. + state = expand( + state, operator; alg="global_krylov", krylovdim=3, cutoff=eps(real(elt))^(1//4) + ) + state = tdvp( + operator, + -4tau, + state; + nsteps=4, + cutoff=1e-5, + updater_kwargs=(; tol=1e-3, krylovdim=5), + ) + state = normalize(state) + end + @test scalartype(state) === elt + # TODO: Use `fourthroot`/`∜` in Julia 1.10 and above. + @test inner(state', operator, state) ≈ reference_energy rtol = 5 * eps(real(elt))^(1//4) + end +end +end diff --git a/test/test_exports.jl b/test/test_exports.jl index ffc1f85..ab457cd 100644 --- a/test/test_exports.jl +++ b/test/test_exports.jl @@ -4,7 +4,7 @@ using Test: @test, @testset @testset "Test exports" begin @test issetequal( names(ITensorTDVP), - [:ITensorTDVP, :TimeDependentSum, :dmrg_x, :linsolve, :tdvp, :to_vec], + [:ITensorTDVP, :TimeDependentSum, :dmrg_x, :expand, :linsolve, :tdvp, :to_vec], ) end end