[1]:
addpath(genpath('../../../../src'))
relTol = 1e-12;
Matrix product states in the thermodynamic limit¶
This notebook briefly demonstrates the implementation of uniform matrix product states (MPS) directly in the thermodynamic limit, using the Tensor backend. It is based on the second chapter of the lecture notes on tangent space methods for uniform MPS by Laurens Vanderstraeten, Jutho Haegeman and Frank Verstraete.
The contents of this notebook mirror that of a tutorial given at the 2020 school on Tensor Network based approaches to Quantum Many-Body Systems held in Bad Honnef, Germany, which can be found here.
1 Normalization¶
We start by considering a uniform MPS in the thermodynamic limit, which is defined as
Here, \(\boldsymbol{v}_L^\dagger\) and \(\boldsymbol{v}_R\) represent boundary vectors at infinity and the \(A^i\) are complex matrices of size \(D \times D\) for every entry of the index \(i\). This allows for the interpretation of the object \(A\) as a three-legged tensor of dimensions \(D\times d \times D\), where we will refer to \(D\) as the bond dimension and \(d\) as the physical dimension. With this object and the diagrammatic language of tensor networks, we can represent the state as
Thus, we initialize and represent a uniform MPS state as follows:
[2]:
%%file createMPS.m
function A = createMPS(D, d)
% Returns a random complex MPS tensor.
%
% Parameters
% ----------
% D : int
% Bond dimension for MPS.
% d : int
% Physical dimension for MPS.
%
% Returns
% -------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
A = Tensor.randnc([D, d, D]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/createMPS.m'.
[3]:
d = 3;
D = 5;
A = createMPS(D, d);
assert(isequal(A.dims, [D, d, D]), 'Generated MPS tensor has incorrect shape.')
assert(~isreal(A), 'MPS tensor should have complex values.')
One of the central objects in any MPS calculation is the transfer matrix, defined in our case as
This object corresponds to an operator acting on the space of \(D \times D\) matrices, and can be interpreted as a 4-leg tensor. We will use the following convention for ordering the legs: 1. top left 2. bottom left 3. top right 4. bottom right
The transfer matrix can be shown to be a completely positive map, such that its leading eigenvalue is a positive number. This eigenvalue should be rescaled to one to ensure a proper normalization of the state in the thermodynamic limit. To perform this normalization, we must therefore find the left and right fixed points \(l\) and \(r\) which correspond to the largest eigenvalues of the eigenvalue equations
Normalizing the state then means rescaling the MPS tensor \(A \leftarrow A / \sqrt{\lambda}\). Additionally, we may fix the normalization of the eigenvectors by requiring that their trace is equal to one:
With these properties in place, the norm of an MPS can be evaluated as
It can be readily seen that the infinite product of transfer matrices reduces to a projector onto the fixed points, so that the norm reduces to the overlap between the boundary vectors and the fixed points. Since there is no effect of the boundary vectors on the bulk properties of the MPS, we can always choose these such that MPS is properly normalized as $ \left `:nbsphinx-math:langle :nbsphinx-math:Psi`(\bar`{A}):nbsphinx-math:middle | :nbsphinx-math:Psi`(A) \right `:nbsphinx-math:rangle `= 1$.
[4]:
%%file createTransfermatrix.m
function E = createTransfermatrix(A)
% Form the transfermatrix of an MPS.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
% Returns
% -------
% E : :class:`Tensor` (D, D, D, D)
% Transfer matrix with 4 legs,
% ordered topLeft-bottomLeft-topRight-bottomRight.
E = contract(A, [-1, 1, -3], conj(A), [-2, 1, -4], 'Rank', [2, 2]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/createTransfermatrix.m'.
[5]:
%%file normalizeMPS.m
function Anew = normalizeMPS(A)
% Normalize an MPS tensor.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% Anew : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Complexity
% ----------
% O(D^6) algorithm,
% diagonalizing (D^2, D^2) matrix.
% create transfer matrix and reshape to appropriate tensor map acting on right fixed point
E = createTransfermatrix(A);
E = permute(E, [1, 2, 4, 3]);
% determine eigenvalue and rescale MPS tensor accordingly
norm = eigsolve(E);
Anew = A / sqrt(norm);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/normalizeMPS.m'.
[6]:
%%file leftFixedPoint.m
function l = leftFixedPoint(A)
% Find left fixed point.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% l : :class:`Tensor` (D, D)
% left fixed point with 2 legs,
% ordered bottom-top.
%
% Complexity
% ----------
% O(D^6) algorithm,
% diagonalizing (D^2, D^2) matrix.
% initialize transfer matrix and reshape to appropriate tensor map acting on left fixed point
E = createTransfermatrix(A);
E = permute(E, [4, 3, 1, 2]);
% find fixed point
[l, ~] = eigsolve(E);
% interpret as matrix
l = repartition(l, [1, 1]);
% make left fixed point hermitian and positive semidefinite explicitly
l = l * trace(l) / abs(trace(l)); % remove possible phase, actually forgot why I had to do this
l = (l + l') / 2; % force hermitian
l = l * sign(trace(l)); % force positive semidefinite
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/leftFixedPoint.m'.
[7]:
%%file rightFixedPoint.m
function r = rightFixedPoint(A)
% Find right fixed point.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% r : :class:`Tensor` (D, D)
% left fixed point with 2 legs,
% ordered top-bottom.
%
% Complexity
% ----------
% O(D^6) algorithm,
% diagonalizing (D^2, D^2) matrix.
% initialize transfer matrix and reshape to appropriate tensor map acting on right fixed point
E = createTransfermatrix(A);
E = permute(E, [1:E.rank(1), flip(1:E.rank(2)) + E.rank(1)]);
% find fixed point
[r, ~] = eigsolve(E);
% interpret to matrix
r = repartition(r, [1, 1]);
% make right fixed point hermitian and positive semidefinite explicitly
r = r * trace(r) / abs(trace(r)); % remove possible phase, actually forgot why I had to do this
r = (r + r') / 2; % force hermitian
r = r * sign(trace(r)); % force positive semidefinite
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/rightFixedPoint.m'.
[8]:
%%file fixedPoints.m
function [l, r] = fixedPoints(A)
% Find normalized fixed points.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% l : :class:`Tensor` (D, D)
% left fixed point with 2 legs,
% ordered bottom-top.
% r : :class:`Tensor` (D, D)
% right fixed point with 2 legs,
% ordered top-bottom.
%
% Complexity
% ----------
% O(D^6) algorithm,
% diagonalizing (D^2, D^2) matrix
% find fixed points
l = leftFixedPoint(A);
r = rightFixedPoint(A);
% calculate trace and normalize
l = l / trace(l * r);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/fixedPoints.m'.
[9]:
A = normalizeMPS(A);
[l, r] = fixedPoints(A);
assert(isapprox(l, l', 'RelTol', relTol), 'left fixed point should be hermitian!')
assert(isapprox(r, r', 'RelTol', relTol), 'right fixed point should be hermitian!')
assert(isapprox(l, contract(A, [1, 2, -2], l, [3, 1], conj(A), [3, 2, -1]), 'RelTol', relTol), 'l should be a left fixed point!')
assert(isapprox(r, contract(A, [-1, 2, 1], r, [1, 3], conj(A), [-2, 2, 3]), 'RelTol', relTol), 'r should be a right fixed point!')
assert(abs(trace(l * r) - 1) < relTol, 'Left and right fixed points should be trace normalized!') % TODO: fix once traces are a thing
2 Gauge freedom¶
While a given MPS tensor \(A\) corresponds to a unique state \(\left | \Psi(A) \right \rangle\), the converse is not true, as different tensors may give rise to the same state. This is easily seen by noting that the gauge transform
leaves the physical state invariant. We may use this freedom in parametrization to impose canonical forms on the MPS tensor \(A\).
We start by considering the left-orthonormal form of an MPS, which is defined in terms of a tensor \(A_L\) that satisfies the condition
We can find the gauge transform \(L\) that brings \(A\) into this form
by decomposing the fixed point \(l\) as \(l = L^\dagger L\), such that
Note that this gauge choice still leaves room for unitary gauge transformations
which can be used to bring the right fixed point \(r\) into diagonal form. Similarly, we can find the gauge transform that brings \(A\) into right-orthonormal form
such that
and the left fixed point \(l\) is diagonal. The routines that bring a given MPS into canonical form by decomposing the corresponding transfer matrix fixed points can be defined as follows:
[10]:
%%file leftOrthonormalize.m
function [L, Al] = leftOrthonormalize(A, l)
% Transform A to left-orthonormal gauge.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% L : :class:`Tensor` (D, D)
% left gauge with 2 legs,
% ordered left-right.
% Al : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% left orthonormal
%
% Complexity
% ----------
% O(D^6) algorithm,
% diagonalizing (D^2, D^2) matrix
% find left fixed point
if nargin < 2 || isempty(l)
l = leftFixedPoint(A);
end
% decompose l = L * L
L = sqrtm(l);
% apply gauge L to A
Al = contract(L, [-1, 1], A, [1, -2, 2], inv(L), [2, -3]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/leftOrthonormalize.m'.
[11]:
%%file rightOrthonormalize.m
function [R, Ar] = rightOrthonormalize(A, r)
% Transform A to right-orthonormal gauge.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% R : :class:`Tensor` (D, D)
% right gauge with 2 legs,
% ordered left-right.
% Ar : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% left orthonormal
%
% Complexity
% ----------
% O(D^6) algorithm,
% diagonalizing (D^2, D^2) dmatrix
% find right fixed point
if nargin < 2 || isempty(r)
r = rightFixedPoint(A);
end
% decompose r = R * R
R = sqrtm(r);
% apply gauge R to A
Ar = contract(inv(R), [-1, 1], A, [1, -2, 2], R, [2, -3]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/rightOrthonormalize.m'.
[12]:
[L, Al] = leftOrthonormalize(A, l);
[R, Ar] = rightOrthonormalize(A, r);
assert(isapprox(R * R', r, 'RelTol', relTol), 'Right gauge doesn"t square to r')
assert(isapprox(L' * L, l, 'RelTol', relTol), 'Left gauge doesn"t sqaure to l')
assert(isapprox(contract(Ar, [-1, 1, 2], conj(Ar), [-2, 1, 2]), R.eye(R.codomain, R.domain), 'RelTol', 1e-10), 'Ar not in right-orthonormal form')
assert(isapprox(contract(Al, [1, 2, -2], conj(Al), [1, 2, -1]), L.eye(L.codomain, L.domain), 'RelTol', 1e-10), 'Al not in left-orthonormal form')
Finally, we can define a mixed gauge for the uniform MPS by choosing one site, the ‘center site’, and bringing all tensors to the left of it in the left-orthonormal form and all the tensors to the right of it in the right-orthonormal form. Defining a new tensor \(A_C\) on the center site, we obtain the form
By contrast, the original representation using the same tensor at every site is commonly referred to as the uniform gauge. The mixed gauge has an intuitive interpretation. Defining \(C = LR\), this tensor then implements the gauge transform that maps the left-orthonormal tensor to the right-orthonormal one, thereby defining the center-site tensor \(A_C\):
This relation is called the mixed gauge condition and allows us to freely move the center tensor \(A_C\) through the MPS, linking the left- and right orthonormal tensors.
Finally we may bring \(C\) into diagonal form by performing a singular value decomposition \(C = USV^\dagger\) and absorbing \(U\) and \(V^\dagger\) into the definition of \(A_L\) and \(A_R\) using the residual unitary gauge freedom
The mixed canonical form with a diagonal \(C\) now allows to straightforwardly write down a Schmidt decomposition of the state across an arbitrary bond in the chain
where the states \(\left | \Psi^i_L(A_L) \right \rangle\) and \(\left | \Psi^i_R(A_R) \right \rangle\) are orthogonal states on half the lattice. The diagonal elements \(C_i\) are exactly the Schmidt numbers of any bipartition of the MPS, and as such determine its bipartite entanglement entropy
[13]:
%%file mixedCanonical.m
function [Al, Ar, C, Ac] = mixedCanonical(A)
% Bring MPS tensor into mixed gauge, such that -Al-C- = -C-Ar- = Ac.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% Al : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% left orthonormal.
% Ac : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% center gauge.
% Ar : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% right orthonormal.
% C : :class:`Tensor` (D, D)
% Center gauge with 2 legs,
% ordered left-right,
% diagonal.
%
% Complexity
% ----------
% O(D^6) algorithm,
% diagonalization of (D^2, D^2) matrix
% Compute left and right orthonormal forms
[L, Al] = leftOrthonormalize(A);
[R, Ar] = rightOrthonormalize(A);
% center matrix C is matrix multiplication of L and R
C = L * R;
% singular value decomposition to diagonalize C
[U, S, V] = tsvd(C, 1, 2);
% absorb corresponding unitaries in Al and Ar
Al = contract(U', [-1, 1], Al, [1, -2, 2], U, [2, -3]);
Ar = contract(V, [-1, 1], Ar, [1, -2, 2], V', [2, -3]);
% normalize center matrix
C = S / norm(S);
% compute center MPS tensor
Ac = contract(Al, [-1, -2, 1], C, [1, -3]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/mixedCanonical.m'.
[14]:
%%file entanglementSpectrum.m
function [S, entropy] = entanglementSpectrum(A)
% Calculate the entanglement spectrum of an MPS.
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% S : double (D, 1)
% Singular values of center matrix,
% representing the entanglement spectrum
% entropy : double
% Entanglement entropy across a leg.
% go to mixed gauge
[~, ~, C, ~] = mixedCanonical(A);
% calculate entropy
S = diag(double(C));
entropy = -sum(S.^2 .* log(S.^2));
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/entanglementSpectrum.m'.
[15]:
% check that gauging is still correct
[Al, Ar, C, Ac] = mixedCanonical(A);
[S, entropy] = entanglementSpectrum(A);
assert(isapprox(contract(Ar, [-1, 1, 2], conj(Ar), [-2, 1, 2]), C.eye(C.codomain, C.domain), 'RelTol', 1e-10), 'Ar not in right-orthonormal form')
assert(isapprox(contract(Al, [1, 2, -2], conj(Al), [1, 2, -1]), C.eye(C.codomain, C.domain), 'RelTol', 1e-10), 'Al not in left-orthonormal form')
LHS = contract(Al, [-1, -2, 1], C, [1, -3]);
RHS = contract(C, [-1, 1], Ar, [1, -2, -3]);
assert(isapprox(LHS, RHS, 'RelTol', relTol) && isapprox(RHS, Ac, 'RelTol', relTol), 'Mixed gauge condition not satisfied!')
3 Truncation of a uniform MPS¶
The mixed canonical form also enables efficient truncatation an MPS. The sum in the above Schmidt decomposition can be truncated, giving rise to a new MPS that has a reduced bond dimension for that bond. This truncation is optimal in the sense that the norm between the original and the truncated MPS is maximized. To arrive at a translation invariant truncated MPS, we can truncate the columns of the absorbed isometries \(U\) and \(V^\dagger\) correspondingly, thereby transforming every tensor \(A_L\) or \(A_R\). The truncated MPS in the mixed gauge is then given by
We note that the resulting state based on this local truncation is not guaranteed to correspond to the MPS with a lower bond dimension that is globally optimal. This would require a variational optimization of the cost function.
[16]:
%%file truncateMPS.m
function [AlTilde, ArTilde, CTilde, AcTilde] = truncateMPS(A, Dtrunc)
% Truncate an MPS to a lower bond dimension.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
% Dtrunc : int
% lower bond dimension
%
% Returns
% -------
% AlTilde : :class:`Tensor` (Dtrunc, d, Dtrunc)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% left orthonormal.
% AcTilde : :class:`Tensor` (Dtrunc, d, Dtrunc)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% center gauge.
% ArTilde : :class:`Tensor` (Dtrunc, d, Dtrunc)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% right orthonormal.
% CTilde : :class:`Tensor` (Dtrunc, Dtrunc)
% Center gauge with 2 legs,
% ordered left-right,
% diagonal.
[Al, Ar, C, Ac] = mixedCanonical(A);
% perform SVD with truncation option
[U, S, V] = tsvd(C, 1, 2, 'TruncDim', Dtrunc);
% reabsorb unitaries
AlTilde = contract(U', [-1, 1], Al, [1, -2, 2], U, [2, -3]);
ArTilde = contract(V, [-1, 1], Ar, [1, -2, 2], V', [2, -3]);
% normalize center matrix
CTilde = S / norm(S);
% compute center MPS tensor
AcTilde = contract(AlTilde, [-1, -2, 1], CTilde, [1, -3]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/truncateMPS.m'.
[17]:
Dtrunc = 3;
[AlTilde, ArTilde, CTilde, AcTilde] = truncateMPS(A, Dtrunc);
assert(AlTilde.dims(1) == Dtrunc && AlTilde.dims(3) == Dtrunc, 'Something went wrong in truncating the MPS')
4 Algorithms for finding canonical forms¶
The success of using MPS for describing physical systems stems from the fact that they provide efficient approximations to a large class of physically relevant states. In one dimension, they have been shown to approximate low-energy states of gapped systems arbitrarily well at only a polynomial cost in the bond dimension \(D\). This means that in principle we can push MPS results for these systems to arbitrary precision as long as we increase the bond dimension enough. However, increasing the bond dimension comes at a numerical cost, as the complexity of any MPS algorithm scales with \(D\). As opposed to the naive routines given above, it is possible to ensure that the complexity of all MPS algorithms scales as \(O(D^3)\), so long as we are a bit careful when implementing them.
As a first example, we can refrain from explicitly contructing the matrices that are used in the eigenvalue problems, and instead pass a function that implements the action of the corresponding operator on a vector to the eigenvalue solver. We demonstrate this for the problem of normalizing an MPS, where instead of explicitly constructing the transfer matrix we now define a function handle which implements its action on the right and left fixed points using an optimal contraction sequence:
[18]:
%%file normalizeMPS.m
function Anew = normalizeMPS(A)
% Normalize an MPS tensor.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% Anew : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Complexity
% ----------
% O(D^3) algorithm,
% D^3 contraction for transfer matrix handle.
% define handle for transfer matrix acting on right vector
handleERight = @(v) contract(A, [-1, 2, 1], conj(A), [-2, 2, 3], v, [1, 3], 'Rank', [1, 1]);
% start from random initial guess for fixed point
r0 = similar([], A, 3, A, 3, 'Conj', [true, false], 'Rank', [1, 1]);
% calculate eigenvalue
lam = eigsolve(handleERight, r0);
Anew = A / sqrt(lam);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/normalizeMPS.m'.
[19]:
%%file leftFixedPoint.m
function l = leftFixedPoint(A)
% Find left fixed point.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% l : :class:`Tensor` (D, D)
% left fixed point with 2 legs,
% ordered bottom-top.
%
% Complexity
% ----------
% O(D^3) algorithm,
% D^3 contraction for transfer matrix handle.
% define handle for transfer matrix acting on left vector
handleELeft = @(v) contract(A, [1, 2, -2], conj(A), [3, 2, -1], v, [3, 1], 'Rank', [1, 1]);
% start from random initial guess for fixed point
l0 = similar([], A, 1, A, 1, 'Conj', [false, true], 'Rank', [1, 1]);
% calculate fixed point
[l, ~] = eigsolve(handleELeft, l0);
% make left fixed point hermitian and positive semidefinite explicitly
l = l * trace(l) / abs(trace(l)); % remove possible phase, actually forgot why I had to do this
l = (l + l') / 2; % force hermitian
l = l * sign(trace(l)); % force positive semidefinite
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/leftFixedPoint.m'.
[20]:
%%file rightFixedPoint.m
function r = rightFixedPoint(A)
% Find right fixed point.
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
% Returns
% -------
% r : :class:`Tensor` (D, D)
% right fixed point with 2 legs,
% ordered top-bottom.
% Complexity
% ----------
% O(D^3) algorithm,
% D^3 contraction for transfer matrix handle.
% define handle for transfer matrix acting on right vector
handleERight = @(v) contract(A, [-1, 2, 1], conj(A), [-2, 2, 3], v, [1, 3], 'Rank', [1, 1]);
% start from random initial guess for fixed point
r0 = similar([], A, 3, A, 3, 'Conj', [true, false], 'Rank', [1, 1]);
% calculate fixed point
[r, ~] = eigsolve(handleERight, r0);
% make right fixed point hermitian and positive semidefinite explicitly
r = r * trace(r) / abs(trace(r)); % remove possible phase, actually forgot why I had to do this
r = (r + r') / 2; % force hermitian
r = r * sign(trace(r)); % force positive semidefinite
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/rightFixedPoint.m'.
[21]:
%%file fixedPoints.m
function [l, r] = fixedPoints(A)
% Find normalized fixed points.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% l : :class:`Tensor` (D, D)
% left fixed point with 2 legs,
% ordered bottom-top.
% r : :class:`Tensor` (D, D)
% right fixed point with 2 legs,
% ordered top-bottom.
%
% Complexity
% ----------
% O(D^6) algorithm,
% diagonalizing (D^2, D^2) matrix
% find fixed points
l = leftFixedPoint(A);
r = rightFixedPoint(A);
% calculate trace and normalize
l = l / trace(l * r);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/fixedPoints.m'.
[22]:
A = createMPS(D, d);
A = normalizeMPS(A);
[l, r] = fixedPoints(A);
assert(isapprox(l, l', 'RelTol', relTol), 'left fixed point should be hermitian!')
assert(isapprox(r, r', 'RelTol', relTol), 'right fixed point should be hermitian!')
assert(isapprox(l, contract(A, [1, 2, -2], l, [3, 1], conj(A), [3, 2, -1]), 'RelTol', relTol), 'l should be a left fixed point!')
assert(isapprox(r, contract(A, [-1, 2, 1], r, [1, 3], conj(A), [-2, 2, 3]), 'RelTol', relTol), 'r should be a right fixed point!')
assert(isapprox(trace(l * r), 1, 'RelTol', relTol), 'Left and right fixed points should be trace normalized!') % TODO: fix once traces are a thing
We can similarly improve both the efficiency and accuracy of the routines bringing a given MPS into its mixed canonical form. While plugging in the more efficient ways of finding the left and right fixed point into the above mixedCanonical routine would reduce its complexity to \(O(D^3)\), this algorithm would still be suboptimal in terms of numerical accuracy. This arises from the fact that, while \(l\) and \(r\) are theoretically known to be positive hermitian matrices, at
least one of them will nevertheless have small eigenvalues, say of order \(\eta\), if the MPS is supposed to provide a good approximation to an actual state. In practice, \(l\) and \(r\) are determined using an iterative eigensolver and will only be accurate up to a specified tolerance \(\epsilon\). Upon taking the ‘square roots’ \(L\) and \(R\), the numerical precision will then decrease to \(\text{min}(\sqrt{\epsilon}, \epsilon/\sqrt{\eta})\). Furthermore, gauge
transforming \(A\) with \(L\) or \(R\) requires the potentially ill-conditioned inversions of \(L\) and \(R\), and will typically yield \(A_L\) and \(A_R\) which violate the orthonormalization condition in the same order \(\epsilon/\sqrt{\eta}\). We can circumvent both these probelems by resorting to so-called single-layer algorithms. These are algorithms that only work on the level of the MPS tensors in the ket layer, and never consider operations for which
contractions with the bra layer are needed. We now demonstrate such a single-layer algorithm for finding canonical forms.
Suppose we are given an MPS tensor \(A\), then from the above discussion we know that bringing it into left canonical form means finding a left-orthonormal tensor \(A_L\) and a matrix \(L\) such that \(L A=A_L L\). The idea is then to solve this equation iteratively, where in every iteration
we start from a matrix \(L^{i}\)
we construct the tensor \(L^{i}A\)
we take a QR decomposition to obtain \(A_L^{i+1} L^{i+1} = L^{i}A\), and
we take \(L^{i+1}\) to the next iteration
The QR decomposition is represented diagrammatically as
This iterative procedure is bound to converge to a fixed point for which \(L^{(i+1)}=L^{(i)}=L\) and \(A_L\) is left orthonormal by construction:
A similar procedure can be used to find a right-orthonormal tensor \(A_R\) and a matrix \(R\) such that \(A R = R A_R\). It is important to note that the convergence of this procedure relies on the fact that the QR decomposition is unique, which is not actually the case in general. However, it can be made unique by imposing that the diagonal elements of the triangular matrix \(R\) must be positive. This extra condition is imposed by calling the Tensor factorization methods
leftorth and rightorth with values 'qrpos' and 'rqpos' for the algorithm argument respectively.
[23]:
%%file rightOrthonormalize.m
function [R, Ar] = rightOrthonormalize(A, R0, tol, maxIter)
% Transform A to right-orthonormal gauge.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
% R0 : :class:`Tensor` (D, D), optional
% Right gauge matrix,
% initial guess.
% tol : float, optional
% convergence criterium,
% norm(R - Rnew) < tol.
% maxIter : int
% maximum amount of iterations.
%
% Returns
% -------
% R : :class:`Tensor` (D, D)
% right gauge with 2 legs,
% ordered left-right.
% Ar : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% right-orthonormal
arguments
A
R0 = []
tol = 1e-14
maxIter = 1e5
end
if isempty(R0)
R0 = similar([], A, 3, A, 3, 'Conj', [true, false], 'Rank', [1, 1]);
end
% Normalize R0
R0 = R0 / norm(R0);
% Initialize loop
i = 1;
[R, Ar] = rightorth(contract(A, [-1, -2, 1], R0, [1, -3]), 1, [2, 3], 'rqpos');
R = R / norm(R);
convergence = distance(R, R0);
% Decompose A*R until R converges
while convergence > tol
% calculate AR and decompose
[Rnew, Ar] = rightorth(contract(A, [-1, -2, 1], R, [1, -3]), 1, [2, 3], 'rqpos');
% normalize new R
Rnew = Rnew / norm(Rnew);
% calculate convergence criterium
convergence = distance(Rnew, R);
R = Rnew;
% check if iterations exceeds maxIter
if i > maxIter
warning('Right-orthonormalization has not converged!')
break
end
i = i + 1;
end
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/rightOrthonormalize.m'.
[24]:
%%file leftOrthonormalize.m
function [L, Al] = leftOrthonormalize(A, L0, tol, maxIter)
% Transform A to left-orthonormal gauge.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
% L0 : :class:`Tensor` (D, D), optional
% Left gauge matrix,
% initial guess.
% tol : float, optional
% convergence criterium,
% norm(R - Rnew) < tol.
% maxIter : int
% maximum amount of iterations.
%
% Returns
% -------
% L : :class:`Tensor` (D, D)
% left gauge with 2 legs,
% ordered left-right.
% Al : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% left-orthonormal
arguments
A
L0 = []
tol = 1e-14
maxIter = 1e5
end
if isempty(L0)
L0 = similar([], A, 1, A, 1, 'Conj', [false, true], 'Rank', [1, 1]);
end
% Normalize L0
L0 = L0 / norm(L0);
% Initialize loop
i = 1;
[Al, L] = leftorth(contract(L0, [-1, 1], A, [1, -2, -3]), [1, 2], 3, 'qrpos');
L = L / norm(L);
convergence = distance(L, L0);
% Decompose L*A until L converges
while convergence > tol
% calculate LA and decompose
[Al, Lnew] = leftorth(contract(L, [-1, 1], A, [1, -2, -3]), [1, 2], 3, 'qrpos');
% normalize new L
Lnew = Lnew / norm(Lnew);
% calculate convergence criterium
convergence = distance(Lnew, L);
L = Lnew;
% check if iterations exceeds maxIter
if i > maxIter
warning('Left-orthonormalization has not converged!')
break
end
i = i + 1;
end
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/leftOrthonormalize.m'.
[25]:
%%file mixedCanonical.m
function [Al, Ar, C, Ac] = mixedCanonical(A, L0, R0, tol, maxIter)
% Bring MPS tensor into mixed gauge, such that -Al-C- = -C-Ar- = Ac.
%
% Parameters
% ----------
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
%
% Returns
% -------
% Al : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% left orthonormal.
% Ac : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% center gauge.
% Ar : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% right orthonormal.
% C : :class:`Tensor` (D, D)
% Center gauge with 2 legs,
% ordered left-right,
% diagonal.
%
% Complexity
% ----------
% O(D^3) algorithm.
arguments
A
L0 = []
R0 = []
tol = 1e-14
maxIter = 1e5
end
tol = max(tol, 1e-14);
if isempty(L0)
L0 = similar([], A, 1, A, 1, 'Conj', [false, true], 'Rank', [1, 1]);
end
if isempty(R0)
R0 = similar([], A, 3, A, 3, 'Conj', [true, false], 'Rank', [1, 1]);
end
% Compute left and right orthonormal forms
[L, Al] = leftOrthonormalize(A, L0, tol, maxIter);
[R, Ar] = rightOrthonormalize(A, R0, tol, maxIter);
% center matrix C is matrix multiplication of L and R
C = L * R;
% singular value decomposition to diagonalize C
[U, S, V] = tsvd(C, 1, 2);
% absorb corresponding unitaries in Al and Ar
Al = contract(U', [-1, 1], Al, [1, -2, 2], U, [2, -3], 'Rank', Al.rank);
Ar = contract(V, [-1, 1], Ar, [1, -2, 2], V', [2, -3], 'Rank', Ar.rank);
% normalize center matrix
C = S / norm(S);
% compute center MPS tensor
Ac = contract(Al, [-1, -2, 1], C, [1, -3], 'Rank', Al.rank);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/mixedCanonical.m'.
[26]:
% check that gauging is still correct
[Al, Ar, C, Ac] = mixedCanonical(A);
assert(isapprox(contract(Ar, [-1, 1, 2], conj(Ar), [-2, 1, 2]), C.eye(C.codomain, C.domain), 'RelTol', relTol), 'Ar not in right-orthonormal form')
assert(isapprox(contract(Al, [1, 2, -2], conj(Al), [1, 2, -1]), C.eye(C.codomain, C.domain), 'RelTol', relTol), 'Al not in left-orthonormal form')
LHS = contract(Al, [-1, -2, 1], C, [1, -3]);
RHS = contract(C, [-1, 1], Ar, [1, -2, -3]);
assert(isapprox(LHS, RHS, 'RelTol', relTol) && isapprox(RHS, Ac, 'RelTol', relTol), 'Mixed gauge condition not satisfied!')
5 Computing expectation values¶
Now that we have seen the different ways to parametrize a given MPS, namely the uniform gauge and the mixed gauge, we wish to use these to compute expectation values of an extensive operator:
If we assume that each \(O_n\) acts on a single site and we are working with a properly normalized MPS, translation invariance dictates that the expectation value of \(O\) is given by the contraction
In the uniform gauge, we can use the fixed points of the transfer matrix to contract everything to the left and to the right of the operator, such that we are left with the contraction
In the mixed gauge, we can locate the center site where the operator is acting, and then contract everything to the left and right to the identity to arrive at the particularly simple expression
[27]:
%%file expVal1Uniform.m
function o = expVal1Uniform(O, A, l, r)
% Calculate the expectation value of a 1-site operator in uniform gauge.
% Parameters
% ----------
% O : :class:`Tensor` (d, d)
% single-site operator.
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
% l : :class:`Tensor` (D, D), optional
% left fixed point of transfermatrix,
% normalized.
% r : :class:`Tensor` (D, D), optional
% right fixed point of transfermatrix,
% normalized.
% Returns
% -------
% o : complex float
% expectation value of O.
arguments
O
A
l = []
r = []
end
% calculate fixed points if not given
if isempty(l) || isempty(r)
[l, r] = fixedPoints(A);
end
% contract expectation value network
o = contract(l, [4, 1], r, [3, 6], A, [1, 2, 3], conj(A), [4, 5, 6], O, [2, 5]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/expVal1Uniform.m'.
[28]:
%%file expVal1Mixed.m
function o = expVal1Mixed(O, Ac)
% Calculate the expectation value of a 1-site operator in mixed gauge.
%
% Parameters
% ----------
% O : :class:`Tensor` (d, d)
% single-site operator.
% Ac : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% center gauged.
%
% Returns
% -------
% o : complex float
% expectation value of O.
% contract expectation value network
o = contract(Ac, [2, 1, 3], conj(Ac), [2, 4, 3], O, [1, 4]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/expVal1Mixed.m'.
[29]:
O = Tensor.randnc([d, d]);
A = createMPS(D, d);
A = normalizeMPS(A);
[Al, Ar, C, Ac] = mixedCanonical(A);
expVal = expVal1Uniform(O, A);
expValMix = expVal1Mixed(O, Ac);
assert(isapprox(expVal, expValMix, 'RelTol', relTol), 'different gauges give different values?')
This procedure can be readily generalized to operators that act on multiple sites. In particular, a two-site operator such as a Hamiltonian term \(h\) can be evaluated as
[30]:
%%file expVal2Uniform.m
function o = expVal2Uniform(O, A, l, r)
% Calculate the expectation value of a 2-site operator in uniform gauge.
%
% Parameters
% ----------
% O : :class:`Tensor` (d, d, d, d)
% two-site operator,
% ordered topLeft-topRight-bottomLeft-bottomRight.
% A : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right.
% l : :class:`Tensor` (D, D), optional
% left fixed point of transfermatrix,
% normalized.
% r : :class:`Tensor` (D, D), optional
% right fixed point of transfermatrix,
% normalized.
%
% Returns
% -------
% o : complex float
% expectation value of O.
arguments
O
A
l = []
r = []
end
% calculate fixed points if not given
if isempty(l) || isempty(r)
[l, r] = fixedPoints(A);
end
% contract expectation value network
o = contract(l, [6, 1], r, [5, 10], A, [1, 2, 3], A, [3, 4, 5], conj(A), [6, 7, 8], conj(A), [8, 9, 10], O, [2, 4, 7, 9]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/expVal2Uniform.m'.
[31]:
%%file expVal2Mixed.m
function o = expVal2Mixed(O, Ac, Ar)
% Calculate the expectation value of a 2-site operator in mixed gauge.
% Parameters
% ----------
% O : :class:`Tensor` (d, d, d, d)
% two-site operator,
% ordered topLeft-topRight-bottomLeft-bottomRight.
% Ac : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% center gauged.
% Ar : :class:`Tensor` (D, d, D)
% MPS tensor with 3 legs,
% ordered left-bottom-right,
% right gauged.
% Returns
% -------
% o : complex float
% expectation value of O.
% contract expectation value network
o = contract(Ac, [4, 2, 1], Ar, [1, 3, 6], conj(Ac), [4, 5, 8], conj(Ar), [8, 7, 6], O, [2, 3, 5, 7]);
end
Created file '/home/leburgel/git/TensorTrack/docs/src/examples/uniformMps/expVal2Mixed.m'.
[32]:
O2 = Tensor.randnc([d, d, d, d]);
expVal = expVal2Uniform(O2, A);
expValGauge = expVal2Mixed(O2, Ac, Ar);
expValGauge2 = expVal2Mixed(O2, Al, Ac);
diff1 = abs(expVal - expValGauge);
diff2 = abs(expVal - expValGauge2);
assert(diff1 < relTol && diff2 < relTol, 'different gauges give different values?')
[ ]:
