Wolfram Language Paclet Repository
Community-contributed installable additions to the Wolfram Language
Matrix Product States
A (MPS) writes the amplitudes of a quantum state on sites as a one-dimensional chain of rank-3 tensors with rank-2 boundaries. The ansatz is exact when the bond dimension is large enough to capture every Schmidt coefficient across every cut, and it is the dominant variational class behind DMRG, TEBD, and time evolution of one-dimensional quantum systems. This tutorial walks through the MPS-specific routines that the paclet exposes at Tier 7: canonical forms, overlaps and norms, Schmidt decomposition, entanglement entropy, and bond truncation. The treatment is practical — every concept is paired with a runnable example built on the same six-site chain so the reader can read the numbers off the page and predict how each operation will reshape the next one.
MatrixProductState
N
In[205]:=
Needs["Wolfram`TensorNetworks`"]
1. What an MPS is
An MPS on sites with physical dimension and bond dimension expresses the coefficients of a state as a contracted chain of site tensors. The bulk tensors carry three indices : a left bond, a right bond, and a physical leg. The two boundary tensors have only one bond and one physical leg each. The total number of parameters in the ansatz is — linear in the chain length, instead of the parameters of the dense wavefunction. The bond dimension controls expressiveness: every Schmidt rank across every bipartition of the chain is upper bounded by , and the ansatz is exact when . For one-dimensional gapped Hamiltonians, however, the Schmidt rank saturates at a constant set by the desired truncation accuracy , regardless of system size — this is the area law that makes MPS competitive for ground-state simulation.
N
d
D
ψ=,…,⋯|⋯〉
∑
σ
1
σ
N
[1]
A
[2]
A
[N]
A
σ
1
σ
N
[i]
A
(,,)
α
i
β
i
σ
i
2·D·d+(N-2)··d
2
D
N
d
D
D
D≥
⌊N/2⌋
d
D(ε)
ε
The paclet builds a random MPS via with the named pattern . The following constructs a six-site chain with bond dimension four and qubit physical legs.
"MPS"[N,D,d]
In[206]:=
BlockRandomSeedRandom[42];mps=
["MPS"[6,4,2]]
 | RandomTensorNetwork |
Out[206]=
TensorNetwork
 |  |
|
The MPS routines in this paclet assume the open-boundary tensor layout used by : every bulk tensor is shaped and each boundary tensor is shaped . Reading the tensor dimensions confirms this for the chain just built.
{leftBond,rightBond,physical}
{bond,physical}
In[207]:=
 | TensorNetworkData |
 | SparseTensorNetwork |
Out[207]=
TensorsSparseArray
,SparseArray
,SparseArray
,SparseArray
,SparseArray
,SparseArray
,Dimensions{{4,2},{4,4,2},{4,4,2},{4,4,2},{4,4,2},{4,2}},Hyperedges{{1,7},{1,2,8},{2,3,9},{3,4,10},{4,5,11},{5,12}},Indices{{,},{,,},{,,},{,,},{,,},{,}},Vertices{1,2,3,4,5,6},Bonds{{,}4,{,}4,{,}4,{,}4,{,}4},Contractions{{{,},},{{,},{,},},{{,},{,},},{{,},{,},},{{,},{,},},{{,},}},ContractionIndices{{1,7},{1,2,8},{2,3,9},{3,4,10},{4,5,11},{5,12}},FreeIndices{7,8,9,10,11,12}
|  |
|
|  |
|
| |
|
| |
|
| |
|
|  |
|
1
1
7
1
1
2
2
2
8
2
2
3
3
3
9
3
3
4
4
4
10
4
4
5
5
5
11
5
5
6
12
6
1
1
1
2
2
2
2
3
3
3
3
4
4
4
4
5
5
5
5
6
1
1
1
2
7
1
1
1
1
2
2
2
2
3
8
2
2
2
2
3
3
3
3
4
9
3
3
3
3
4
4
4
4
5
10
4
4
4
4
5
5
5
5
6
11
5
5
5
5
6
12
6
In[208]:=
mps["Hypergraph"]
Out[208]=
The chain has two rank-2 boundary tensors at the ends and four rank-3 bulk tensors in the middle, with bond indices linking neighbors and one free physical leg per site. The bond indices form a hyperedge between two adjacent tensors, the physical legs are free legs that are not contracted. Reading the network's property confirms this: the first six entries on the integer index list are bond indices linking sites 1-2, 2-3 and so on, and the next six are the physical indices unique to each site. A tensor-ring (periodic) topology is requested with ; that variant promotes the boundary tensors to rank-3 by giving the first and last sites an extra bond index that links them, closing the chain into a ring. Tensor rings are useful for thermodynamic-limit calculations and translationally invariant ansatze, but the algorithms documented in the remainder of this tutorial — canonicalization, overlap, Schmidt decomposition and truncation — target the open-boundary layout. For periodic chains, contract the ring directly or transform to an open MPS with a boundary virtual leg before invoking these primitives.
"Hyperedges"
"Boundary""Periodic"
2. Canonical forms
An MPS has a gauge freedom: inserting any invertible matrix and its inverse on a bond does not change the contracted state but does change the individual site tensors. This redundancy is not a bug; it is the structural reason algorithms like DMRG can converge quickly and stably. Canonical forms fix the gauge in a controlled way, putting the chain into a state where one site holds the wavefunction's spectral content and every other site is an isometry. Algorithms then operate on the small spectral block instead of the whole chain. Different gauges of an MPS represent the same physical state but expose different structure for downstream algorithms. The function produces three of them. The form makes every site tensor left-isometric, meaning the matrix obtained by grouping its left bond and physical leg satisfies .A=I. The form is the mirror image: every site is right-isometric, . The mixed form is left-canonical for sites strictly below and right-canonical for sites strictly above ; site is the orthogonality center carrying the wavefunction's spectral weight and is the local site updated by DMRG. Internally each form is produced by a sweep of singular value decompositions: a left sweep absorbs the residual factor into the next tensor and keeps as the new left-isometric site, and the right sweep is the transpose of this procedure. The orthogonality center after one full sweep sits at the last site, which is why is equivalent to .
"Left"
†
A
"Right"
A.=I
†
A
{"Mixed",k}
k
k
k
S·
†
V
U
"Left"
{"Mixed",N}
In[209]:=
canon=
[mps,"Left"]
| MPSCanonicalForm |
Out[209]=
TensorNetwork
| |
|
In[210]:=
| TensorNetworkData |
| SparseTensorNetwork |
Out[210]=
TensorsSparseArray
,SparseArray
,SparseArray
,SparseArray
,SparseArray
,SparseArray
,Dimensions{{2,2},{2,4,2},{4,4,2},{4,4,2},{4,4,2},{4,2}},Hyperedges{{1,7},{1,2,8},{2,3,9},{3,4,10},{4,5,11},{5,12}},Indices{{,},{,,},{,,},{,,},{,,},{,}},Vertices{1,2,3,4,5,6},Bonds{{,}2,{,}4,{,}4,{,}4,{,}4},Contractions{{{,},},{{,},{,},},{{,},{,},},{{,},{,},},{{,},{,},},{{,},}},ContractionIndices{{1,7},{1,2,8},{2,3,9},{3,4,10},{4,5,11},{5,12}},FreeIndices{7,8,9,10,11,12}
|  |
|
| |
|
| |
|
| |
|
| |
|
|  |
|
1
1
7
1
1
2
2
2
8
2
2
3
3
3
9
3
3
4
4
4
10
4
4
5
5
5
11
5
5
6
12
6
1
1
1
2
2
2
2
3
3
3
3
4
4
4
4
5
5
5
5
6
1
1
1
2
7
1
1
1
1
2
2
2
2
3
8
2
2
2
2
3
3
3
3
4
9
3
3
3
3
4
4
4
4
5
10
4
4
4
4
5
5
5
5
6
11
5
5
5
5
6
12
6
In[211]:=
| MPSCanonicalQ |
Out[211]=
True
The companion predicate verifies the isometry condition site by site, with a default tolerance of on the Frobenius distance from the identity. The same predicate accepts the right and mixed forms.
-10
10
In[212]:=
[mps,"Right"],"Right",
[mps,{"Mixed",3}],{"Mixed",3}
 | MPSCanonicalQ |
 | MPSCanonicalForm |
| MPSCanonicalQ |
| MPSCanonicalForm |
Out[212]=
{True,True}
Both options and feed into the singular value decomposition that performs each canonicalization step. Asking for a left-canonical form with a smaller bond truncates as it canonicalizes, dropping the smallest Schmidt values at every cut.
"MaxBond"D
"Tolerance"ϵ
In[213]:=
MapDimensions,
[mps,"Left","MaxBond"2]["Tensors"]
| MPSCanonicalForm |
Out[213]=
{{2,2},{2,2,2},{2,2,2},{2,2,2},{2,2,2},{2,2}}
Every bond has been reduced from four down to two. The truncation happens SVD by SVD as the sweep proceeds: at each bond the SVD returns Schmidt values in descending order, the largest are kept, and the remainder is dropped together with the small singular modes of its product with the next site. The discarded weight contributes a local error and the accumulated error after the full sweep is bounded by the sum of the discarded squared Schmidt values across all bonds, the standard truncation-error estimate. The same truncation can also be reached by passing a positive : singular values below the cutoff are dropped, which is the variant to use when you want a relative accuracy target rather than a hard bond cap. The two cutoffs compose — if both are set, the algorithm keeps at most singular values and additionally drops any below , whichever is more aggressive at the given bond.
D
max
"Tolerance"
"MaxBond"
"Tolerance"
3. Overlap and norm
Overlaps are the workhorse subroutine of nearly every MPS algorithm: expectation values of local operators, energies in DMRG, fidelities in compression, and infidelity estimates in TEBD all reduce to inner products of two MPS. The inner product between two MPS is computed by . The second argument is conjugated, matching the bra convention. The implementation contracts the two chains side-by-side, sweeping a transfer-matrix from left to right; the cost grows linearly in the chain length and as in the bond dimension. The transfer matrix at site is the partial contraction of the chains over the physical leg and the left bond, leaving the right bond open. The sweep multiplies these transfer matrices in sequence and the final scalar is the overlap. Self-overlap returns the squared norm.
〈|〉
ψ
1
ψ
2
3
D
i
In[214]:=
[mps,mps],
[mps]
| MPSOverlap |
| MPSNorm |
Out[214]=
{69.3437,8.32729}
The two numbers satisfy =〈ψ|ψ〉, because is defined as . To rescale a state to unit norm without rebuilding it from scratch use , which divides each site tensor by so that the global norm is one.
2
ψ
Sqrt[Abs[MPSOverlap[mps,mps]]]
1/N
ψ
In[215]:=
nmps=
[mps];
[nmps]
| MPSNormalize |
| MPSNorm |
Out[216]=
1.
Normalization is effectively idempotent: a second application changes nothing up to floating-point error. Distributing the rescaling factor across every site instead of dumping it on the first tensor keeps the matrix elements of each site small, which avoids overflow in long chains and keeps later overlaps numerically well-conditioned. Normalizing before computing Schmidt values or entanglement entropy keeps the probabilities on the unit simplex without per-call rescaling. Note that the overlap implementation never expands the chains into their dense -element vectors — for the six-site qubit chain above the dense state has 64 amplitudes, but the transfer-matrix sweep evaluates a chain of contractions whose intermediate sizes are at most = 16. This is the reason MPS arithmetic stays tractable when is large.
p
i
N
d
2
D
N
4. Schmidt decomposition and entanglement entropy
Splitting the chain across a bond at site realizes a Schmidt decomposition: there is a basis of left states on sites and right states on sites such that with the Schmidt values non-negative and the partial probabilities = summing to one. The function returns the normalized list.
k
|〉
L
i
1…k
|〉
R
i
k+1…N
ψ=|〉|〉
r
∑
i=1
λ
i
L
i
R
i
λ
i
p
i
2
λ
i
In[217]:=
svals=
[nmps,3];{svals,Total[svals^2]}
| MPSSchmidtValues |
Out[218]=
{{0.920345,0.305964,0.21164,0.120663},1.}
Internally the function builds the mixed-canonical form centered at site , reshapes the center tensor into the bipartition matrix, and reads the singular values directly. Because the Schmidt values are normalized, the list is exactly the bond's reduced density spectrum . The decay rate of these values across is the fingerprint of the state's entanglement structure: an exponential decay indicates a state easily approximated by a small bond dimension, while a slow power-law decay is the signature of critical (gapless) physics that requires a larger for the same accuracy.
k
p
i
i
D
The von Neumann entanglement entropy of the same bipartition is returned by . For a product state at bond it is zero; for a maximally entangled cut with bond dimension it is .
S=-log
∑
i
p
i
p
i
k
D
logD
In[219]:=
| MPSEntanglementEntropy |
Out[219]=
0.563041
The entanglement is below the ≈ 1.386 ceiling because the random tensors did not saturate every Schmidt direction. Plotting the entropy across all interior bonds gives the entanglement profile of the state, often called the entanglement spectrum's first moment.
log(4)
In[220]:=
Table
[nmps,k],{k,1,5}
| MPSEntanglementEntropy |
Out[220]=
{0.201076,0.477574,0.563041,0.657778,0.398148}
The arch shape is typical of a finite chain: the boundary bonds carry less entanglement than the bulk because there are fewer sites on the short side of the cut, so the Schmidt rank is constrained more by the small block than by the bond dimension. For ground states of one-dimensional gapped Hamiltonians, the area law caps the bulk plateau by a constant, regardless of system size, which is the structural reason MPS work — a fixed bond dimension is enough to capture the entanglement of arbitrarily long chains. Critical (gapless) one-dimensional systems instead show , the celebrated logarithmic CFT scaling, where is the central charge of the universality class and is the conformal block length of the cut. A diagnostic plot of the entanglement profile against on a converged DMRG state is therefore a quick way to detect criticality without computing dynamical structure factors.
S∼logL(k)
c
6
c
L(k)
logk
5. Bond truncation and compression
Variational compression of an MPS — reducing every bond down to a target dimension while losing as little fidelity as possible — is exposed as . The standard recipe is a left sweep that performs an SVD at every bond, keeps the top singular values, and absorbs the remainder into the next tensor; the result is a left-canonical MPS by construction, optionally renormalized.
D
max
D
max
In[221]:=
truncated=
[nmps,2];Map[Dimensions,truncated["Tensors"]]
| MPSTruncate |
Out[222]=
{{2,2},{2,2,2},{2,2,2},{2,2,2},{2,2,2},{2,2}}