MPS and MPO Basics

This tutorial continues from ITensor Basics. This tutorial will walk you through how to define MPS, MPO, and density matrices in the Hilbert space, and how to compute expectation values and reduced density matrices.

For the original ITensorMPS tutorial on MPS and MPO objects, see the MPS and MPO examples page.

For MPS/MPO definitions and contraction pictures, see Tensor Networks in Physics.

A many-body wavefunction on N sites can be written as

\[|\psi\rangle = \sum_{\sigma_1,\ldots,\sigma_N} c_{\sigma_1\cdots\sigma_N} |\sigma_1,\ldots,\sigma_N\rangle.\]

A matrix product state represents the coefficient tensor

\[c_{\sigma_1\cdots\sigma_N}\]

as a product of smaller tensors:

\[c_{\sigma_1\cdots\sigma_N} = A^{\sigma_1}_1 A^{\sigma_2}_2 \cdots A^{\sigma_N}_N.\]

This is useful because many physical states can be represented accurately using moderate bond dimensions.

MPS from Tensors

Setup

We load ITensors for basic tensor operations and import ITensorMPS only to inspect the object stored inside the ProcessTensors.jl wrappers.

We use ProcessTensors.jl for the public MPS/MPO interface used throughout the rest of these tutorials.

using ITensors
import ITensorMPS
import LinearAlgebra
using ProcessTensors

roundreal(x; digits=6) = round(real(x); digits=digits)

roundreal (generic function with 1 method)

Decomposing a dense tensor into an MPS

Let us first build a dense four-site wavefunction. Each site has dimension 2, so the full coefficient tensor has size 2 × 2 × 2 × 2. We fill it with the numbers 1, 2, ..., 16, then convert it into an MPS.

mps_sites = [Index(2, "Qubit,Site,n=$n") for n in 1:4]

coeffs = reshape(ComplexF64.(1:16), 2, 2, 2, 2)

ψ_dense = MPS(coeffs, mps_sites)

println(ψ_dense)

4-element MPS{Hilbert}
  site dims: 2, 2, 2, 2
  link dims: 2, 4, 2
  maxlinkdim: 4
  combiners: none

  tensors:
    [1] ((dim=2|id=610|"Qubit,Site,n=1"), (dim=2|id=699|"Link,l=1"))
    [2] ((dim=2|id=699|"Link,l=1"), (dim=2|id=507|"Qubit,Site,n=2"), (dim=4|id=830|"Link,l=2"))
    [3] ((dim=4|id=830|"Link,l=2"), (dim=2|id=520|"Qubit,Site,n=3"), (dim=2|id=0|"Link,l=3"))
    [4] ((dim=2|id=0|"Link,l=3"), (dim=2|id=881|"Qubit,Site,n=4"))

Each ITensor core in this MPS contains physical indices (the mps_sites), and link indices connecting the cores. One can access these indices from the MPS object using the siteinds and linkinds functions.

The MPS stores the same state as the dense tensor, but in factorized form. To check this, we contract all MPS tensors back into one full tensor and compare the coefficients.

function dense_array_from_mps(ψ)
    full_tensor = ITensor(1.0)
    for n in eachindex(ψ)
        full_tensor *= ψ[n]
    end
    return Array(full_tensor, siteinds(ψ)...)
end

coeffs_from_mps = dense_array_from_mps(ψ_dense)

basis_index = (1, 2, 2, 1)

println("Dense coefficient at $basis_index: ", coeffs[basis_index...])
println("MPS coefficient at $basis_index:   ", coeffs_from_mps[basis_index...])

@assert isapprox(coeffs[basis_index...], coeffs_from_mps[basis_index...]; atol=1e-10)

Dense coefficient at (1, 2, 2, 1): 7.0 + 0.0im
MPS coefficient at (1, 2, 2, 1):   6.9999999999999964 + 0.0im

Physical site indices with siteinds

For the sake of illustration, from now on, let us focus on a physical system and define everything in the context of spin-1/2 chains. The first ingredient we need to define a spin system wavefunction are the physical spin indices. siteinds("S=1/2", N) creates the physical Hilbert-space indices for a spin chain of length N.

N = 6
sites = siteinds("S=1/2", N)

println("Spin-chain site indices:")
foreach(println, sites)

println("First spin site:")
println("  index = ", sites[1])
println("  dim   = ", dim(sites[1]))
println("  tags  = ", tags(sites[1]))

@assert length(sites) == N
@assert all(s -> dim(s) == 2, sites)

Spin-chain site indices:
(dim=2|id=691|"S=1/2,Site,n=1")
(dim=2|id=882|"S=1/2,Site,n=2")
(dim=2|id=772|"S=1/2,Site,n=3")
(dim=2|id=716|"S=1/2,Site,n=4")
(dim=2|id=698|"S=1/2,Site,n=5")
(dim=2|id=391|"S=1/2,Site,n=6")
First spin site:
  index = (dim=2|id=691|"S=1/2,Site,n=1")
  dim   = 2
  tags  = "S=1/2,Site,n=1"

A product-state MPS

Consider the Néel state

\[|\psi_{\mathrm{N\acute eel}}\rangle = |\uparrow\downarrow\uparrow\downarrow\uparrow\downarrow\rangle.\]

In ITensorMPS syntax, this is very direct.

neel_config = ["Up", "Dn", "Up", "Dn", "Up", "Dn"]

ψ_neel = MPS(sites, neel_config)

println(ψ_neel)

6-element MPS{Hilbert}
  site dims: 2, 2, 2, 2, 2, 2
  link dims: 1, 1, 1, 1, 1
  maxlinkdim: 1
  combiners: none

  tensors:
    [1] ((dim=2|id=691|"S=1/2,Site,n=1"), (dim=1|id=527|"Link,l=1"))
    [2] ((dim=1|id=527|"Link,l=1"), (dim=2|id=882|"S=1/2,Site,n=2"), (dim=1|id=246|"Link,l=2"))
    [3] ((dim=1|id=246|"Link,l=2"), (dim=2|id=772|"S=1/2,Site,n=3"), (dim=1|id=453|"Link,l=3"))
    [4] ((dim=1|id=453|"Link,l=3"), (dim=2|id=716|"S=1/2,Site,n=4"), (dim=1|id=875|"Link,l=4"))
    [5] ((dim=1|id=875|"Link,l=4"), (dim=2|id=698|"S=1/2,Site,n=5"), (dim=1|id=461|"Link,l=5"))
    [6] ((dim=1|id=461|"Link,l=5"), (dim=2|id=391|"S=1/2,Site,n=6"))

ProcessTensors.jl prints this as an MPS{Hilbert}.

This means:

• the object is an ordinary Hilbert-space MPS,
• the underlying ITensorMPS object is stored in .core,
• common MPS functions are forwarded to .core.

println("Wrapper type: ", typeof(ψ_neel))
println("Core type:    ", typeof(ψ_neel.core))

@assert ψ_neel isa MPS{Hilbert}
@assert ψ_neel.core isa ITensorMPS.MPS
@assert maxlinkdim(ψ_neel) == 1

Wrapper type: MPS{Hilbert, Nothing}
Core type:    ITensorMPS.MPS

Let us inspect one local MPS tensor and see if it contains the expected site and link indices.

middle = 3
println("Indices of ψ_neel[$middle]:")
println(inds(ψ_neel[middle]))

@assert siteind(ψ_neel, middle) in collect(inds(ψ_neel[middle]))

Indices of ψ_neel[3]:
((dim=1|id=246|"Link,l=2"), (dim=2|id=772|"S=1/2,Site,n=3"), (dim=1|id=453|"Link,l=3"))

Entanglement entropy of a pure MPS

For a pure state $|\psi\rangle,$ the entanglement entropy across a bipartition A|B is

\[S_A = -\operatorname{Tr}(\rho_A \log \rho_A), \qquad \rho_A = \operatorname{Tr}_B |\psi\rangle\langle\psi|.\]

In MPS language, this entropy is computed from the singular values across a bond.

For the product Néel state, the entanglement entropy across every cut should be zero.

S_neel = [entropy(ψ_neel, b) for b in 1:(N - 1)]

println("Bond entropies of the Néel state:")
println(roundreal.(S_neel))

@assert maximum(abs.(S_neel)) < 1e-12

Bond entropies of the Néel state:
[0.0, 0.0, 0.0, 0.0, 0.0]

Let us also build a small GHZ state,

\[|\mathrm{GHZ}\rangle = \frac{1}{\sqrt{2}} \left( |\uparrow\uparrow\cdots\uparrow\rangle + |\downarrow\downarrow\cdots\downarrow\rangle \right).\]

Every bipartition of this state has entropy log(2).

function ghz_state(sites)
    N = length(sites)
    all_ups = fill("Up", N)
    all_dns = fill("Dn", N)
    ghz = (MPS(sites, all_ups) + MPS(sites, all_dns)) / sqrt(2)
    return ghz
end

ψ_ghz = ghz_state(sites)

S_ghz = [entropy(ψ_ghz, b) for b in 1:(N - 1)]

println("Bond entropies of the GHZ state:")
println(roundreal.(S_ghz))

@assert all(S -> isapprox(S, log(2); atol=1e-10), S_ghz)

Bond entropies of the GHZ state:
[0.693147, 0.693147, 0.693147, 0.693147, 0.693147]

MPO and expectation values

Local observables

The local magnetization is

\[\langle S^z_j \rangle = \langle \psi | S^z_j | \psi \rangle.\]

For the Néel state, the answer should alternate between +1/2 and -1/2. The function expect performs the tensor-network contractions for all local one-site observables.

mz = expect(ψ_neel, "Sz")
mx = expect(ψ_neel, "Sx")

println("⟨Sz_j⟩:")
println(roundreal.(mz))

println("⟨Sx_j⟩:")
println(roundreal.(mx))

@assert isapprox.(real.(mz), [0.5, -0.5, 0.5, -0.5, 0.5, -0.5]; atol=1e-12) |> all
@assert all(abs.(mx) .< 1e-12)

⟨Sz_j⟩:
[0.5, -0.5, 0.5, -0.5, 0.5, -0.5]
⟨Sx_j⟩:
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

The same expectation value can also be obtained by building the local operator explicitly as an MPO. For example, the one-site operator

\[O = S^z_3\]

can be written as an OpSum, converted to an MPO, and contracted with the state.

function local_observable_mpo(opname::String, n::Integer, sites)
    os = OpSum()
    os += 1.0, opname, n
    return MPO(os, sites)
end

Sz3_mpo = local_observable_mpo("Sz", 3, sites)

sz3_from_mpo = inner(ψ_neel', Sz3_mpo, ψ_neel)

println("⟨Sz₃⟩ from local MPO: ", roundreal(sz3_from_mpo))
println("⟨Sz₃⟩ from expect:    ", roundreal(mz[3]))

@assert isapprox(sz3_from_mpo, mz[3]; atol=1e-12)

⟨Sz₃⟩ from local MPO: 0.5
⟨Sz₃⟩ from expect:    0.5

Equivalently, we can apply the MPO to the state first:

\[|\phi\rangle = S^z_3 |\psi\rangle,\]

then compute

\[\langle \psi | \phi \rangle.\]

ϕ = apply(Sz3_mpo, ψ_neel; cutoff=0.0)

sz3_from_apply = inner(ψ_neel, ϕ)

println("⟨ψ|Sz₃|ψ⟩ using apply: ", roundreal(sz3_from_apply))

@assert isapprox(sz3_from_apply, mz[3]; atol=1e-12)

⟨ψ|Sz₃|ψ⟩ using apply: 0.5

@assert ϕ isa MPS{Hilbert}

Building Hamiltonians with OpSum

Now consider a simple spin-chain Hamiltonian:

\[H = h_x \sum_{j=1}^{N} S^x_j + J_z \sum_{j=1}^{N-1} S^z_j S^z_{j+1}.\]

The OpSum object lets us write this in a physics-friendly way.

function spin_chain_opsum(N; hx=0.7, Jz=1.0)
    os = OpSum()
    for n in 1:N
        os += hx, "Sx", n
    end
    for n in 1:(N - 1)
        os += Jz, "Sz", n, "Sz", n + 1
    end
    return os
end

hx = 0.7
Jz = 1.0

H_os = spin_chain_opsum(N; hx, Jz)

println("Hamiltonian OpSum:")
println(H_os)

Hamiltonian OpSum:
sum(
  0.7 Sx(1,)
  0.7 Sx(2,)
  0.7 Sx(3,)
  0.7 Sx(4,)
  0.7 Sx(5,)
  0.7 Sx(6,)
  1.0 Sz(1,) Sz(2,)
  1.0 Sz(2,) Sz(3,)
  1.0 Sz(3,) Sz(4,)
  1.0 Sz(4,) Sz(5,)
  1.0 Sz(5,) Sz(6,)
)

Matrix product operators

A matrix product operator represents an operator as a tensor network:

\[\hat{O} = \sum_{\boldsymbol{\sigma},\boldsymbol{\sigma}'} W^{\sigma_1,\sigma'_1}_1 W^{\sigma_2,\sigma'_2}_2 \cdots W^{\sigma_N,\sigma'_N}_N |\boldsymbol{\sigma}\rangle \langle \boldsymbol{\sigma}'|.\]

Constructing MPOs by hand can be difficult because the operator may need nontrivial internal bond dimensions.

ITensorMPS has an efficient algorithm for constructing an MPO from an OpSum.

H_mpo = MPO(H_os, sites)

println(H_mpo)

println("MPO wrapper type: ", typeof(H_mpo))
println("MPO core type:    ", typeof(H_mpo.core))
println("MPO link dimensions: ", collect(linkdims(H_mpo)))

@assert H_mpo isa MPO{Hilbert}
@assert H_mpo.core isa ITensorMPS.MPO

6-element MPO{Hilbert}
  site dims: 2, 2, 2, 2, 2, 2
  link dims: 3, 3, 3, 3, 3
  maxlinkdim: 3
  combiners: none

  tensors:
    [1] ((dim=3|id=885|"Link,l=1"), (dim=2|id=691|"S=1/2,Site,n=1")', (dim=2|id=691|"S=1/2,Site,n=1"))
    [2] ((dim=3|id=885|"Link,l=1"), (dim=3|id=349|"Link,l=2"), (dim=2|id=882|"S=1/2,Site,n=2")', (dim=2|id=882|"S=1/2,Site,n=2"))
    [3] ((dim=3|id=349|"Link,l=2"), (dim=3|id=583|"Link,l=3"), (dim=2|id=772|"S=1/2,Site,n=3")', (dim=2|id=772|"S=1/2,Site,n=3"))
    [4] ((dim=3|id=583|"Link,l=3"), (dim=3|id=136|"Link,l=4"), (dim=2|id=716|"S=1/2,Site,n=4")', (dim=2|id=716|"S=1/2,Site,n=4"))
    [5] ((dim=3|id=136|"Link,l=4"), (dim=3|id=202|"Link,l=5"), (dim=2|id=698|"S=1/2,Site,n=5")', (dim=2|id=698|"S=1/2,Site,n=5"))
    [6] ((dim=3|id=202|"Link,l=5"), (dim=2|id=391|"S=1/2,Site,n=6")', (dim=2|id=391|"S=1/2,Site,n=6"))

MPO wrapper type: MPO{Hilbert, Nothing}
MPO core type:    ITensorMPS.MPO
MPO link dimensions: [3, 3, 3, 3, 3]

Energy expectation value

The energy of a state is

\[E = \langle \psi | H | \psi \rangle.\]

In code:

E_neel = inner(ψ_neel', H_mpo, ψ_neel)

println("Energy of the Néel state: ", roundreal(E_neel))

Energy of the Néel state: -1.25

We can predict this value analytically.

In the Néel state,

\[\langle S^x_j \rangle = 0\]

and each nearest-neighbor pair is anti-aligned:

\[\langle S^z_j S^z_{j+1} \rangle = -\frac{1}{4}.\]

Therefore

\[E = J_z (N-1)\left(-\frac{1}{4}\right).\]

E_expected = -Jz * (N - 1) / 4

println("Expected energy: ", E_expected)

@assert isapprox(real(E_neel), E_expected; atol=1e-12)

Expected energy: -1.25

Density matrices, reduced density matrices and entropy

We now move from pure states to density matrices. A pure state has density matrix

\[\rho = |\psi\rangle\langle\psi|.\]

A classical mixture of two pure states is

\[\rho = p |\psi_1\rangle\langle\psi_1| + (1-p) |\psi_2\rangle\langle\psi_2|.\]

ProcessTensors.jl provides to_dm for exactly this purpose.

A pure density MPO

ρ_neel = to_dm(ψ_neel)

println(ρ_neel)
println("Tr(ρ_neel) = ", roundreal(tr(ρ_neel)))

@assert ρ_neel isa MPO{Hilbert}
@assert isapprox(real(tr(ρ_neel)), 1.0; atol=1e-12)

6-element MPO{Hilbert}
  site dims: 2, 2, 2, 2, 2, 2
  link dims: 1, 1, 1, 1, 1
  maxlinkdim: 1
  combiners: none

  tensors:
    [1] ((dim=2|id=691|"S=1/2,Site,n=1")', (dim=2|id=691|"S=1/2,Site,n=1"), (dim=1|id=423|"Link,l=1"))
    [2] ((dim=2|id=882|"S=1/2,Site,n=2")', (dim=2|id=882|"S=1/2,Site,n=2"), (dim=1|id=861|"Link,l=2"), (dim=1|id=423|"Link,l=1"))
    [3] ((dim=2|id=772|"S=1/2,Site,n=3")', (dim=2|id=772|"S=1/2,Site,n=3"), (dim=1|id=70|"Link,l=3"), (dim=1|id=861|"Link,l=2"))
    [4] ((dim=2|id=716|"S=1/2,Site,n=4")', (dim=2|id=716|"S=1/2,Site,n=4"), (dim=1|id=121|"Link,l=4"), (dim=1|id=70|"Link,l=3"))
    [5] ((dim=2|id=698|"S=1/2,Site,n=5")', (dim=2|id=698|"S=1/2,Site,n=5"), (dim=1|id=179|"Link,l=5"), (dim=1|id=121|"Link,l=4"))
    [6] ((dim=2|id=391|"S=1/2,Site,n=6")', (dim=2|id=391|"S=1/2,Site,n=6"), (dim=1|id=179|"Link,l=5"))

Tr(ρ_neel) = 1.0

A mixed density MPO

Let us mix two product states:

\[|\psi_1\rangle = |\uparrow\downarrow\uparrow\downarrow\cdots\rangle,\]

and

\[|\psi_2\rangle = |\downarrow\uparrow\downarrow\uparrow\cdots\rangle.\]

Define

\[\rho = p|\psi_1\rangle\langle\psi_1| + (1-p)|\psi_2\rangle\langle\psi_2|.\]

flipped_pattern = ["Dn", "Up", "Dn", "Up", "Dn", "Up"]
ψ_flipped = MPS(sites, flipped_pattern)

p = 0.3

ρ_mix = to_dm([ψ_neel, ψ_flipped]; coeffs=[p, 1 - p])

println("Tr(ρ_mix) = ", roundreal(tr(ρ_mix)))

@assert ρ_mix isa MPO{Hilbert}
@assert isapprox(real(tr(ρ_mix)), 1.0; atol=1e-12)

Tr(ρ_mix) = 1.0

One-site reduced density matrix

Let us compute the reduced density matrix at one site.

For a pure state,

\[\rho_j = \operatorname{Tr}_{\bar{j}} |\psi\rangle\langle\psi|.\]

For an MPS, we can do this by constructing the pure density matrix and contracting the bra/ket indices of the sites we wish to trace out.

function one_site_reduced_density_matrix(ρ, sites, site_number::Integer=1)
    ρ_tensor = ITensor(1.0)
    for n in eachindex(ρ)
        ρ_tensor *= ρ[n]
        if n != site_number
            ρ_tensor *= delta(sites[n], prime(sites[n]))
        end
    end
    return Array(ρ_tensor, prime(sites[site_number]), sites[site_number])
end

ρ1_red = one_site_reduced_density_matrix(ρ_mix, sites, 1)

@assert isapprox(LinearAlgebra.tr(ρ1_red), 1.0 + 0im; atol=1e-12)
display(roundreal.(ρ1_red))

Physical entropy of the one-site mixed state

The physical von Neumann entropy is

\[S(\rho) = -\operatorname{Tr}(\rho \log \rho).\]

Since ρ1_red is a 2 × 2 dense matrix, we can compute its eigenvalues directly.

function von_neumann_entropy(ρ::AbstractMatrix)
    vals = real.(LinearAlgebra.eigvals(LinearAlgebra.Hermitian(ρ)))
    return -sum(λ -> λ > 1e-14 ? λ * log(λ) : 0.0, vals)
end

ρ1_mix = Matrix(ρ1_red)

S_ρ1_mix = von_neumann_entropy(ρ1_mix)

S_expected = -(p * log(p) + (1 - p) * log(1 - p))

println("S(ρ1_mix) = ", roundreal(S_ρ1_mix))
println("Expected binary entropy = ", roundreal(S_expected))

@assert isapprox(S_ρ1_mix, S_expected; atol=1e-12)

S(ρ1_mix) = 0.610864
Expected binary entropy = 0.610864

MPO bond entropy as a compression diagnostic

We can still call entropy on the density MPO.

This computes an SVD-based bond entropy of the MPO tensor network. It tells us how much operator-space structure the density MPO carries across a cut.

middle_bond = div(N, 2)
S_operator_bond = entropy(ρ_mix, middle_bond)

println("MPO bond entropy of ρ_mix across bond $middle_bond:")
println(roundreal(S_operator_bond))

MPO bond entropy of ρ_mix across bond 3:
0.566257

Summary

In this tutorial, we learned the basic ITensorMPS workflow used throughout ProcessTensors.jl for Hilbert-space states and operators.

We saw that:

• dense coefficient tensors can be decomposed into MPS form,
• siteinds creates physical spin-chain site indices,
• MPS creates Hilbert-space many-body states,
• OpSum expresses many-body operators compactly,
• MPO converts an OpSum into a matrix product operator,
• expect computes local expectation values,
• correlation_matrix computes two-point functions,
• entropy(ψ, b) computes pure-state bond entanglement,
• to_dm constructs Hilbert-space density MPOs.

The main convention to remember is:

Reuse the correct indices. ITensors makes contractions safe because tensor legs are labelled objects, not anonymous array axes.

In the next tutorial, Liouville-Space Basics, we take a density MPO and vectorize it into Liouville space.