diff --git a/Project.toml b/Project.toml
index edb849908c..b15ffe42ee 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.3.18"
 [deps]
 BitIntegers = "c3b6d118-76ef-56ca-8cc7-ebb389d030a1"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
diff --git a/src/mps/mps.jl b/src/mps/mps.jl
index 3b9120a02a..8b990293b4 100644
--- a/src/mps/mps.jl
+++ b/src/mps/mps.jl
@@ -148,15 +148,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 = NDTensors.random_unitary(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 = NDTensors.random_unitary(ElT, chi, dim(sites[j]) * dim(l[j]))
     T = reshape(O, (chi, dim(sites[j]), dim(l[j])))
diff --git a/test/mps.jl b/test/mps.jl
index 0fb41baaf0..e88751db10 100644
--- a/test/mps.jl
+++ b/test/mps.jl
@@ -197,6 +197,16 @@ include("util.jl")
     @test norm(phic[4]) ≈ 1.0
   end
 
+  @testset "randomMPS bond dimensions" begin
+    phi = randomMPS(ComplexF64, sites, 50)
+    expected_dims = [2, 4, 8, 16, 32, 16, 8, 4, 2]
+
+    for i in 1:9
+      l = linkind(phi, i)
+      @test dim(l) == expected_dims[i]
+    end
+  end
+  
   @testset "randomMPS with chi>1" for linkdims in [1, 4]
     phi = randomMPS(Float32, sites; linkdims)
     @test LinearAlgebra.promote_leaf_eltypes(phi) === Float32
@@ -482,7 +492,8 @@ include("util.jl")
 
     ϕ2 = +(ψ1, ψ2; alg="directsum")
     for j in 1:8
-      @test linkdim(ϕ2, j) == χ1 + χ2
+      #@test linkdim(ϕ2, j) == χ1 + χ2
+      @test linkdim(ϕ2, j) == linkdim(ψ1, j) + linkdim(ψ2, j)
     end
     @test inner(ϕ2, ψ1) + inner(ϕ2, ψ2) ≈ inner(ϕ2, ϕ2)
   end