MPS and MPO#
Overview#
In this notebook we will introduce the construction and manipulation of MPS and MPO. We will rely on Op, BasisSet and Model from the previous tutorial.
The spin-boson model#
In this notebook we will consider the MPS and MPO based on a 1-mode spin-boson model.
The spin degree of freedom is labeled as a str "spin" and the boson degree of freedom is labeled as a str "boson". We consider \(\epsilon=0\), \(\Delta=1\), \(\omega=1\) and \(g=0.5\)
[1]:
from renormalizer import Op, BasisHalfSpin, BasisSHO, Model
epsilon = 0
delta = 1
omega = 1
g = 0.5
2024-08-24 09:07:18,241[INFO] Use NumPy as backend
2024-08-24 09:07:18,242[INFO] numpy random seed is 9012
2024-08-24 09:07:18,243[INFO] random seed is 1092
2024-08-24 09:07:18,254[INFO] Git Commit Hash: 78516d8ee697e5cfb88a1cf2c55ad3299ff3b640
2024-08-24 09:07:18,255[INFO] use 64 bits
First, construct all terms in the Hamiltonian in a list or OpSum.
[2]:
ham_terms = [
Op("sigma_z", "spin", epsilon),
Op("sigma_x", "spin", delta),
Op(r"b^\dagger b", "boson", omega),
Op("sigma_z", "spin", g) * Op(r"b^\dagger+b", "boson"),
]
ham_terms
[2]:
[Op('sigma_z', ['spin'], 0.0),
Op('sigma_x', ['spin'], 1.0),
Op('b^\\dagger b', ['boson', 'boson'], 1.0),
Op('sigma_z b^\\dagger+b', ['spin', 'boson'], 0.5)]
Then, define the basis for the model. Note that the boson mode is described by simple harmonic oscillator eigenbasis (BasisSHO) and the number of basis is truncated to 8. This means that the maximum number of bosons for the mode is 8
[3]:
basis = [BasisHalfSpin("spin"), BasisSHO("boson", omega=omega, nbas=8)]
basis
[3]:
[BasisHalfSpin(dof: spin, nbas: 2),
BasisSHO(dof: boson, x0: 0.0, omega: 1, nbas: 8)]
Checkout some of the operators defined by the basis
[4]:
basis[0].op_mat("Z")
[4]:
array([[ 1., 0.],
[ 0., -1.]])
[5]:
basis[1].op_mat(r"b^\dagger b")
[5]:
array([[0., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 2., 0., 0., 0., 0., 0.],
[0., 0., 0., 3., 0., 0., 0., 0.],
[0., 0., 0., 0., 4., 0., 0., 0.],
[0., 0., 0., 0., 0., 5., 0., 0.],
[0., 0., 0., 0., 0., 0., 6., 0.],
[0., 0., 0., 0., 0., 0., 0., 7.]])
Lastly, build our model with Model
[6]:
model = Model(basis, ham_terms)
Construct MPS#
The matrix product state data structure is implemented in the Mps class
[7]:
from renormalizer import Mps
Most commonly, MPS is initialized by two ways. The first is random initialization, which is usually followed by ground state/excited state search by the DMRG algorithm. The second is Hartree product state, which is usually followed by time evolution.
Random Initialization#
Mps.random requires the input of a model. Additionally, the user must also specify the targeted total quantum number qntot and the maximum bond dimension m_max.
Since the spin-boson model does not conserve any quantum number, all states/operators are considered to have quantum number 0 in Renormalizer. Thus, the total quantum number is set to 0.
The maximum bond dimension is set to 2 since a Schmidt decomposition for the 1-mode spin-boson model yields only 2 singular values.
[8]:
mps = Mps.random(model, qntot=0, m_max=2)
mps
[8]:
<class 'renormalizer.mps.mps.Mps'> with 2 sites
The individual local matrices for the MPS can be accessed via simple indexing
[9]:
mps[0], mps[1]
[9]:
(<Matrix at 0x7f40dfe54fd0 (1, 2, 2) float64>,
<Matrix at 0x7f40dfe54d30 (2, 8, 1) float64>)
Matrix can be considered as a wrapper for np.ndarray with more metadata. The indices of the array are left virtual index, the physical index and the right virtual index. The underlaying matrix can be accessed via the array attribute.
[10]:
mps[0].array, mps[1].array
[10]:
(array([[[-0.99496648, -0.10020826],
[-0.10020826, 0.99496648]]]),
array([[[-0.43283409],
[-0.01476036],
[-0.21888015],
[-0.16222924],
[-0.19864132],
[ 0.36843999],
[-0.2858778 ],
[ 0.31116499]],
[[ 0.17636599],
[ 0.2216647 ],
[ 0.45632911],
[ 0.22315835],
[-0.04071718],
[-0.02460395],
[ 0.19383492],
[-0.07963948]]]))
Conceptually, Mps can be considered as a sparse representation for a wavefunction/vector
[11]:
mps.todense()
[11]:
array([ 0.41298208, -0.00752657, 0.17205047, 0.13905035, 0.20172165,
-0.36411992, 0.26501497, -0.3016182 , 0.2188518 , 0.22202805,
0.47596577, 0.23829179, -0.02060673, -0.06140084, 0.22150656,
-0.11041992])
Mps has a lot of utilities to calculate the properties of the wavefunction. For more of them, please refer to the API reference.
[12]:
mps.norm
[12]:
1.0
[13]:
mps.calc_1site_rdm()
[13]:
{0: array([[0.55402862, 0.31394311],
[0.31394311, 0.44597138]]),
1: array([[ 0.21845031, 0.0454829 , 0.17521973, 0.10957589, 0.07879761,
-0.16381269, 0.15792354, -0.14872851],
[ 0.0454829 , 0.04935311, 0.1043828 , 0.05186089, -0.00609355,
-0.01089213, 0.04718602, -0.02224617],
[ 0.17521973, 0.1043828 , 0.25614478, 0.13734241, 0.02489821,
-0.0918717 , 0.15102549, -0.10444965],
[ 0.10957589, 0.05186089, 0.13734241, 0.07611798, 0.02313905,
-0.06526232, 0.08963362, -0.06825227],
[ 0.07879761, -0.00609355, 0.02489821, 0.02313905, 0.04111626,
-0.0721856 , 0.04889473, -0.05856753],
[-0.16381269, -0.01089213, -0.0918717 , -0.06526232, -0.0721856 ,
0.13635338, -0.11009792, 0.11660507],
[ 0.15792354, 0.04718602, 0.15102549, 0.08963362, 0.04889473,
-0.11009792, 0.11929809, -0.10439207],
[-0.14872851, -0.02224617, -0.10444965, -0.06825227, -0.05856753,
0.11660507, -0.10439207, 0.1031661 ]])}
[14]:
# virtual bond dimension
mps.bond_dims
[14]:
[1, 2, 1]
[15]:
# physical bond dimension
mps.pbond_dims
[15]:
[2, 8]
Hartree product state#
Mps.hartree_product_state also requires the input of a model. The state of each degrees of freedom can be specified through the condition argument. The degrees of freedom that are not specified through condition will be set to the ground state, i.e., \([1, 0, 0, \dots, 0]\).
[16]:
import numpy as np
mps2 = Mps.hartree_product_state(model, condition={"spin":[1/np.sqrt(2), 1/np.sqrt(2)]})
mps2
[16]:
<class 'renormalizer.mps.mps.Mps'> with 2 sites
[17]:
mps2[0].array, mps2[1].array
[17]:
(array([[[0.70710678],
[0.70710678]]]),
array([[[1.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.]]]))
Construct MPO#
As described in the previous section, Renormalizer is able to exactly construct the most compact MPO exactly. By feeding the model into Mpo, the Hamiltonian is constructed as MPO
[18]:
from renormalizer import Mpo
mpo = Mpo(model)
mpo[0], mpo[1]
2024-08-24 09:07:18,420[DEBUG] # of operator terms: 3
2024-08-24 09:07:18,420[DEBUG] Input operator terms: 3
2024-08-24 09:07:18,422[DEBUG] After combination of the same terms: 3
2024-08-24 09:07:18,422[DEBUG] symbolic mpo algorithm: qr
[18]:
(<Matrix at 0x7f40dfe45f70 (1, 2, 2, 3) float64>,
<Matrix at 0x7f40dfe45ee0 (3, 8, 8, 1) float64>)
Here the logging outputs that there are 3 terms in the Hamiltonian, since \(\hat \sigma_z\) with coefficient 0 is dropped automatically.
Similar to MPS, it is sometimes conceptionally convenient to consider MPO as a sparse representation of quantum operator/matrix. The interface of MPO is also similar to MPS.
[19]:
# converting to dense operator/matrix
mpo.todense().shape
[19]:
(16, 16)
[20]:
# virtual bond dimension
mpo.bond_dims
[20]:
[1, 3, 1]
Operators other than the Hamiltonian can be constructed by providing the corresponding Op or OpSum. Note that model still has to be provided, since it contains the basis set information.
[21]:
mpo2 = Mpo(model, Op("Z", "spin"))
2024-08-24 09:07:18,446[DEBUG] # of operator terms: 1
2024-08-24 09:07:18,447[DEBUG] Input operator terms: 1
2024-08-24 09:07:18,448[DEBUG] After combination of the same terms: 1
MPO/MPS manipulation#
Two MPSs can be added together, resulting in another MPS
[22]:
mps + mps2
[22]:
<class 'renormalizer.mps.mps.Mps'> with 2 sites
MPO can apply on MPS, resulting in another MPS
[23]:
mpo @ mps
[23]:
<class 'renormalizer.mps.mps.Mps'> with 2 sites
This is equivalent to mpo.apply(mps)
[24]:
mpo.apply(mps).distance(mpo @ mps)
[24]:
0.0
Calculating expectation value between MPS and MPO \(\langle \rm{MPS}|\rm{MPO}|MPS\rangle\) using the expectation function
[25]:
mps.expectation(mpo)
[25]:
2.89309731239322
Calculating transition amplitudes \(\langle \rm{MPS2}|\rm{MPO}|MPS\rangle\) is also straight-forward
[26]:
mps.expectation(mpo, self_conj=mps2)
[26]:
0.3656142066119614