Compressing MPS

Overview

In this notebook we will show how to compress the wavefuntion stored in an MPS via singular value decomposition (SVD).

The model

In this tutorial we will skip the process of building the model from scratch and use one of the pre-built models in Renormalizer that is frequently used for testing. The model we use is a Frenkel-Holstein model with 3 electronic sites and 2 vibration modes for each electronic site

\[\hat H = \sum_{i=1}^3\sum_{j=1}^3 J_{ij} \hat a^\dagger_i \hat a_{i+1} + \sum_{i=1}^3\sum_{k=1,2}\frac{1}{2} (\hat p_{ik}^2 + \omega_k^2 \hat x_{ik}^2) + \sum_{i=1}^3\sum_{k=1,2} \hat a^\dagger_i \hat a_i \hat x_{ik}\]

Therefore, the model has 9 sites in total. 3 of them are electronic sites and the rest are vibrational sites.

from renormalizer.tests.parameter import holstein_model

2024-08-24 09:07:11,933[INFO] Use NumPy as backend
2024-08-24 09:07:11,934[INFO] numpy random seed is 9012
2024-08-24 09:07:11,935[INFO] random seed is 1092
2024-08-24 09:07:11,946[INFO] Git Commit Hash: 78516d8ee697e5cfb88a1cf2c55ad3299ff3b640
2024-08-24 09:07:11,948[INFO] use 64 bits

holstein_model.basis

[BasisSimpleElectron(dof: 0, nbas: 2, qn: [[0], [1]]),
 BasisSHO(dof: (0, 0), x0: 0.0, omega: 0.0004852952677876329, nbas: 4),
 BasisSHO(dof: (0, 1), x0: 0.0, omega: 0.007087607302666907, nbas: 4),
 BasisSimpleElectron(dof: 1, nbas: 2, qn: [[0], [1]]),
 BasisSHO(dof: (1, 0), x0: 0.0, omega: 0.0004852952677876329, nbas: 4),
 BasisSHO(dof: (1, 1), x0: 0.0, omega: 0.007087607302666907, nbas: 4),
 BasisSimpleElectron(dof: 2, nbas: 2, qn: [[0], [1]]),
 BasisSHO(dof: (2, 0), x0: 0.0, omega: 0.0004852952677876329, nbas: 4),
 BasisSHO(dof: (2, 1), x0: 0.0, omega: 0.007087607302666907, nbas: 4)]

holstein_model.ham_terms

[Op('a^\\dagger a', [0, 0], 0.10016074648883302, [[1], [-1]]),
 Op('a^\\dagger a', [0, 1], -0.003674932217565499, [[1], [-1]]),
 Op('a^\\dagger a', [0, 2], -0.007349864435130998, [[1], [-1]]),
 Op('a^\\dagger a', [1, 0], -0.003674932217565499, [[1], [-1]]),
 Op('a^\\dagger a', [1, 1], 0.10016074648883302, [[1], [-1]]),
 Op('a^\\dagger a', [1, 2], -0.011024796652696497, [[1], [-1]]),
 Op('a^\\dagger a', [2, 0], -0.007349864435130998, [[1], [-1]]),
 Op('a^\\dagger a', [2, 1], -0.011024796652696497, [[1], [-1]]),
 Op('a^\\dagger a', [2, 2], 0.10016074648883302, [[1], [-1]]),
 Op('p^2', [(0, 0)], 0.5),
 Op('x^2', [(0, 0)], 1.1775574846853516e-07),
 Op('p^2', [(0, 1)], 0.5),
 Op('x^2', [(0, 1)], 2.5117088638408635e-05),
 Op('p^2', [(1, 0)], 0.5),
 Op('x^2', [(1, 0)], 1.1775574846853516e-07),
 Op('p^2', [(1, 1)], 0.5),
 Op('x^2', [(1, 1)], 2.5117088638408635e-05),
 Op('p^2', [(2, 0)], 0.5),
 Op('x^2', [(2, 0)], 1.1775574846853516e-07),
 Op('p^2', [(2, 1)], 0.5),
 Op('x^2', [(2, 1)], 2.5117088638408635e-05),
 Op('a^\\dagger a x', [0, 0, (0, 0)], -7.097609983192488e-06, [[1], [-1], [0]]),
 Op('a^\\dagger a x', [0, 0, (0, 1)], -0.00044069941383179025, [[1], [-1], [0]]),
 Op('a^\\dagger a x', [1, 1, (1, 0)], -7.097609983192488e-06, [[1], [-1], [0]]),
 Op('a^\\dagger a x', [1, 1, (1, 1)], -0.00044069941383179025, [[1], [-1], [0]]),
 Op('a^\\dagger a x', [2, 2, (2, 0)], -7.097609983192488e-06, [[1], [-1], [0]]),
 Op('a^\\dagger a x', [2, 2, (2, 1)], -0.00044069941383179025, [[1], [-1], [0]])]

Compressing MPS

Firstly construct a random MPS and the Hamiltonian MPO

from renormalizer import Mps, Mpo

mps = Mps.random(holstein_model, qntot=1, m_max=5)
mps.bond_dims

[1, 2, 5, 5, 5, 5, 5, 5, 5, 1]

mpo = Mpo(holstein_model)
mpo.bond_dims

2024-08-24 09:07:12,029[DEBUG] # of operator terms: 27
2024-08-24 09:07:12,029[DEBUG] Input operator terms: 27
2024-08-24 09:07:12,031[DEBUG] After combination of the same terms: 27
2024-08-24 09:07:12,032[DEBUG] symbolic mpo algorithm: qr

[1, 4, 5, 4, 5, 5, 4, 3, 3, 1]

After applying the MPO on the MPS, the resulting MPS has a larger bond dimension. More specifically, the bond dimension is the product of the corresponding bond dimension of MPS and MPO

mps2 = mpo @ mps
mps2.bond_dims

[1, 8, 25, 20, 25, 25, 20, 15, 15, 1]

import numpy as np
np.array(mps.bond_dims) * np.array(mpo.bond_dims)  # the same as the bond dimension of `mps2`

array([ 1,  8, 25, 20, 25, 25, 20, 15, 15,  1])

We then try to compress the MPS. We first make a copy for future reference since the compress is in place

mps3 = mps2.copy()

Then compress mps2 by SVD truncation to a fixed bond dimension of \(M=5\). Note that as the product of an MPO and an MPS, mps2 is not canonicalised, so it has to be canonicalised before SVD compression.

After the compression, the size of the matrices in mps2 is greatly reduced.

m_trunc = 5
mps3.canonicalise().compress(m_trunc)
mps3.bond_dims

2024-08-24 09:07:12,101[DEBUG] size before/after compress: 47.2KiB/4.2KiB, ratio: 11.177777777777777

[1, 2, 5, 5, 5, 5, 5, 5, 4, 1]

We can see that the information loss, measured by the difference of energy expectations, is relatively small.

e1 = mps3.expectation(mpo)
e2 = mps2.expectation(mpo)
e1, e2, abs((e1 - e2) / e2)

(0.0031058187888397606, 0.003148758678336175, 0.013637084922336527)

As a result of the random generation of the MPS, the compress accuracy is not very impressive. Much higher efficiency can be expected for the real wavefunction of quantum systems.

Compress Configuration

It is sometimes desirable to perform truncation based on the magnitude of the singular values. Also, in production level calculations, MPSs are frequently compressed according to a fixed strategy. We can set the compress_config attribute for MPS to control the compression behavior

from renormalizer.utils import CompressConfig, CompressCriteria

mps3 = mps2.copy()
mps3.compress_config = CompressConfig(CompressCriteria.threshold, threshold=1e-5)

Here the compress criteria is set to the singular value threshold, which is set to 1e-5.

CompressCriteria is an Enum class for all possible compression strategies. Currently there’re three different strategies.

[s for s in dir(CompressCriteria) if not s.startswith("__")]

['both', 'fixed', 'threshold']

Using CompressCriteria might seem to be a overkill compared to an implementation using simple strings. However, when it comes to time evolution configurations, since there’re a lot of different algorithms and some of them have rather long names, the Enum class can help manage the different algorithms. So here the compression configuration uses the same implementation. This is also for forward compatibility.

We next compress mps3 using the compress configuration associated with the mps3 object. \(10^{-5}\) is a relatively tight threshold, so the bond dimension is higher than 5.

mps3.canonicalise().compress()
mps3.bond_dims

2024-08-24 09:07:12,184[DEBUG] size before/after compress: 47.2KiB/26.5KiB, ratio: 1.7784325279905715

[1, 2, 8, 15, 17, 18, 11, 8, 4, 1]

We can expect that the accuracy is also higher

e1 = mps3.expectation(mpo)
e2 = mps2.expectation(mpo)
e1, e2, abs((e1 - e2) / e2)

(0.003148758678336171, 0.003148758678336175, 1.377307419517337e-15)