Start adding infinite TDVP interface (#38)

* Start adding infinite TDVP interface by sharing code with VUMPS and swapping out the solver. VUMPS can be run with a TDVP time step of `-Inf`. Currently sequential TDVP gives incorrect results, and parallel TDVP is broken.
* Add better formatting Github Action.

Author: mtfishman

.github/workflows/format_check.yml

@@ -1,42 +1,35 @@
-name: format-check
-
+name: Format check
 on:
   push:
-    branches:
-      - 'main'
-      - 'release-'
-    tags: '*'
+    branches: [master]
+    tags: [v*]
   pull_request:
 
 jobs:
-  build:
-    runs-on: ${{ matrix.os }}
-    strategy:
-      matrix:
-        julia-version: [1.7.1]
-        julia-arch: [x86]
-        os: [ubuntu-latest]
+  format:
+    name: "Format Check"
+    runs-on: ubuntu-latest
     steps:
-      - uses: julia-actions/setup-julia@latest
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
         with:
-          version: ${{ matrix.julia-version }}
-
-      - uses: actions/checkout@v1
+          version: 1
       - name: Install JuliaFormatter and format
-        # This will use the latest version by default but you can set the version like so:
-        #
-        # julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))'
         run: |
           julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))'
           julia  -e 'using JuliaFormatter; format(".", verbose=true)'
-      - name: Format check
+      - name: Check format
         run: |
           julia -e '
           out = Cmd(`git diff --name-only`) |> read |> String
           if out == ""
               exit(0)
           else
-              @error "Some files have not been formatted !!!"
+              @error "The following files have not been formatted:"
               write(stdout, out)
+              out_diff = Cmd(`git diff`) |> read |> String
+              @error "Diff:"
+              write(stdout, out_diff)
               exit(1)
+              @error ""
           end'

.github/workflows/format_suggestions.yml

@@ -0,0 +1,27 @@
+name: Format suggestions
+
+on:
+  pull_request:
+
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+
+jobs:
+  format:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: 1
+      - run: |
+          julia  -e 'using Pkg; Pkg.add("JuliaFormatter")'
+          julia  -e 'using JuliaFormatter; format("."; verbose=true)'
+      - uses: reviewdog/action-suggester@v1
+        with:
+          tool_name: JuliaFormatter
+          fail_on_error: true
+          filter_mode: added

examples/vumps/vumps_ising.jl

@@ -9,9 +9,10 @@ maxdim = 5 # Maximum bond dimension
 cutoff = 1e-12 # Singular value cutoff when increasing the bond dimension
 max_vumps_iters = 10 # Maximum number of iterations of the VUMPS algorithm at a fixed bond dimension
 outer_iters = 10 # Number of times to increase the bond dimension
+time_step = -Inf # -Inf corresponds to VUMPS
 
 # Parameters of the transverse field Ising model
-model_params = (J=1.0, h=0.8)
+model_params = (J=1.0, h=0.9)
 
 ##############################################################################
 # CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
@@ -38,14 +39,14 @@ H = InfiniteITensorSum(model, s; model_params...)
 
 vumps_kwargs = (tol=1e-5, maxiter=max_vumps_iters)
 subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
-ψ = vumps(H, ψ; vumps_kwargs...)
+ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
 
 # Alternate steps of running VUMPS and increasing the bond dimension
 @time for _ in 1:outer_iters
   println("\nIncrease bond dimension")
   global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
   println("Run VUMPS with new bond dimension")
-  global ψ = vumps(H, ψ; vumps_kwargs...)
+  global ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
 end
 
 # Check translational invariance

src/ITensorInfiniteMPS.jl

@@ -18,7 +18,7 @@ using HDF5
 using QuadGK
 
 using ITensors.NDTensors: eachdiagblock
-using KrylovKit: eigsolve, linsolve
+using KrylovKit: eigsolve, linsolve, exponentiate
 
 import Base: getindex, length, setindex!, +, -, *
 
@@ -67,6 +67,7 @@ export Cell,
   reference,
   subspace_expansion,
   translatecell,
+  tdvp,
   vumps,
   ⊕,
   ⊗,

src/vumps_mpo.jl renamed to src/broken/vumps_mpo.jl

@@ -1,3 +1,18 @@
+#################################################################
+# XXX: This version of VUMPS using MPOs is not working. It was
+# the initial implementation of VUMPS I (@mtfishman) tried to make, but
+# there was a bug and it was temporarily abandoned in favor
+# of the simpler local Hamiltonian implementation found in
+# `src/vumps_localham.jl`. Help getting this working would be highly
+# appreciated so we can support more generic Hamiltonians.
+# Ideally, the implementation would only overload the minimal amount
+# of code needed to be modified (like computing quasi-left and
+# right environments and updates to `AC` and `C`), and share the high
+# level code with the existing VUMPS code in `src/vumps_localham.jl`.
+#
+# Please reach out on this issue: https://github.com/ITensor/ITensorInfiniteMPS.jl/issues/39
+# if you are interested in implementing this feature.
+#################################################################
 
 # Make an ITensorMap representing the transfer matrix T|v> = |Tv>
 function transfer_matrix(ψ::InfiniteMPS)

src/vumps_localham.jl

@@ -184,28 +184,30 @@ function right_environment(hᴿ, ψ; tol=1e-15)
   return Hᴿ
 end
 
-function vumps_iteration(args...; multisite_update_alg="sequential", kwargs...)
+function tdvp_iteration(args...; multisite_update_alg="sequential", kwargs...)
   if multisite_update_alg == "sequential"
-    return vumps_iteration_sequential(args...; kwargs...)
+    return tdvp_iteration_sequential(args...; kwargs...)
   elseif multisite_update_alg == "parallel"
-    return vumps_iteration_parallel(args...; kwargs...)
+    return tdvp_iteration_parallel(args...; kwargs...)
   else
     error(
       "Multisite update algorithm multisite_update_alg = $multisite_update_alg not supported, use \"parallel\" or \"sequential\"",
     )
   end
 end
 
-function vumps_iteration_sequential(
+function tdvp_iteration_sequential(
+  solver::Function,
   ∑h::InfiniteITensorSum,
   ψ::InfiniteCanonicalMPS;
   (ϵᴸ!)=fill(1e-15, nsites(ψ)),
   (ϵᴿ!)=fill(1e-15, nsites(ψ)),
-  eigsolve_tol=(x -> x / 100),
+  time_step,
+  solver_tol=(x -> x / 100),
 )
   Nsites = nsites(ψ)
   ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
-  krylov_tol = eigsolve_tol(ϵᵖʳᵉˢ)
+  _solver_tol = solver_tol(ϵᵖʳᵉˢ)
   ψᴴ = dag(ψ)
   ψ′ = ψᴴ'
   # XXX: make this prime the center sites
@@ -267,24 +269,23 @@ function vumps_iteration_sequential(
     for k in 1:Nsites
       𝕙ᴸ[k] = left_environment_cell(ψ, ψ̃, hᴸ, k)
     end
-    Hᴸ = left_environment(hᴸ, 𝕙ᴸ, ψ; tol=krylov_tol)
+    Hᴸ = left_environment(hᴸ, 𝕙ᴸ, ψ; tol=_solver_tol)
     for k in 2:Nsites
       hᴿ[k] = hᴿ[k + 1] * ψ.AR[k + 1] * ψ̃.AR[k + 1] + hᴿ[k]
     end
-    Hᴿ = right_environment(hᴿ, ψ; tol=krylov_tol)
+    Hᴿ = right_environment(hᴿ, ψ; tol=_solver_tol)
 
-    Cvalsₙ₋₁, Cvecsₙ₋₁, Cinfoₙ₋₁ = eigsolve(
-      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n - 1), ψ.C[n - 1], 1, :SR; ishermitian=true, tol=krylov_tol
+    Cvalsₙ₋₁, Cvecsₙ₋₁, Cinfoₙ₋₁ = solver(
+      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n - 1), time_step, ψ.C[n - 1], _solver_tol
     )
-    Cvalsₙ, Cvecsₙ, Cinfoₙ = eigsolve(
-      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), ψ.C[n], 1, :SR; ishermitian=true, tol=krylov_tol
+    Cvalsₙ, Cvecsₙ, Cinfoₙ = solver(Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.C[n], _solver_tol)
+    Avalsₙ, Avecsₙ, Ainfoₙ = solver(
+      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol
     )
-    Avalsₙ, Avecsₙ, Ainfoₙ = eigsolve(
-      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), ψ.AL[n] * ψ.C[n], 1, :SR; ishermitian=true, tol=krylov_tol
-    )
-    C̃[n - 1] = Cvecsₙ₋₁[1]
-    C̃[n] = Cvecsₙ[1]
-    Ãᶜ[n] = Avecsₙ[1]
+
+    C̃[n - 1] = Cvecsₙ₋₁
+    C̃[n] = Cvecsₙ
+    Ãᶜ[n] = Avecsₙ
 
     function ortho_overlap(AC, C)
       AL, _ = polar(AC * dag(C), uniqueinds(AC, C))
@@ -324,16 +325,17 @@ function vumps_iteration_sequential(
   return ψ, (eᴸ, eᴿ)
 end
 
-function vumps_iteration_parallel(
+function tdvp_iteration_parallel(
+  solver::Function,
   ∑h::InfiniteITensorSum,
   ψ::InfiniteCanonicalMPS;
   (ϵᴸ!)=fill(1e-15, nsites(ψ)),
   (ϵᴿ!)=fill(1e-15, nsites(ψ)),
-  eigsolve_tol=(x -> x / 100),
+  solver_tol=(x -> x / 100),
 )
   Nsites = nsites(ψ)
   ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
-  krylov_tol = ϵᵖʳᵉˢ / 100
+  _solver_tol = solver_tol(ϵᵖʳᵉˢ)
   ψᴴ = dag(ψ)
   ψ′ = ψᴴ'
   # XXX: make this prime the center sites
@@ -381,22 +383,21 @@ function vumps_iteration_parallel(
   for n in 2:Nsites
     hᴸ[n] = hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n]
   end
-  Hᴸ = left_environment(hᴸ, ψ; tol=krylov_tol)
+  Hᴸ = left_environment(hᴸ, ψ; tol=_solver_tol)
 
   for n in 2:Nsites
     hᴿ[n] = hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
   end
-  Hᴿ = right_environment(hᴿ, ψ; tol=krylov_tol)
+  Hᴿ = right_environment(hᴿ, ψ; tol=_solver_tol)
 
   C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
   Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
   for n in 1:Nsites
-    Cvalsₙ, Cvecsₙ, Cinfoₙ = eigsolve(
-      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), ψ.C[n], 1, :SR; ishermitian=true, tol=krylov_tol
-    )
-    Avalsₙ, Avecsₙ, Ainfoₙ = eigsolve(
-      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), ψ.AL[n] * ψ.C[n], 1, :SR; ishermitian=true, tol=krylov_tol
+    Cvalsₙ, Cvecsₙ, Cinfoₙ = solver(Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.C[n], _solver_tol)
+    Avalsₙ, Avecsₙ, Ainfoₙ = solver(
+      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol
     )
+
     C̃[n] = Cvecsₙ[1]
     Ãᶜ[n] = Avecsₙ[1]
   end
@@ -421,28 +422,32 @@ function vumps_iteration_parallel(
   return InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ), (eᴸ, eᴿ)
 end
 
-function vumps(
+function tdvp(
+  solver::Function,
   ∑h,
   ψ;
   maxiter=10,
   tol=1e-8,
   outputlevel=1,
   multisite_update_alg="sequential",
-  eigsolve_tol=(x -> x / 100),
+  solver_tol=(x -> x / 100),
+  time_step,
 )
   N = nsites(ψ)
   (ϵᴸ!) = fill(tol, nsites(ψ))
   (ϵᴿ!) = fill(tol, nsites(ψ))
   outputlevel > 0 &&
     println("Running VUMPS with multisite_update_alg = $multisite_update_alg")
   for iter in 1:maxiter
-    ψ, (eᴸ, eᴿ) = vumps_iteration(
+    ψ, (eᴸ, eᴿ) = tdvp_iteration(
+      solver,
       ∑h,
       ψ;
       (ϵᴸ!)=(ϵᴸ!),
       (ϵᴿ!)=(ϵᴿ!),
       multisite_update_alg=multisite_update_alg,
-      eigsolve_tol=eigsolve_tol,
+      solver_tol=solver_tol,
+      time_step=time_step,
     )
     ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
     maxdimψ = maxlinkdim(ψ[0:(N + 1)])
@@ -459,6 +464,38 @@ function vumps(
   return ψ
 end
 
+function vumps_solver(M, time_step, v₀, solver_tol)
+  λ⃗, v⃗, info = eigsolve(M, v₀, 1, :SR; ishermitian=true, tol=solver_tol)
+  return λ⃗[1], v⃗[1], info
+end
+
+return function tdvp_solver(M, time_step, v₀, solver_tol)
+  v, info = exponentiate(M, time_step, v₀; ishermitian=true, tol=solver_tol)
+  return nothing, v, info
+end
+
+function vumps(
+  args...; time_step=-Inf, eigsolve_tol=(x -> x / 100), solver_tol=eigsolve_tol, kwargs...
+)
+  @assert isinf(time_step) && time_step < 0
+  println("Using VUMPS solver with time step $time_step")
+  return tdvp(vumps_solver, args...; time_step=time_step, solver_tol=solver_tol, kwargs...)
+end
+
+function tdvp(args...; time_step, solver_tol=(x -> x / 100), kwargs...)
+  solver = if !isinf(time_step)
+    println("Using TDVP solver with time step $time_step")
+    tdvp_solver
+  elseif time_step < 0
+    # Call VUMPS instead
+    println("Using VUMPS solver with time step $time_step")
+    vumps_solver
+  else
+    error("Time step $time_step not supported.")
+  end
+  return tdvp(solver, args...; time_step=time_step, solver_tol=solver_tol, kwargs...)
+end
+
 ##################################################################
 # Old functionality, only used for testing