Skip to content

Commit c01cba6

Browse files
mtfishmanmtfishman
authored
Looser bounds in expand (#164)
1 parent cc227c2 commit c01cba6

File tree

2 files changed

+68
-4
lines changed

2 files changed

+68
-4
lines changed

‎src/solvers/expand.jl‎

Lines changed: 16 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,17 @@ function expand_cutoff_doctring()
3838
"""
3939
end
4040

41+
function expand_projection_cutoff_doctring()
42+
return """
43+
The projection_cutoff is used to control when to expand the
44+
basis and defaults to half the precision of the scalar type
45+
of the input state, i.e. ~1e-8 for `Float64`. Specifically,
46+
the basis on a link is only expanded if the norm of the
47+
density matrix projected into the null space is larger than
48+
the projection_cutoff.
49+
"""
50+
end
51+
4152
function expand_warning_doctring()
4253
return """
4354
!!! warning
@@ -56,7 +67,7 @@ function expand_citation_docstring()
5667
end
5768

5869
"""
59-
expand(state::MPS, references::Vector{MPS}; alg="orthogonalize", cutoff)
70+
expand(state::MPS, references::Vector{MPS}; alg="orthogonalize", cutoff, projection_cutoff)
6071
6172
Given an MPS `state` and a collection of MPS `references`,
6273
returns an MPS which is equal to `state`
@@ -67,6 +78,8 @@ See [^global_expansion] for more details.
6778
6879
$(expand_cutoff_doctring())
6980
81+
$(expand_projection_cutoff_doctring())
82+
7083
$(expand_warning_doctring())
7184
7285
$(expand_citation_docstring())
@@ -76,6 +89,7 @@ function expand(
7689
state::MPS,
7790
references::Vector{MPS};
7891
cutoff = (√(eps(real(scalartype(state))))),
92+
projection_cutoff = (√(eps(real(scalartype(state))))),
7993
)
8094
n = length(state)
8195
state = orthogonalize(state, n)
@@ -97,7 +111,7 @@ function expand(
97111
# Apply projectorⱼ
98112
ρⱼ_projected = apply(apply(projectorⱼ, ρⱼ), projectorⱼ)
99113
expanded_basisⱼ = basisⱼ
100-
if norm(ρⱼ_projected) > 10^3 * eps(real(scalartype(state)))
114+
if norm(ρⱼ_projected) > projection_cutoff
101115
# Diagonalize projected density matrix ρⱼ_projected
102116
# to compute reference_basisⱼ, which spans part of right basis
103117
# of references which is orthogonal to right basis of state

‎test/base/test_solvers/test_expand.jl‎

Lines changed: 52 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
@eval module $(gensym())
2-
using ITensors: scalartype
2+
using ITensors: ITensor, Index, scalartype
33
using ITensorMPS:
4-
OpSum, MPO, MPS, expand, inner, linkdims, maxlinkdim, random_mps, siteinds, tdvp
4+
OpSum, MPO, MPS, apply, expand, inner, linkdims, maxlinkdim, random_mps, siteinds, tdvp
55
using ITensorMPS.Experimental: dmrg
66
using LinearAlgebra: normalize
77
using StableRNGs: StableRNG
@@ -91,5 +91,55 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
9191
# TODO: Use `fourthroot`/`∜` in Julia 1.10 and above.
9292
@test inner(state', operator, state) ≈ reference_energy rtol = 5 * eps(real(elt))^(1 // 4)
9393
end
94+
if elt == Float64
95+
@testset "Regression test issue #160" begin
96+
# Regression test for https://github.com/ITensor/ITensorMPS.jl/issues/160
97+
s = siteinds("S=1/2", 4)
98+
i = Index(2, "Link,l=1")
99+
j = Index(3, "Link,l=2")
100+
k = Index(2, "Link,l=3")
101+
l = Index(4, "HLink,l=1")
102+
m = Index(6, "HLink,l=2")
103+
n = Index(4, "HLink,l=3")
104+
ψ_data = [
105+
[-0.8946112605757434 - 3.6983198682267044e-19im 0.12770990068113935 + 0.3385178997037456im;;; -0.09831866026969982 + 9.596903243670168e-18im -0.0858106795611182 - 0.22745653126540122im],
106+
[-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],
107+
[-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],
108+
[-0.9270540711350147 + 0.0im 0.13234125887211426 + 0.3507941567246226im; 0.3749276586116248 + 0.0im 0.3272298001989081 + 0.8673810631261137im;;;],
109+
]
110+
ψ = MPS(
111+
[
112+
ITensor(ψ_data[1], s[1], i),
113+
ITensor(ψ_data[2], i, s[2], j),
114+
ITensor(ψ_data[3], j, s[3], k),
115+
ITensor(ψ_data[4], k, s[4]),
116+
]
117+
)
118+
H_data = [
119+
[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],
120+
[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],
121+
[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],
122+
[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;;;;],
123+
]
124+
H = MPO(
125+
[
126+
ITensor(H_data[1], s[1]', s[1], l),
127+
ITensor(H_data[2], l, s[2]', s[2], m),
128+
ITensor(H_data[3], m, s[3]', s[3], n),
129+
ITensor(H_data[4], n, s[4]', s[4]),
130+
]
131+
)
132+
133+
Hψ = apply(H, ψ)
134+
ψ_expanded = expand(ψ, [Hψ]; alg = "orthogonalize")
135+
@test inner(ψ_expanded, ψ) ≈ 1.0
136+
ψ_expanded2 = expand(ψ, [normalize(Hψ)]; alg = "orthogonalize")
137+
@test inner(ψ_expanded2, ψ) ≈ 1
138+
ψ_expanded3 = expand(ψ, H; alg = "global_krylov", krylovdim = 1)
139+
@test inner(ψ_expanded3, ψ) ≈ 1
140+
ψ_expanded4 = expand(ψ, H; alg = "global_krylov", krylovdim = 1, cutoff = 1.0e-6)
141+
@test inner(ψ_expanded4, ψ) ≈ 1
142+
end
143+
end
94144
end
95145
end

0 commit comments

Comments
 (0)