[ITensorMPS] Bond dimensions produced by randomCircuitMPS

src/mps/mps.jl

@@ -127,7 +127,7 @@ function randomizeMPS!(M::MPS, sites::Vector{<:Index}, linkdim=1)
   end
 end
 
-function randomCircuitMPS(
+function myRandomCircuitMPS(
   ::Type{ElT}, sites::Vector{<:Index}, linkdim::Int; kwargs...
 ) where {ElT<:Number}
   _rmatrix(::Type{Float64}, n, m) = NDTensors.random_orthog(n, m)
@@ -144,15 +144,28 @@ function randomCircuitMPS(
 
   l = Vector{Index}(undef, N)
 
+  # Compute the bond dimension for each link.
+  # It should be the the minimum of the product of local
+  # Hilbert space dimensions taken from either end of the chain
+  # or the `linkdim`.
+  maxdims = Vector{Int}(undef, N - 1)
+  maxdims[1] = min(dim(sites[1]), linkdim)
+  for i in 2:(N - 1)
+    maxdims[i] = min(maxdims[i - 1] * dim(sites[i]), linkdim)
+  end
+  maxdims[N - 1] = min(dim(sites[N]), maxdims[N - 1])
+  for i in (N - 2):-1:1
+    maxdims[i] = min(dim(sites[i + 1]) * maxdims[i + 1], maxdims[i])
+  end
+
   d = dim(sites[N])
-  chi = min(linkdim, d)
+  chi = maxdims[N - 1]
   l[N - 1] = Index(chi, "Link,l=$(N-1)")
   O = _rmatrix(ElT, chi, d)
   M[N] = itensor(O, l[N - 1], sites[N])
 
   for j in (N - 1):-1:2
-    chi *= dim(sites[j])
-    chi = min(linkdim, chi)
+    chi = maxdims[j - 1]
     l[j - 1] = Index(chi, "Link,l=$(j-1)")
     O = _rmatrix(ElT, chi, dim(sites[j]) * dim(l[j]))
     T = reshape(O, (chi, dim(sites[j]), dim(l[j])))

test/mps.jl

@@ -197,6 +197,16 @@ include("util.jl")
     @test norm(phic[4]) ≈ 1.0
   end
 
+  @testset "randomCircuitMPS bond dimensions" begin
+    phi = randomCircuitMPS(ComplexF64, sites, 32)
+    expected_dims = [2, 4, 8, 16, 32, 16, 8, 4, 2]
+
+    for i in 1:9
+      l = getfirst(x -> hastags(x, "Link,l=$(i)"), inds(mps[i]))
+      @test dim(l) == expected_dims[i]
+    end
+  end
+
   @testset "inner different MPS" begin
     phi = randomMPS(sites)
     psi = randomMPS(sites)