Algorithms for MPS and MPO

Summing MPS

• sum(MPS A1, MPS A2, Args args = Args::global()) -> MPS
  

  Return the sum of the MPS A1 and A2. The returned MPS will have an orthogonality center on site 1. Before being returned, the MPS representing the sum will be compressed using truncation parameters provided in the named arguments args.

  The input MPS must have the same site indices. The link indices of the output MPS will have the same tags as the first input MPS.

  Show Example

  auto A3 = sum(A1,A2,{"MaxDim",500,"Cutoff",1E-8});
  

• sum(vector<MPS> terms, Args args = Args::global()) -> MPS

  Returns the sum of all the MPS provided in the vector terms as a single MPS, using the truncation accuracy parameters (such as "Cutoff" or "MaxDim") provided in the named arguments args to control the accuracy of the sum.

  This function uses a hierarchical, tree-like algorithm which first sums pairs of MPS, then pairs of pairs, etc. so that the largest bond dimensions are only reached toward the end of the process for maximum efficiency. Therefore using this algorithm can be much faster than calling the above two-argument sum function to sum the terms one at a time.

  Show Example

  auto terms = vector<MPS>(4);
  terms.at(0) = A0;
  terms.at(1) = A1;
  terms.at(2) = A2;
  terms.at(3) = A3;
  
  auto res = sum(terms,{"Cutoff",1E-8});
  

Measuring Properties of MPS

• expect(MPS psi, SiteSet sites, string A) -> vector<Real>
expectC(MPS psi, SiteSet sites, string A) -> vector<Cplx>
expect(MPS psi, SiteSet sites, string A, vector<int> site_list) -> vector<Real>
expectC(MPS psi, SiteSet sites, string A, vector<int> site_list) -> vector<Cplx>

  Compute the expected value @@\langle\psi| \hat{A} |\psi\rangle@@ for an operator A on every site, returning a vector of the results. The operator can be any operator name recognized by the SiteSet, such as "Sz" for the SpinHalf SiteSet. The function expectC returns a vector of complex values and should be used when complex results are expected.

  You can pass a vector site_list of integers as an optional last argument to obtain expectation values only for sites in that vector.

  Show Example

  auto exsz = expect(psi,sites,"Sz");
  
  //Measure only sites 3,5,8
  auto exsz = expect(psi,sites,"Sz",{{3,5,8}});
  

• expect(MPS psi, SiteSet sites, vector<string> ops) -> vector<vector<Real>>
expectC(MPS psi, SiteSet sites, vector<string> ops) -> vector<vector<Cplx>> expect(MPS psi, SiteSet sites, vector<string> ops, vector<int> site_list) -> vector<vector<Real>>
expectC(MPS psi, SiteSet sites, vector<string> ops, vector<int> site_list) -> vector<vector<Cplx>>

  Compute the expected value of a set of operators ops for each site of the MPS psi. Similar to calling the single-operator version of expect above but is more efficient. The returned value is a vector of vectors, one for each operator, containing the expected values on each site.

  You can pass a vector site_list of integers as an optional last argument to obtain expectation values only for sites in that vector.

  Show Example

  auto ex = expect(psi,sites,{"Sz","Sx","Sy"});
  
  //Measure only sites 3,5,8
  auto ex = expect(psi,sites,{"Sz","Sx","Sy"},{{3,5,8}});
  

• correlationMatrix(MPS psi, SiteSet sites, string A, string B) -> vector<vector<Real>> correlationMatrixC(MPS psi, SiteSet sites, string A, string B) -> vector<vector<Cplx>>

  Given an MPS psi, a SiteSet, and two strings denoting operators computes the two-point correlation function matrix @@M_{ij} = \langle\psi| \hat{A}_i \hat{B}_j |\psi\rangle@@
  using efficient MPS techniques. Returns the matrix M.

  Show Example

  auto czz = correlationMatrix(psi,sites,"Sz","Sz");
  auto cpm = correlationMatrix(psi,sites,"S+","S-");
  

Inner Products and Expectation Values

• inner(MPS y, MPS x) -> Real
innerC(MPS y, MPS x) -> Cplx

  Compute the exact inner product @@\langle y|x \rangle@@ of two MPS (the tensors of y will get conjugated). If the inner product is expected to be a complex number use innerC.

  The algorithm used scales as @@m^3 d@@ where @@m@@ is a typical link dimension of the MPS and @@d@@ is the site dimension.

  Note that if x and y don't have the same site indices, this function will attempt to make them match.

• inner(MPS y, MPO A, MPS x) -> Real
innerC(MPS y, MPO A, MPS x) -> Cplx

  Compute the exact inner product @@\langle y|A|x \rangle@@ of two MPS y and x with respect to an MPO A (the tensors of y will get conjugated).

  The algorithm used scales as @@m^3\, k\,d + m^2\, k^2\, d^2@@ where @@m@@ is typical link dimension of the MPS, @@k@@ is the typical MPO dimension, and @@d@@ is the site dimension.

  Note that A and x must share a set of site indices. If the remaining site indices of A are not shared with y, this function will attempt to match them (i.e. it has the same behavior as inner(y,Ax) if Ax was the exact application of MPO A to MPS x).

• inner(MPS y, MPO B, MPO A, MPS x) -> Real
innerC(MPS y, MPO B, MPO A, MPS x) -> Cplx

  Compute the exact inner product @@\langle y|BA|x \rangle@@ of two MPS y and x with respect to two MPOs B and A (MPS y will get conjugated).

  MPO A must share a set of site indices with MPS x, and the other set of site indices with MPO B. If the remaining set of site indices of B are not shared with y, with function will attempt to make them match.

  The algorithm used scales as @@m^3\, k^2\,d + m^2\, k^3\, d^2@@ where @@m@@ is typical bond dimension of the MPS, @@k@@ is the typical MPO dimension, and @@d@@ is the site dimension.

• inner(MPO B, MPS y, MPO A, MPS x) -> Real
innerC(MPO B, MPS y, MPO A, MPS x) -> Cplx

  Compute the exact inner product @@\langle By|A|x \rangle@@ (i.e. the inner product of of @@B|y \rangle@@ and @@A|x \rangle@@ ). MPO B and MPS y will get conjugated.

  MPO A must share a set of site indices with MPS x, and MPO B must share a set of site indices with MPS y. If the remaining site indices of A and B do not match with each other, the function will attempt to make them match.

  The algorithm used scales as @@m^3\, k^2\,d + m^2\, k^3\, d^2@@ where @@m@@ is typical bond dimension of the MPS, @@k@@ is the typical MPO dimension, and @@d@@ is the site dimension.

Tracing an MPO

• trace(MPO A) -> Real

  traceC(MPO A) -> Cplx

  Trace over the site indices of the MPO.

• trace(MPO A, MPO B) -> Real

  traceC(MPO A, MPO B) -> Cplx

  Return the trace of the operator that would result from performing the exact contraction of MPO A with MPO B. For each j, A(j) and B(j) must share one or two site indices.

  Note that neither A or B will get conjugated by this function.

Multiplying MPOs

• nmultMPO(MPO A, MPO B, Args args = Args::global()) -> MPO

  Multiply MPOs A and B, returning the results MPO. MPO tensors are multiplied one at a time from left to right and the resulting tensors are compressed using the truncation parameters (such as "Cutoff" and "MaxDim") provided through the named arguments args.

  For each j, MPO tensors A(j) and B(j) must share a single site index. MPO C will contain the site indices not shared by MPOs A and B. In addition, the link indices of MPO C will have the same tags as the link indices of the MPO A.

  Show Example

  auto sites = SiteSet(10,2);
  
  // Make trivial MPOs
  auto A = MPO(sites);
  auto B = MPO(sites);
  
  // Prime MPO A to ensure only one set of site indices are shared
  auto C = nmultMPO(prime(A),B,{"MaxDim",500,"Cutoff",1E-8});
  
  auto s3 = sites(3);
  Print(hasInds(C(3),{s3,prime(s3,2)})); //print: true
  

Applying MPO to MPS

• applyMPO(MPO A, MPS x, Args args = Args::global()) -> MPS

  Apply an MPO A to an MPS x, resulting in an approximation to the MPS y:
  @@|y\rangle = A |x\rangle@@ .
  The resulting MPS is returned. The algorithm used is chosen with the parameter "Method" in the named arguments args.

  MPO A and MPS x must share a set of site indices. The links of the output MPS will have the same tags as the links of the input MPS x. The site indices of the output MPS will be the site indices of A that are not shared with x.

  The default algorithm used is the "density matrix" algorithm, chosen by setting the parameter "Method" to "DensityMatrix". If the input MPS has a typical bond dimension of @@m@@ and the MPO has typical bond dimension @@k@@ , this algorithm scales as @@m^3 k^2 + m^2 k^3@@ .

  No approximation is made when applying the MPO, but after applying it the resulting MPS is compressed using the truncation parameters provided in the named arguments args.

  An alternative algorithm can be chosen by setting the parameter "Method" to "Fit". This is a sweeping algorithm that iteratively optimizes the resulting MPS @@|y\rangle@@ (analogous to DMRG). The scaling of the "Fit" method is @@m^3 k + m^2 k^2@@ where @@m@@ is the typical MPS bond dimension and @@k@@ is the typical MPO bond dimension. Thus this algorithm has better scaling in the MPO bond dimension @@k@@ compared to the "DensityMatrix" method, but is not guaranteed to converge (depending on the input MPO and MPS). To ensure convergence it is helpful to use an initial guess for the result MPS that is close to the actual one. Also it can be helpful (if appropriate) if the MPO is close to the identity, such as when time evolving using a small time step. The number of sweeps can be chosen with the parameter "Nsweep".

  It is recommended to try the default "DensityMatrix" first because it is more reliable. Then, the "Fit" method can be tried if higher performance is required.

  Named arguments recognized:

  • "Method" — (default: "DensityMatrix") algorithm used for applying the MPO to the MPS. Currently available options are

    • "DensityMatrix"
    • "Fit"

  • "Cutoff" — (default: 1E-13) truncation error cutoff for compressing resulting MPS

  • "MaxDim" — maximum bond dimension of resulting compressed MPS

  • "Verbose" — (default: false) if true, prints extra output

  • "Normalize" — (default: false) choose whether or not to normalize the output wavefunction

  • "Nsweep" — (default: 1) sets the number of sweeps of the "Fit" algorithm

  Show Example

  //Use the method "DensityMatrix"
  auto y1 = applyMPO(A,x,{"Method=","DensityMatrix","MaxDim=",100,"Cutoff=",1E-8});
  
  //Use the method "Fit" with 5 sweeps
  auto y2 = applyMPO(A,x,{"Method=","Fit","MaxDim=",100,"Cutoff=",1E-8,"Nsweep=",5});
  

• applyMPO(MPO A, MPS x, MPS x0, Args args = Args::global()) -> MPS
  

  Similar to applyMPO above, but accepts a guess for the output wavefunction (the guess wavefunction x0 is not overwritten).

  MPO A and MPS x must share a set of site indices. The site indices of x0 will be made to match the site indices of A that are not shared by x. The links of the output MPS will have the same tags as the links of the guess MPS x0.

  Currently, this version of applyMPO only accepts "Fit" for the parameter "Method". Choosing a good guess state x0 can improve the convergence of the "Fit" method.

  Show Example

  auto sites = SiteSet(10,2);
  
  // Make trivial MPO and random MPS
  auto A = MPO(sites);
  auto x = randomMPS(sites);
  
  // Some other random starting state
  auto x0 = randomMPS(sites);
  
  //Use the method "Fit" with 5 sweeps and a guess state x0
  auto y = applyMPO(A,x,x0,{"Method=","Fit","MaxDim=",100,"Cutoff=",1E-8,"Nsweep=",2});
  

• errorMPOProd(MPS y, MPO A, MPS x) -> Real

  Computes, without approximation, the difference @@||\, |y\rangle - A |x\rangle ||^2@@ , where A is an MPO that shares a set of site indices with MPS x. This is especially useful for testing methods for applying an MPO to an MPS.

  A and x need to share a set of site indices. The function will attempt to match the remaining site indices of A with the site indices of y.

  Show Example

  auto sites = SiteSet(10,2);
  
  // Make trivial MPO and random MPS
  auto A = MPO(sites);
  auto x = randomMPS(sites);
  
  //Approximate A|x>
  auto y = applyMPO(A,x,{"MaxDim=",200,"Cutoff=",1E-12});
  
  //Check 
  Print(errorMPOProd(y,A,x)); //should be close to zero
  

This page current as of version 3.0.0

Back to Classes
Back to Main