GitHub

MPSKit.jl

Efficient and versatile tools for working with matrix product states

Table of contents

Installation

MPSKit.jl is a part of the general registry, and can be installed via the package manager as:

pkg> add MPSKit

Key Features

  • Construction and manipulation of Matrix Product States (MPS)
  • Calculation of observables and expectation values
  • Various optimization methods for obtaining MPS fixed points
  • Support for both finite and infinite MPS
  • Support for wide variety of symmetries, including Abelian, non-Abelian, fermionic and anyonic symmetries

Usage

To get started with MPSKit, we recommend also including TensorKit.jl and MPSKitModels.jl. The former defines the tensor backend which is used throughout MPSKit, while the latter includes some common operators and models.

using TensorOperations
using TensorKit
using MPSKit
using MPSKitModels
using LinearAlgebra: norm

Finite Matrix Product States

Finite MPS are characterised by a set of tensors, one for each site, which each have 3 legs. They can be constructed by specifying the virtual spaces and the physical spaces, i.e. the dimensions of each of the legs. These are then contracted to form the MPS. In MPSKit, they are represented by FiniteMPS, which can be constructed either by passing in the tensors directly, or by specifying the dimensions of the legs.

d = 2 # physical dimension
D = 5 # virtual dimension
L = 10 # number of sites

mps = FiniteMPS(L, ComplexSpace(d), ComplexSpace(D)) # random MPS with maximal bond dimension D
10-site FiniteMPS (ComplexF64, TensorKit.ComplexSpace):
┌ CL[11]: TensorMap(ℂ^1 ← ℂ^1)
├── AL[10]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^1)
├── AL[9]: TensorMap((ℂ^4 ⊗ ℂ^2) ← ℂ^2)
├── AL[8]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^4)
├── AL[7]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
├── AL[6]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
├── AL[5]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
├── AL[4]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
├── AL[3]: TensorMap((ℂ^4 ⊗ ℂ^2) ← ℂ^5)
├── AL[2]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^4)
└── AL[1]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^2)

The FiniteMPS object then handles the gauging of the MPS, which is necessary for many of the algorithms. This is done automatically when needed, and the user can access the gauged tensors by getting and setting the AL, AR, CR/CL and AC fields, which each represent a vector of these tensors.

al = mps.AL[3] # left gauged tensor of the third site
@tensor E[a; b] := al[c, d, b] * conj(al[c, d, a])
@show isapprox(E, id(left_virtualspace(mps, 3)))
true
ar = mps.AR[3] # right gauged tensor of the third site
@tensor E[a; b] := ar[a, d, c] * conj(ar[b, d, c])
@show isapprox(E, id(right_virtualspace(mps, 2)))
true

As the mps will be kept in a gauged form, updating a tensor will also update the gauged tensors. For example, we can set the tensor of the third site to the identity, and the gauged tensors will be updated accordingly.

mps.CR[3] = id(domain(mps.CR[3]))
println(mps)
┌── AR[10]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^1)
├── AR[9]: TensorMap((ℂ^4 ⊗ ℂ^2) ← ℂ^2)
├── AR[8]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^4)
├── AR[7]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
├── AR[6]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
├── AR[5]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
├── AR[4]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
│ CL[4]: TensorMap(ℂ^5 ← ℂ^5)
├── AL[3]: TensorMap((ℂ^4 ⊗ ℂ^2) ← ℂ^5)
├── AL[2]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^4)
└── AL[1]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^2)

These objects can then be used to compute observables and expectation values. For example, the expectation value of the identity operator at the third site, which is equal to the norm of the MPS, can be computed as:

N1 = LinearAlgebra.norm(mps)
N2 = expectation_value(mps, 3 => id(physicalspace(mps, 3)))
println("‖mps‖ = $N1")
println("<mps|𝕀₃|mps> = $N2")
‖mps‖ = 2.2360679774997894
<mps|𝕀₃|mps> = 5.000000000000001 + 0.0im

Finally, the MPS can be optimized in order to determine groundstates of given Hamiltonians. Using the pre-defined models in MPSKitModels, we can construct the groundstate for the transverse field Ising model:

H = transverse_field_ising(; J=1.0, g=0.5)
find_groundstate!(mps, H, DMRG(; maxiter=10))
E0 = expectation_value(mps, H)
println("<mps|H|mps> = $(sum(real(E0)) / length(mps))")
[ Info: DMRG init:	obj = -2.733816237983e+00	err = 4.8070e-01
[ Info: DMRG   1:	obj = -9.765501602864e+00	err = 2.4968668312e-02	time = 0.50 sec
[ Info: DMRG   2:	obj = -9.765503468882e+00	err = 1.0879268589e-04	time = 0.32 sec
[ Info: DMRG   3:	obj = -9.765503492654e+00	err = 7.5131230665e-06	time = 0.02 sec
[ Info: DMRG   4:	obj = -9.765503493140e+00	err = 3.3658030475e-06	time = 0.01 sec
[ Info: DMRG   5:	obj = -9.765503493261e+00	err = 1.5349083930e-06	time = 0.01 sec
[ Info: DMRG   6:	obj = -9.765503493294e+00	err = 8.4153324601e-07	time = 0.01 sec
[ Info: DMRG   7:	obj = -9.765503493304e+00	err = 5.0942677133e-07	time = 0.01 sec
[ Info: DMRG   8:	obj = -9.765503493307e+00	err = 3.1298953463e-07	time = 0.01 sec
[ Info: DMRG   9:	obj = -9.765503493308e+00	err = 1.9413796258e-07	time = 0.01 sec
┌ Warning: DMRG cancel 10:	obj = -9.765503493309e+00	err = 1.2114747492e-07	time = 0.94 sec
└ @ MPSKit ~/work/MPSKit.jl/MPSKit.jl/src/algorithms/groundstate/dmrg.jl:54
<mps|H|mps> = -0.9765503493308904

Infinite Matrix Product States

Similarly, an infinite MPS can be constructed by specifying the tensors for the unit cell, characterised by the spaces (dimensions) thereof.

d = 2 # physical dimension
D = 5 # virtual dimension
mps = InfiniteMPS(d, D) # random MPS
single site InfiniteMPS:
│   ⋮
│ CR[1]: TensorMap(ℂ^5 ← ℂ^5)
├── AL[1]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
│   ⋮

The InfiniteMPS object then handles the gauging of the MPS, which is necessary for many of the algorithms. This is done automatically upon creation of the object, and the user can access the gauged tensors by getting and setting the AL, AR, CR/CL and AC fields, which each represent a (periodic) vector of these tensors.

al = mps.AL[1] # left gauged tensor of the third site
@tensor E[a; b] := al[c, d, b] * conj(al[c, d, a])
@show isapprox(E, id(left_virtualspace(mps, 1)))
true
ar = mps.AR[1] # right gauged tensor of the third site
@tensor E[a; b] := ar[a, d, c] * conj(ar[b, d, c])
@show isapprox(E, id(right_virtualspace(mps, 2)))
true

As regauging the MPS is not possible without recomputing all the tensors, setting a single tensor is not supported. Instead, the user should construct a new mps object with the desired tensor, which will then be gauged upon construction.

als = 3 .* mps.AL
mps = InfiniteMPS(als)
single site InfiniteMPS:
│   ⋮
│ CR[1]: TensorMap(ℂ^5 ← ℂ^5)
├── AL[1]: TensorMap((ℂ^5 ⊗ ℂ^2) ← ℂ^5)
│   ⋮

These objects can then be used to compute observables and expectation values. For example, the norm of the MPS, which is equal to the expectation value of the identity operator can be computed by:

N1 = norm(mps)
N2 = expectation_value(mps, 1 => id(physicalspace(mps, 1)))
println("‖mps‖ = $N1")
println("<mps|𝕀₁|mps> = $N2")
‖mps‖ = 0.9999999999999999
<mps|𝕀₁|mps> = 0.9999999999999996 + 0.0im
Normalization of infinite MPS

Because infinite MPS cannot sensibly be normalized to anything but $1$, the norm of an infinite MPS is always set to be $1$ at construction. If this were not the case, any observable computed from the MPS would either blow up to infinity or vanish to zero.

Finally, the MPS can be optimized in order to determine groundstates of given Hamiltonians. Using the pre-defined models in MPSKitModels, we can construct the groundstate for the transverse field Ising model:

H = transverse_field_ising(; J=1.0, g=0.5)
mps, = find_groundstate(mps, H, VUMPS(; maxiter=10))
E0 = expectation_value(mps, H)
println("<mps|H|mps> = $(sum(real(E0)) / length(mps))")
[ Info: VUMPS init:	obj = -5.021876093998e-01	err = 4.5075e-01
[ Info: VUMPS   1:	obj = -1.038920289822e+00	err = 1.0636274549e-01	time = 0.00 sec
┌ Warning: ignoring imaginary component 1.1738374042363775e-6 from total weight 1.0888813436641487: operator might not be hermitian?
│   α = -0.10529135664495481 + 1.1738374042363775e-6im
│   β₁ = 0.6939250475606353
│   β₂ = 0.8324928463129875
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.7400000109531188e-6 from total weight 2.6230979783673254: operator might not be hermitian?
│   α = 2.1015354940385427 + 1.7400000109531188e-6im
│   β₁ = 0.8324928463129875
│   β₂ = 1.330844556004044
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.786125635887259e-6 from total weight 3.2318109568483693: operator might not be hermitian?
│   α = 2.3415739860724325 + 1.786125635887259e-6im
│   β₁ = 1.330844556004044
│   β₂ = 1.7861931856059805
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.31521591536371e-6 from total weight 1.3576410523133957: operator might not be hermitian?
│   α = 0.9220754287035805 + 1.31521591536371e-6im
│   β₁ = 0.6571534675491175
│   β₂ = 0.7490763985029655
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 1.9047184168022313e-6 from total weight 3.382197752602955: operator might not be hermitian?
│   α = 3.01972789143682 + 1.9047184168022313e-6im
│   β₁ = 0.7490763985029655
│   β₂ = 1.3264198613535667
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 2.2681652829537047e-6 from total weight 3.5912447010893533: operator might not be hermitian?
│   α = 3.0329402312390186 + 2.2681652829537047e-6im
│   β₁ = 1.3264198613535667
│   β₂ = 1.3924519410865892
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
┌ Warning: ignoring imaginary component 2.370950433157326e-6 from total weight 4.325466307153559: operator might not be hermitian?
│   α = 3.859949253623568 + 2.370950433157326e-6im
│   β₁ = 1.3924519410865892
│   β₂ = 1.368038057047419
└ @ KrylovKit ~/.julia/packages/KrylovKit/xccMN/src/factorizations/lanczos.jl:170
[ Info: VUMPS   2:	obj = -1.063542783231e+00	err = 1.0474740843e-03	time = 0.08 sec
[ Info: VUMPS   3:	obj = -1.063544409917e+00	err = 6.5701602717e-06	time = 0.00 sec
[ Info: VUMPS   4:	obj = -1.063544409973e+00	err = 7.1750093965e-07	time = 0.00 sec
[ Info: VUMPS   5:	obj = -1.063544409973e+00	err = 6.0977966998e-08	time = 0.00 sec
[ Info: VUMPS   6:	obj = -1.063544409973e+00	err = 6.0533436026e-09	time = 0.00 sec
[ Info: VUMPS   7:	obj = -1.063544409973e+00	err = 6.4589017660e-10	time = 0.00 sec
[ Info: VUMPS conv 8:	obj = -1.063544409973e+00	err = 7.6260743024e-11	time = 0.10 sec
<mps|H|mps> = -1.063544409973251

Additional Resources

For more detailed information on the functionality and capabilities of MPSKit, refer to the Manual section, or have a look at the Examples page.

Keep in mind that the documentation is still a work in progress, and that some features may not be fully documented yet. If you encounter any issues or have questions, please check the library's issue tracker on the GitHub repository and open a new issue.

Publications using MPSKit

Below you can find a list of publications that have made use of MPSKit. If you have used this package and wish to have your publication added to this list, please open a pull request or an issue on the GitHub repository.

  • R. Belyansky et al., “High-Energy Collision of Quarks and Hadrons in the Schwinger Model: From Tensor Networks to Circuit QED,” 2023, doi: 10.48550/ARXIV.2307.02522.
  • L. Devos, L. Vanderstraeten, and F. Verstraete, “Haldane gap in the SU(3) [3 0 0] Heisenberg chain,” Phys. Rev. B, vol. 106, no. 15, p. 155103, Oct. 2022, doi: 10.1103/PhysRevB.106.155103.
  • J. C. Halimeh, M. V. Damme, T. V. Zache, D. Banerjee, and P. Hauke, “Achieving the quantum field theory limit in far-from-equilibrium quantum link models,” Quantum, vol. 6, p. 878, Dec. 2022, doi: 10.22331/q-2022-12-19-878.
  • J. C. Halimeh, D. Trapin, M. Van Damme, and M. Heyl, “Local measures of dynamical quantum phase transitions,” Phys. Rev. B, vol. 104, no. 7, p. 075130, Aug. 2021, doi: 10.1103/PhysRevB.104.075130.
  • M. Hauru, M. Van Damme, and J. Haegeman, “Riemannian optimization of isometric tensor networks,” SciPost Physics, vol. 10, no. 2, p. 040, Feb. 2021, doi: 10.21468/SciPostPhys.10.2.040.
  • M. Van Damme, R. Vanhove, J. Haegeman, F. Verstraete, and L. Vanderstraeten, “Efficient matrix product state methods for extracting spectral information on rings and cylinders,” Phys. Rev. B, vol. 104, no. 11, p. 115142, Sep. 2021, doi: 10.1103/PhysRevB.104.115142.
  • M. Van Damme, T. V. Zache, D. Banerjee, P. Hauke, and J. C. Halimeh, “Dynamical quantum phase transitions in spin-$S$ $\text{U}(1)$ quantum link models,” Phys. Rev. B, vol. 106, no. 24, p. 245110, Dec. 2022, doi: 10.1103/PhysRevB.106.245110.
  • E. L. Weerda and M. Rizzi, “Fractional quantum Hall states with variational Projected Entangled-Pair States: a study of the bosonic Harper-Hofstadter model,” 2023, doi: 10.48550/ARXIV.2309.12811.
  • C. Yu and J.-W. Lee, “Closing of the Haldane gap in a spin-1 XXZ chain,” J. Korean Phys. Soc., vol. 79, no. 9, pp. 841–845, Nov. 2021, doi: 10.1007/s40042-021-00283-z.
  • Y. Zhang, A. Hulsch, H.-C. Zhang, W. Tang, L. Wang, and H.-H. Tu, “Universal Scaling of Klein Bottle Entropy near Conformal Critical Points,” Phys. Rev. Lett., vol. 130, no. 15, p. 151602, Apr. 2023, doi: 10.1103/PhysRevLett.130.151602.
  • Gertian Roose, Laurens Vanderstraeten, Jutho Haegeman, and Nick Bultinck. Anomalous domain wall condensation in a modified ising chain. Phys. Rev. B, 99: 195132, May 2019. 10.1103/​PhysRevB.99.195132.

https:/​/​doi.org/​10.1103/​PhysRevB.99.195132

  • Roose, G., Bultinck, N., Vanderstraeten, L. et al. Lattice regularisation and entanglement structure of the Gross-Neveu model. J. High Energ. Phys. 2021, 207 (2021). https://doi.org/10.1007/JHEP07(2021)207
  • Roose, G., Haegeman, J., Van Acoleyen, K. et al. The chiral Gross-Neveu model on the lattice via a Landau-forbidden phase transition. J. High Energ. Phys. 2022, 19 (2022). https://doi.org/10.1007/JHEP06(2022)019