Output eigenvalue from ITensorTDVP.dmrg and dmrg_x (#71)

ITensor · May 6, 2024 · 4669573 · 4669573 · mtfishman · May 6, 2024
1 parent 8ecdb83
commit 4669573
Show file tree

Hide file tree

Showing 11 changed files with 70 additions and 17 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "ITensorTDVP"
 uuid = "25707e16-a4db-4a07-99d9-4d67b7af0342"
 authors = ["Matthew Fishman <mfishman@flatironinstitute.org> and contributors"]
-version = "0.2.6"
+version = "0.3.0"
 
 [deps]
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"

diff --git a/README.md b/README.md
@@ -19,6 +19,21 @@ pkg> add ITensorTDVP
 
 ## News
 
+### ITensorTDVP.jl v0.3 Release Notes
+
+#### Breaking changes
+
+- `ITensorTDVP.dmrg` and `ITensorTDVP.dmrg_x` now output a tuple containing the eigenvalue and eigenvector, while before they just output the eigenvector. You should update code like this:
+```julia
+psi = dmrg_x(H, psi0; nsweeps=10, maxdim=100, cutoff=1e-6)
+psi = ITensorTDVP.dmrg(H, psi0; nsweeps=10, maxdim=100, cutoff=1e-6)
+```
+to:
+```julia
+energy, psi = dmrg_x(H, psi0; nsweeps=10, maxdim=100, cutoff=1e-6)
+energy, psi = ITensorTDVP.dmrg(H, psi0; nsweeps=10, maxdim=100, cutoff=1e-6)
+```
+
 ### ITensorTDVP.jl v0.2 Release Notes
 
 #### Breaking changes

diff --git a/examples/01_tdvp.jl b/examples/01_tdvp.jl
@@ -36,11 +36,11 @@ function main()
 
   e2, ϕ2 = ITensors.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10)
 
-  @show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2)
+  @show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2), e2
 
-  ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1)
+  e3, ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1)
 
-  @show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3)
+  @show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3), e3
   return nothing
 end
 

diff --git a/examples/02_dmrg-x.jl b/examples/02_dmrg-x.jl
@@ -36,11 +36,11 @@ function main()
     nsweeps=10, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1
   )
 
-  ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
+  e, ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
 
   @show inner(ψ', H, ψ) / inner(ψ, ψ)
   @show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2
-  @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
+  @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ), e
   @show inner(H, ϕ, H, ϕ) - inner(ϕ', H, ϕ)^2
   return nothing
 end

diff --git a/examples/Project.toml b/examples/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342"
+ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/src/dmrg.jl b/src/dmrg.jl
@@ -21,8 +21,7 @@ function dmrg_solver(
       maxiter=solver_maxiter,
       verbosity=solver_verbosity,
     )
-    psi = vecs[1]
-    return psi, info
+    return vecs[1], (; info, eigval=vals[1])
   end
   return solver
 end
@@ -36,9 +35,13 @@ function dmrg(
   solver_krylovdim=default_solver_krylovdim(eigsolve),
   solver_maxiter=default_solver_maxiter(eigsolve),
   solver_verbosity=default_solver_verbosity(),
+  (observer!)=default_observer!(),
   kwargs...,
 )
   reverse_step = false
+  info_ref! = Ref{Any}()
+  info_observer! = values_observer(; info=info_ref!)
+  observer! = compose_observers(observer!, info_observer!)
   psi = alternating_update(
     dmrg_solver(
       eigsolve;
@@ -52,7 +55,8 @@ function dmrg(
     H,
     psi0;
     reverse_step,
+    observer!,
     kwargs...,
   )
-  return psi
+  return info_ref![].eigval, psi
 end
diff --git a/src/dmrg_x.jl b/src/dmrg_x.jl
@@ -5,12 +5,19 @@ function dmrg_x_solver(PH, t, psi0; current_time, outputlevel)
   H = contract(PH, ITensor(true))
   D, U = eigen(H; ishermitian=true)
   u = uniqueind(U, H)
+  u′ = uniqueind(D, U)
   max_overlap, max_ind = findmax(abs, array(psi0 * dag(U)))
   U_max = U * dag(onehot(eltype(U), u => max_ind))
-  return U_max, nothing
+  D_max = D[u′ => max_ind, u => max_ind]
+  return U_max, (; eigval=D_max)
 end
 
-function dmrg_x(PH, psi0::MPS; reverse_step=false, kwargs...)
-  psi = alternating_update(dmrg_x_solver, PH, psi0; reverse_step, kwargs...)
-  return psi
+function dmrg_x(
+  PH, psi0::MPS; reverse_step=false, (observer!)=default_observer!(), kwargs...
+)
+  info_ref = Ref{Any}()
+  info_observer! = values_observer(; info=info_ref)
+  observer! = compose_observers(observer!, info_observer!)
+  psi = alternating_update(dmrg_x_solver, PH, psi0; reverse_step, observer!, kwargs...)
+  return info_ref[].eigval, psi
 end
diff --git a/src/update_observer.jl b/src/update_observer.jl
@@ -6,3 +6,25 @@ using ITensors.ITensorMPS: ITensorMPS
 function update_observer!(observer::ITensorMPS.AbstractObserver; kwargs...)
   return ITensorMPS.measure!(observer; kwargs...)
 end
+
+struct ValuesObserver{Values<:NamedTuple}
+  values::Values
+end
+function update_observer!(observer::ValuesObserver; kwargs...)
+  for key in keys(observer.values)
+    observer.values[key][] = kwargs[key]
+  end
+  return observer
+end
+values_observer(; kwargs...) = ValuesObserver(NamedTuple(kwargs))
+
+struct ComposedObservers{Observers<:Tuple}
+  observers::Observers
+end
+compose_observers(observers...) = ComposedObservers(observers)
+function update_observer!(observer::ComposedObservers; kwargs...)
+  for observerᵢ in observer.observers
+    update_observer!(observerᵢ; kwargs...)
+  end
+  return observer
+end
diff --git a/test/test_dmrg.jl b/test/test_dmrg.jl
@@ -21,12 +21,14 @@ using Test: @test, @testset
   psi = randomMPS(elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=20)
   nsweeps = 10
   maxdim = [10, 20, 40, 100]
-  psi = ITensorTDVP.dmrg(
+  e, psi = ITensorTDVP.dmrg(
     H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1
   )
+  @test inner(psi', H, psi) ≈ e
   e2, psi2 = ITensors.dmrg(H, psi; nsweeps, maxdim, cutoff, outputlevel=0)
   @test ITensors.scalartype(psi2) == elt
   @test e2 isa real(elt)
+  @test e ≈ e2 rtol = √(eps(real(elt))) * 10
   @test inner(psi', H, psi) ≈ inner(psi2', H, psi2) rtol = √(eps(real(elt))) * 10
 end
 end
diff --git a/test/test_dmrg_x.jl b/test/test_dmrg_x.jl
@@ -34,13 +34,15 @@ using Test: @test, @testset
   dmrg_x_kwargs = (
     nsweeps=20, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0
   )
-  ϕ = dmrg_x(ProjMPO(H), ψ; nsite=2, dmrg_x_kwargs...)
+  e, ϕ = dmrg_x(ProjMPO(H), ψ; nsite=2, dmrg_x_kwargs...)
   @test ITensors.scalartype(ϕ) == elt
+  @test inner(ϕ', H, ϕ) / inner(ϕ, ϕ) ≈ e
   @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = √eps(real(elt))
   @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1
   @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = √eps(real(elt))
-  ϕ̃ = dmrg_x(ProjMPO(H), ϕ; nsite=1, dmrg_x_kwargs...)
+  ẽ, ϕ̃ = dmrg_x(ProjMPO(H), ϕ; nsite=1, dmrg_x_kwargs...)
   @test ITensors.scalartype(ϕ̃) == elt
+  @test inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) ≈ ẽ
   @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) rtol = 1e-1
   scale(::Type{Float32}) = 10^2
   scale(::Type{Float64}) = 10^5

diff --git a/test/test_linsolve.jl b/test/test_linsolve.jl
@@ -22,7 +22,7 @@ using Random: Random
   state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
   Random.seed!(1234)
   x_c = randomMPS(elt, s, state; linkdims=2)
-  x_c = dmrg(H, x_c; nsweeps=10, cutoff=1e-6, maxdim=20, outputlevel=0)
+  e, x_c = dmrg(H, x_c; nsweeps=10, cutoff=1e-6, maxdim=20, outputlevel=0)
   @test ITensors.scalartype(x_c) == elt
   # Compute `b = H * x_c`
   b = apply(H, x_c; cutoff=1e-8)