From 68d18c95a0317ed07b08ec61d995b6b54fbefd39 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Thu, 23 Oct 2025 11:29:47 -0400
Subject: [PATCH] Looser bounds in expand

---
 src/solvers/expand.jl                 | 18 ++++++++-
 test/base/test_solvers/test_expand.jl | 54 ++++++++++++++++++++++++++-
 2 files changed, 68 insertions(+), 4 deletions(-)

diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl
index 0e09ed3..9b7259e 100644
--- a/src/solvers/expand.jl
+++ b/src/solvers/expand.jl
@@ -38,6 +38,17 @@ function expand_cutoff_doctring()
     """
 end
 
+function expand_projection_cutoff_doctring()
+    return """
+    The projection_cutoff is used to control when to expand the
+    basis and defaults to half the precision of the scalar type
+    of the input state, i.e. ~1e-8 for `Float64`. Specifically,
+    the basis on a link is only expanded if the norm of the
+    density matrix projected into the null space is larger than
+    the projection_cutoff.
+    """
+end
+
 function expand_warning_doctring()
     return """
     !!! warning
@@ -56,7 +67,7 @@ function expand_citation_docstring()
 end
 
 """
-    expand(state::MPS, references::Vector{MPS}; alg="orthogonalize", cutoff)
+    expand(state::MPS, references::Vector{MPS}; alg="orthogonalize", cutoff, projection_cutoff)
 
 Given an MPS `state` and a collection of MPS `references`,
 returns an MPS which is equal to `state`
@@ -67,6 +78,8 @@ See [^global_expansion] for more details.
 
 $(expand_cutoff_doctring())
 
+$(expand_projection_cutoff_doctring())
+
 $(expand_warning_doctring())
 
 $(expand_citation_docstring())
@@ -76,6 +89,7 @@ function expand(
         state::MPS,
         references::Vector{MPS};
         cutoff = (√(eps(real(scalartype(state))))),
+        projection_cutoff = (√(eps(real(scalartype(state))))),
     )
     n = length(state)
     state = orthogonalize(state, n)
@@ -97,7 +111,7 @@ function expand(
         # Apply projectorⱼ
         ρⱼ_projected = apply(apply(projectorⱼ, ρⱼ), projectorⱼ)
         expanded_basisⱼ = basisⱼ
-        if norm(ρⱼ_projected) > 10^3 * eps(real(scalartype(state)))
+        if norm(ρⱼ_projected) > projection_cutoff
             # 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
diff --git a/test/base/test_solvers/test_expand.jl b/test/base/test_solvers/test_expand.jl
index 5faf066..0638ee3 100644
--- a/test/base/test_solvers/test_expand.jl
+++ b/test/base/test_solvers/test_expand.jl
@@ -1,7 +1,7 @@
 @eval module $(gensym())
-using ITensors: scalartype
+using ITensors: ITensor, Index, scalartype
 using ITensorMPS:
-    OpSum, MPO, MPS, expand, inner, linkdims, maxlinkdim, random_mps, siteinds, tdvp
+    OpSum, MPO, MPS, apply, expand, inner, linkdims, maxlinkdim, random_mps, siteinds, tdvp
 using ITensorMPS.Experimental: dmrg
 using LinearAlgebra: normalize
 using StableRNGs: StableRNG
@@ -91,5 +91,55 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         # TODO: Use `fourthroot`/`∜` in Julia 1.10 and above.
         @test inner(state', operator, state) ≈ reference_energy rtol = 5 * eps(real(elt))^(1 // 4)
     end
+    if elt == Float64
+        @testset "Regression test issue #160" begin
+            # Regression test for https://github.com/ITensor/ITensorMPS.jl/issues/160
+            s = siteinds("S=1/2", 4)
+            i = Index(2, "Link,l=1")
+            j = Index(3, "Link,l=2")
+            k = Index(2, "Link,l=3")
+            l = Index(4, "HLink,l=1")
+            m = Index(6, "HLink,l=2")
+            n = Index(4, "HLink,l=3")
+            ψ_data = [
+                [-0.8946112605757434 - 3.6983198682267044e-19im 0.12770990068113935 + 0.3385178997037456im;;; -0.09831866026969982 + 9.596903243670168e-18im -0.0858106795611182 - 0.22745653126540122im],
+                [-0.8910784777311369 - 0.15423796632294098im 0.07170978595132114 + 0.37416085631695717im; 0.28118596675443397 + 0.048673375447632186im 0.10671963134029532 + 0.5568576166502048im;;; 0.0574469949727872 + 0.012881670393400376im 0.020337396850339782 + 0.14414434717491714im; -0.5227159163029482 - 0.11721159907883279im 0.030267642036559694 + 0.21452644770567245im;;; -0.05787580381295044 - 0.004553341135313366im -0.02645607106024196 - 0.091233732757172im; -0.4916242418311841 - 0.038679460765149255im -0.039370853144135375 - 0.13578326382500433im],
+                [-0.9076467289406478 + 0.0im 0.13496759758350133 + 0.35775609151581605im; 0.26511388984285783 + 0.0im 0.2313864094611697 + 0.6133310423992547im; 0.3254105727468315 + 0.0im 0.18794511667614308 + 0.49818130260948507im;;; -0.07783131185640604 - 6.718564066141929e-7im -0.05458384118196475 - 0.14468805481202962im; 0.6555262202240988 - 2.9032095377810825e-12im -0.09357940472115472 - 0.24804890035886243im; -0.7511508912806617 - 1.8739618375970262e-6im -0.076006129623389 - 0.20148297391321887im],
+                [-0.9270540711350147 + 0.0im 0.13234125887211426 + 0.3507941567246226im; 0.3749276586116248 + 0.0im 0.3272298001989081 + 0.8673810631261137im;;;],
+            ]
+            ψ = MPS(
+                [
+                    ITensor(ψ_data[1], s[1], i),
+                    ITensor(ψ_data[2], i, s[2], j),
+                    ITensor(ψ_data[3], j, s[3], k),
+                    ITensor(ψ_data[4], k, s[4]),
+                ]
+            )
+            H_data = [
+                [0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 1.0681415022205298 + 0.0im;;;; 4.133529502045017 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im;;;; 0.0 + 0.0im 2.922846741130685 + 0.0im;;; 2.922846741130685 + 0.0im 0.0 + 0.0im;;;; 1.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 1.0 + 0.0im],
+                [1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 1.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0681415022205298 + 0.0im;;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.7071067811865475 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im; 0.7071067811865475 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;;; 0.0 + 0.0im 0.0 + 0.0im; 0.7071067811865475 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.7071067811865475 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 2.922846741130685 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 1.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im],
+                [1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 1.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0681415022205298 + 0.0im;;;; 0.0 + 0.0im 0.0 + 0.0im; -1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im -1.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im -1.0 + 0.0im; -1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; -1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 1.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im],
+                [1.0 + 0.0im 0.0 + 0.0im; -1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im -1.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 1.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im; -1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0681415022205298 + 0.0im;;;;],
+            ]
+            H = MPO(
+                [
+                    ITensor(H_data[1], s[1]', s[1], l),
+                    ITensor(H_data[2], l, s[2]', s[2], m),
+                    ITensor(H_data[3], m, s[3]', s[3], n),
+                    ITensor(H_data[4], n, s[4]', s[4]),
+                ]
+            )
+
+            Hψ = apply(H, ψ)
+            ψ_expanded = expand(ψ, [Hψ]; alg = "orthogonalize")
+            @test inner(ψ_expanded, ψ) ≈ 1.0
+            ψ_expanded2 = expand(ψ, [normalize(Hψ)]; alg = "orthogonalize")
+            @test inner(ψ_expanded2, ψ) ≈ 1
+            ψ_expanded3 = expand(ψ, H; alg = "global_krylov", krylovdim = 1)
+            @test inner(ψ_expanded3, ψ) ≈ 1
+            ψ_expanded4 = expand(ψ, H; alg = "global_krylov", krylovdim = 1, cutoff = 1.0e-6)
+            @test inner(ψ_expanded4, ψ) ≈ 1
+        end
+    end
 end
 end