Probabilistic Modeling with Matrix Product States

TL;DR

This paper introduces an efficient, gradient-free training algorithm for quantum circuit-inspired probabilistic models, demonstrating their usefulness for classical sequence modeling tasks like parity learning.

Contribution

It presents a novel, exactly solvable, density matrix renormalization group-based algorithm for training quantum circuit models on classical data.

Findings

Effective modeling of classical datasets using circuit-based models

Successful application to the parity learning problem

Supports the usefulness of quantum-inspired models in classical machine learning

Abstract

Inspired by the possibility that generative models based on quantum circuits can provide a useful inductive bias for sequence modeling tasks, we propose an efficient training algorithm for a subset of classically simulable quantum circuit models. The gradient-free algorithm, presented as a sequence of exactly solvable effective models, is a modification of the density matrix renormalization group procedure adapted for learning a probability distribution. The conclusion that circuit-based models offer a useful inductive bias for classical datasets is supported by experimental results on the parity learning problem.

Equations

$$P_{\psi }(e):=|⟨\psi ,e⟩{|}^2.$$

$${\psi }_{\pi }:=\sum _{x\in \text{pazocal}X}\sqrt{\pi (x)}|x⟩.$$

$$\pi (x)=P_{{\psi }_{\pi }}\left(|x⟩\right)\text{ for all x∈pazocalX.}$$

$${\psi }_{\ast }:=\underset{\psi \in \text{pazocal}M}{\mathrm{argmin}}\Vert \psi -{\psi }_{\pi }\Vert .$$

$$\hat{\pi }(x):=\frac{1}{n}\sum _{i=1}^n{\delta }_{X_i}(x).$$

$${\psi }_{\hat{\pi }}=\sum _{x\in \text{pazocal}X}\sqrt{\hat{\pi }(x)}|x⟩.$$

$$\hat{\psi }:=\underset{\psi \in \text{pazocal}M}{\mathrm{argmin}}\Vert \psi -{\psi }_{\hat{\pi }}\Vert .$$

$${\psi }_t\in \text{pazocal}{H}_{t+1}\text{ and }{\psi }_{t+1}=\underset{\{\psi \in \text{pazocal}{H}_{t+1}:\Vert \psi \Vert =1\}}{\mathrm{argmin}}\Vert {\psi }_{\hat{\pi }}-\psi \Vert .$$

$$\text{pazocal}{H}_{\mathrm{eff}}$$

$$P:=\alpha {\alpha }^{\ast }.$$

$$P({\psi }_{\hat{\pi }})$$

$$=\alpha {\alpha }^{\ast }\left(\sum _{x\in \text{pazocal}X}\sqrt{\hat{\pi }(x)}|x⟩\right)$$

$$=\alpha \left(\sum _{x\in \text{pazocal}X}\sqrt{\hat{\pi }(x)}{\alpha }^{\ast }(|x⟩)\right).$$

$$\sum _{x\in \text{pazocal}X}\sqrt{\hat{\pi }(x)}{\alpha }^{\ast }(|x⟩).$$

$$|e_{i_1}{e}_{i_2}\cdots e_{i_N}⟩:=|e_{i_1}⟩\otimes |e_{i_2}⟩\otimes \cdots \otimes |e_{i_N}⟩$$

$$\text{pazocal}{H}_{t+1}:={\alpha }_{t+1}(\text{pazocal}{H}_{\mathrm{eff},t+1})$$

$$\alpha :\text{pazocal}{H}_{\mathrm{eff}}\to \text{pazocal}H$$

$${\psi }_t=$$

$$\alpha$$

$$⟨\alpha (\varphi ),\alpha ({\varphi }^{\prime })⟩=$$

$${\alpha }^{\ast }$$

$$⟨\eta ,\alpha (\varphi )⟩=$$

$$\varphi =\sum _{x\in \text{pazocal}X}\sqrt{\hat{\pi }(x)}{\alpha }^{\ast }(|x⟩).$$

$${\alpha }^{\ast }\left(|x_i⟩\right)=$$

$$\frac{\varphi }{\Vert \varphi \Vert }=$$

$${\psi }_{t+1}=$$

$$\underset{\text{ Sweep 1}}{\underbrace{1,2,3,\dots ,N-1,N,N-1,\dots ,3,2,}}\underset{\text{ Sweep 2}}{\underbrace{1,2,\dots ,N-1,N,N-1,\dots ,2,}}1,2,\dots$$

$$\text{pazocal}{H}_{t+1}:={\alpha }_{t+1}(\text{pazocal}{H}_{\mathrm{eff},t+1})$$

$${\tilde{\psi }}_{t+1}:=\underset{\{\psi \in \text{pazocal}{H}_{t+1}:\Vert \psi \Vert =1\}}{\mathrm{argmin}}\Vert {\psi }_{\hat{\pi }}-\psi \Vert .$$

$${\tilde{\psi }}_{t+1}⇝{\psi }_{t+1}$$

$${\psi }_t=$$

Full text

1. Introduction

The possibility of exponential speedups for certain linear algebra operations has inspired a wave of research into quantum algorithms for machine learning purposes [3]. Many of these exponential speedups hinge on assumptions of fault tolerant quantum devices and efficient data preparation, which are unlikely to be realized in the near future. Focus has thus shifted to hybrid quantum-classical algorithms which involve optimizing the parameters of a variational quantum circuit to prepare a desired quantum state and have the potential to be implemented on near-term intermediate scale quantum devices [16].

Hybrid quantum-classical algorithms have been found to solve difficult eigenvalue problems [15] and to perform hard combinatorial optimization [7]. A number of recent works consider unsupervised learning within the hybrid quantum-classical framework [8, 5, 9, 12, 2, 6, 11].

In the context of machine learning, as emphasized in [16], it is less clear if variational hybrid quantum-classical algorithms offer advantages over existing purely classical algorithms. Density estimation, which attempts to learn a probability distribution from training data, has been suggested as an area to look for advantages [2] because a quantum advantage has been identified in the ability of quantum circuits to sample from certain probability distributions that are hard to sample classically [19]. In high-dimensional density estimation relevant to machine learning, expressive power is only part of the story and indeed algorithms in high-dimensional regime rely crucially on their inductive bias. Do the highly expressive probability distributions implied by quantum circuits offer a useful inductive bias for modeling high-dimensional classical data? We address this question in this paper.

We work within the confines of a classically tractable subset of quantum states modeled by tensor networks, which can be thought of as those states that can be prepared by shallow quantum circuits. Mathematically, tensor networks are a graphical calculus for describing interrelated matrix factorizations for which there exist polylogarithmic algorithms for a restricted set of linear algebra computations. We propose an unsupervised training algorithm for a generative model inspired by the density matrix renormalization group (DMRG) procedure. The training dynamics take place on the unit sphere of a Hilbert space, where in contrast to many variational methods, a state is modified in a sequence of deterministic steps that do not involve gradients. The efficient access to certain vector operations afforded by the tensor network ansatz allows us to implement our algorithm in a purely classical fashion.

We experimentally probe the inductive bias of the model by training on the dataset $P_{20}$ consisting of bitstrings of length $20$ having an even number of $1$ bits. The algorithm rapidly learns the uniform distribution on $P_{20}$ to high precision, indicating that the tensor network quantum circuit model provides a useful inductive bias for this classical dataset, which can be frustrating to learn for other models, such as RBMs.

In an effort to improve accessibility, we avoid the language of quantum-many body physics and quantum information and explain the algorithm and results in terms of elementary linear algebra and statistics. While this means some motivational material is omitted, we believe it sharpens the exposition. One exception is the visual language of tensor networks where the benefits of simplifying tensor contractions outweigh the costs of using elementary, but cumbersome, notation. We refer readers unfamiliar with tensor network notation to [20, 18, 4, 14] or to the many other surveys.

The organization of the paper is as follows. In section 2 we state the optimization problem at the population level and propose a finite-sample estimator. In sections 3 and 4 we describe an abstract discrete-time dynamical system evolving on the unit sphere of Hilbert space which optimizes our empirical objective by exactly solving an effective problem in a sequence of isometrically embedded Hilbert subspaces. In section 5 we provide a concrete realization of this dynamical system for a class of tensor networks called matrix product states. Section 6 outlines experiments demonstrating that the proposed iterative solver successfully the learns parity language using limited data.

Acknowledgements

The authors thank Tai-Danae Bradley, Giuseppe Carleo, Joseph Hirsh, Maxim Kontsevich, Miles Stoudenmire, Jackie Shadlen, and Yiannis Vlassopoulos for many helpful conversations.

2. The problem formulation

Recall that a unit vector $\psi$ in a finite-dimensional Hilbert space $\pazocal{H}$ defines a probability distribution $P_{\psi}$ on any orthonormal basis by setting the probability of each basis vector $e$ to be

[TABLE]

We refer to the probability distribution $P_{\psi}$ in Equation (1) as the Born distribution induced by $\psi$.

Let $\pi$ be a probability distribution on a finite set $\pazocal{X}$ and fix a field of scalars, either $\mathbb{R}$ or $\mathbb{C}$. Let $\pazocal{H}$ be the free vector space on the set $\pazocal{X}$. Use $|x\rangle$ to denote the vector in $\pazocal{H}$ corresponding to the element $x\in\pazocal{X}$. The space $\pazocal{H}$ has a natural inner product defined by declaring the vectors $\{|x\rangle:x\in\pazocal{X}\}$ to be an orthonormal basis.

Define a unit vector $\psi_{\pi}\in\pazocal{H}$ by

[TABLE]

Notice that $\psi_{\pi}$ realizes $\pi$ as a Born distribution:

[TABLE]

The formula for $\psi_{\pi}$ as written in Equation (2) involves perfect knowledge of $\pi$ and unrestricted access to the Hilbert space $\pazocal{H}$. This paper is concerned with situations when knowledge about $\pi$ is limited to a finite number of training examples, and $\psi$ is restricted to some tractable subset $\pazocal{M}$ of the unit sphere.

At the population level, the problem to be solved is to find the closest approximation $\psi_{\ast}$ to $\psi_{\pi}$ within $\pazocal{M}$,

[TABLE]

We assume access to a sequence $(X_{i})_{i=1}^{n}$ of samples drawn independently from $\pi$, giving rise to the associated empirical distribution

[TABLE]

It is natural to define the following estimator whose Born distribution coincides with the empirical distribution

[TABLE]

We are thus led to consider the following optimization problem.

Problem**.**

Given a sequence $\{X_{i}\}_{i=1}^{n}$ of i.i.d. samples drawn from $\pi$ and a subset $\pazocal{M}\subseteq\ \{\psi\in\pazocal{H}:\|\psi\|=1\}$ of the unit sphere in $\pazocal{H}$, find

[TABLE]

In the case where $\pazocal{X}$ consists of strings, the associated Hilbert space $\pazocal{H}$ has a dimension that is exponential in the string length. The model hypothesis class $\pazocal{M}$ should be chosen so that the induced Born distribution $P_{\widehat{\psi}}$ offers a useful inductive bias for modeling high-dimensional probability distributions over the space of sequences. We note, parenthetically, that the plug-in estimator $\|\psi-\psi_{\widehat{\pi}}\|$ is a biased estimator of the population objective $\|\psi-\psi_{\pi}\|$.

3. Outline of our approach to solving the problem

We present an algorithm that, given a fixed realization of data $(x_{1},\ldots,x_{n})\in\pazocal{X}^{n}$ and an initial state $\psi_{0}\in\pazocal{M}$, produces a deterministic sequence $\{\psi_{t}\}_{t\geq 0}$ of unit vectors in $\pazocal{M}$. The algorithm is a variation of the density matrix renormalization group (DMRG) procedure which we call exact single-site DMRG in which each step produces a vector closer to $\psi_{\widehat{\pi}}$. The sequence is defined inductively as follows: given $\psi_{t}$, the inductive step defines a subspace $\pazocal{H}_{t+1}$ of $\pazocal{H}$, which also contains $\psi_{t}$. Then $\psi_{t+1}$ is defined to be the vector in $\pazocal{H}_{t+1}$ closest to $\psi_{\widehat{\pi}}$. Inspired by ideas from the Renormalization Group we provide an analytic formula for $\psi_{t+1}$. The fact that the distance to the target vector $\psi_{\widehat{\pi}}$ decreases after each iteration follows as a simple consequence of the following facts

[TABLE]

See Figure 1.

4. Effective versions of the problem

Each proposal subspace $\pazocal{H}_{t}$ mentioned in the previous section will be defined as the image of an “effective” space. We begin with a general description of an effective space.

Let $\alpha:{\pazocal{H}_{\rm eff}}\to\pazocal{H}$ be an isometric embedding of a Hilbert space ${\pazocal{H}_{\rm eff}}$ into $\pazocal{H}$. We refer to ${\pazocal{H}_{\rm eff}}$ as the effective Hilbert space. The isometry $\alpha$ and its adjoint map $\alpha^{\ast}$ are summarized by the following diagram,

[TABLE]

The composition $\alpha^{*}\alpha=\operatorname{id}{{}_{\pazocal{H}_{\rm eff}}}$ is the identity on ${\pazocal{H}_{\rm eff}}$. The composition in the other order $\alpha\alpha^{*}$ is an orthogonal projection onto $\alpha({\pazocal{H}_{\rm eff}})$ which is a subspace of $\pazocal{H}$ isometrically isomorphic to ${\pazocal{H}_{\rm eff}}$. Call this orthogonal projection $P$

[TABLE]

The effective version of the problem formulated in Section 2 is to find the unit vector $\psi\in\alpha({\pazocal{H}_{\rm eff}})$ in the image of the effective Hilbert space that is closest to $\psi_{\widehat{\pi}}$. This effective problem is solved exactly in two simple steps. The first step is orthogonal projection: $P(\psi_{\widehat{\pi}})$ is the vector in $\alpha({\pazocal{H}_{\rm eff}})$ closest to $\psi_{\widehat{\pi}}$. The second step is to normalize $P(\psi_{\widehat{\pi}})$, which may not be a unit vector, to obtain the unit vector in $\alpha({\pazocal{H}_{\rm eff}})$ closest to $\psi_{\widehat{\pi}}$.

Therefore, the analytic solution to the effective problem is $P(\psi_{\widehat{\pi}})/\|P(\psi_{\widehat{\pi}})\|$ where

[TABLE]

In the exact single-site DMRG algorithm, the space $\alpha({\pazocal{H}_{\rm eff}})$ is contained within our model hypothesis class $\pazocal{M}$. We also offer a multi-site DMRG algorithm in the appendix. In this multi-site algorithm, the analytic solution to the effective problem in $\alpha({\pazocal{H}_{\rm eff}})$ does not lie in $\pazocal{M}$ so the solution to the effective problem needs to undergo an additional “model repair” step.

Before going on to the details of the algorithm, it might be helpful to look more closely at the solution to the effective problem. For each training example $x_{i}$, call the vector $\alpha^{\ast}(|x_{i}\rangle)\in{\pazocal{H}_{\rm eff}}$ an effective data point. Then, the argument of $\alpha$ in (10) becomes the weighted sum of effective data

[TABLE]

The effective data are not necessarily mutually orthogonal and so the vector in (11) will not be a unit vector. One may normalize to obtain a unit vector in ${\pazocal{H}_{\rm eff}}$ and then apply $\alpha$ to obtain the analytic solution to the effective problem. Normalizing in ${\pazocal{H}_{\rm eff}}$ and then applying $\alpha$ is the same as applying $\alpha$ and then normalizing in $\pazocal{H}$ since $\alpha$ is an isometry.

5. The exact single-site DMRG algorithm

Now specialize to the case that $\pi$ is a probability distribution on a set $\pazocal{X}$ of sequences. Suppose that $\pazocal{X}=A^{N}$ consists of sequences of length $N$ in fixed alphabet $A=\{e_{1},\ldots,e_{d}\}$. The Hilbert space $\pazocal{H}$, defined as the free Hilbert space on $\pazocal{X}$, has a natural tensor product structure $V^{\otimes N}$ where $V$ is the free Hilbert space on the alphabet $A$. We refer to $V$ as the site space. So in this situation, the vectors $\{|e_{1}\rangle,\ldots,|e_{d}\rangle\}$ are an orthonormal basis for the $d$-dimensional site space $V$ and the vectors

[TABLE]

are an orthonormal basis for the $d^{N}$ dimensional space $\pazocal{H}=V^{\otimes N}$. We choose as model hypothesis class the subset $\pazocal{M}\subseteq\pazocal{H}$ consisting of normalized elements in $\pazocal{H}$ that have a low rank matrix product state (MPS) factorization. Vectors in this model hypothesis class have efficient representations, even in cases where the Hilbert space $\pazocal{H}$ is of exponentially high dimension. For simplicity of presentation, we consider matrix product states with a single fixed bond space $W$, although everything that follows could be adapted to work with tensor networks without loops having arbitrary bond spaces.

The exact single-site DMRG algorithm begins with an initial vector $\psi_{0}\in\pazocal{M}$ and produces $\psi_{1},\psi_{2},\ldots$ inductively by solving an effective problem

[TABLE]

which we now describe. Let us drop the subscript $t+1$ from the isometry $\alpha_{t+1}$ and the effective Hilbert space $\pazocal{H}_{\mathrm{eff},t+1}$ in the relevant effective problem—just be aware that the embedding

[TABLE]

will change from step to step. The map $\alpha$ is defined using an MPS factorization of $\psi_{t}$ in mixed canonical form relative to a fixed site which varies at each step according to a predetermined schedule. For the purposes of illustration, the third site is the fixed site in the pictures below.

[TABLE]

The effective space is ${\pazocal{H}_{\rm eff}}=W\otimes V\otimes W$ and the isometric embedding $\alpha:W\otimes V\otimes W\to V^{\otimes N}$ is defined for any $\phi\in W\otimes V\otimes W$ by replacing the tensor at the fixed site of $\psi_{t}$ with $\phi$:

[TABLE]

To see that $\alpha$ is an isometry, use the gauge condition that the MPS factorization of $\psi_{t}$ is in mixed canonical form relative to the fixed site, as illustrated below:

[TABLE]

The adjoint map $\alpha^{*}:V^{\otimes N}\to W\otimes V\otimes W$ has a clean pictorial depiction as well.

[TABLE]

To see that $\alpha^{*}$ as pictured above is, in fact, the adjoint of $\alpha$, note that for any $\eta\in\pazocal{H}$ and any $\phi\in{\pazocal{H}_{\rm eff}}$, both $\langle\eta,\alpha(\phi)\rangle$ and $\langle\alpha^{*}(\eta),\phi\rangle$ result in the same tensor contraction:

[TABLE]

In the picture above, begin with the blue tensors. Contracting with the yellow tensor gives $\alpha(\phi)$ and then contracting with the red tensor gives $\langle\eta,\alpha(\phi)\rangle$. On the other hand, first contracting with the red tensor yields $\alpha^{*}(\eta)$ resulting in $\langle\alpha^{*}(\eta),\phi\rangle$ after contracting with the yellow tensor.

Now, Equation (10) describes an analytic solution for the vector in $\pazocal{H}_{t+1}:=\alpha(W\otimes V\otimes W)$ closest to $\psi_{\widehat{\pi}}$. Namely, $\alpha(\phi/\|\phi\|)$ where

[TABLE]

For each sample $|x_{i}\rangle=|e_{i_{1}}e_{i_{2}}\cdots e_{i_{N}}\rangle$, the effective data point $\alpha^{\ast}(|x_{i}\rangle)\in V\otimes W\otimes V$ is given by the contraction

[TABLE]

Once the effective form $\alpha^{\ast}(|x\rangle)$ of each distinct training example $|x\rangle$ has been computed, weighted by $\sqrt{\widehat{\pi}(x)}$, summed, and normalized, one obtains an expression for the unit vector $\phi/\|\phi\|\in W\otimes V\otimes W$, depicted as follows,

[TABLE]

Finally, apply the map $\alpha$ to get $\psi_{t+1}$:

[TABLE]

To complete the description of the exact single-site DMRG algorithm, we need to choose a schedule in which to update the tensors. We use the following schedule, organized into back-and-forth sweeps, for the fixed site at each step

[TABLE]

A schedule that proceeds by moving the fixed site one position at a time allows us to take advantage of two efficiencies resulting in an algorithm that is linear in both the number of training examples $n$ and the number of sites $N$. One efficiency is that most of the calculations of the effective data in Equation (21) used to compute $\psi_{t+1}$ can be reused when computing $\psi_{t+2}$. The second efficiency is that when inserting the updated tensor in Equation (22), it can be done so that the resulting MPS factorization of $\psi_{t+1}$ as pictured in Equation (23) will be in mixed canonical form relative to a site adjacent to the updated tensor, which avoids a costly gauge fixing step.

6. Experiments

This section considers the problem of unsupervised learning of the parity language $P_{N}$, which consists of bitstrings of length $N$ containing an even number of 1 bits. Of the total $2^{N-1}$ such bitstrings, we reserved random disjoint subsets of size $2\%$ for training, cross-validation and testing purposes. A NLL of $N-1$ corresponds to the entropy of the uniform distribution on $P_{N}$. If the model memorizes the training set, it will assign to it a negative-log-likelihood (NLL) of $N-1+\log_{2}(0.02)$ corresponding to the entropy of the uniform distribution on the training data. A NLL of $N$ corresponds to the entropy of the uniform distribution on all bitstrings of length $N$. The measure of generalization performance is the gap $\epsilon$ between the NLL of the training and testing data. We performed exact single-site DMRG over the real number field using the $P_{20}$ dataset for different choices of bond dimension. Training was terminated according to an early stopping criterion as determined by distance between the MPS state and the state of the cross-validation sample. The NLL as a function of bond dimension reported in Fig. 2 displays the expected bias-variance tradeoff, with optimal model complexity occurring at bond dimension $3$ with corresponding generalization gap $\epsilon=0.0237$.

7. Discussion

A number of recent works have explored the parity dataset using restricted Boltzmann machines (RBMs) and found it to be difficult to learn, even in experiments that train using the entire dataset [17, 13]. Recall that an RBM is a universal approximator of distributions on $\{0,1\}^{N}$, given sufficiently many hidden units. Ref. [13] proved that any probability distribution on $\{0,1\}^{N}$ can be approximated within $\epsilon$ in KL-divergence by an RBM with $m\geq 2^{(N-1)(1-\epsilon)+0.1}$ hidden units. For $P_{20}$ this bound works out to be about $4\times 10^{5}$ hidden nodes. It would be interesting to know whether it could be learned with significantly fewer.

It is not difficult to train a feedforward neural network to classify bitstrings by parity using labelled data, but we do not know if there are unsupervised generative neural models that do well learning $P_{N}$. It is reasonable to expect that recurrent models whose training involve conditional probabilities $\pi(x_{1},\ldots,x_{k}|x_{k+1},\ldots,x_{N})$ might be frustrated by $P_{N}$ since the conditional distributions contain no information: any bitstring of length less than $N$ has the same number of completions in $P_{N}$ as not in $P_{N}$.

The reader may be interested in [1, 10] where quantum models are used to learn classical data. Those works considered quantum Boltzman machines which were shown to learn the distribution more effectively than their classical counterparts using the same dataset. The complexity of classically simulating a QBM scales exponentially with the number of sites in contrast to the tensor network algorithms presented here, which scale linearly in the number of sites.

Appendix A Multi-site DMRG

For completeness we now describe a related multi-site DMRG algorithm. The model class $\pazocal{M}$ now consists of normalized vectors with matrix product factorizations, with possibly different bond spaces having dimension less than a fixed upper bound. The algorithm begins with an initial vector $\psi_{0}\in\pazocal{M}$ and produces $\psi_{1},\psi_{2},\ldots$ inductively. The inductive step is similar in that we solve an effective problem in the image of an effective Hilbert space

[TABLE]

to find the unit vector in $\pazocal{H}_{t+1}$ that is closest to the target state $\psi_{\widehat{\pi}}$, which we now denote with a tilde:

[TABLE]

In multi-site DMRG, as opposed to single-site DMRG, the image of the effective space $\pazocal{H}_{t+1}$ is not contained in the MPS model hypothesis class $\pazocal{M}$. So, the solution $\widetilde{\psi}_{t+1}$ to the effective problem must undergo a “model repair” step

[TABLE]

to produce a vector $\psi_{t+1}\in\pazocal{M}$. In summary:

• •

Use $\psi_{t}$ to define an isometric embedding $\alpha_{t+1}:{\pazocal{H}_{\rm eff}}\to\pazocal{H}$ with $\psi_{t}\in\pazocal{H}_{t+1}:=\alpha_{t+1}({\pazocal{H}_{\rm eff}}).$

• •

Let $\widetilde{\psi}_{t+1}$ be the unit vector in $\pazocal{H}_{t+1}$ closest to $\psi_{\widehat{\pi}}$.

• •

Perform a model repair of $\widetilde{\psi}_{t+1}$ to obtain a vector ${\psi}_{t+1}\in\pazocal{M}.$ There are multiple ways to do the model repair.

In order to define the effective problem in the inductive step of multi-site DMRG, one uses an MPS factorization of $\psi_{t}$ in mixed canonical gauge relative to an interval of $r$-sites. In the picture below, the interval consists of the two sites $3$ and $4$.

[TABLE]

The effective Hilbert space ${\pazocal{H}_{\rm eff}}=W_{L}\otimes V^{\otimes r}\otimes W_{R}$ where $W_{L}$ and $W_{R}$ are the bond spaces to the left and right of the fixed interval of sites, and $r$ is the length of the chosen interval. The map $\alpha:W_{L}\otimes V^{\otimes r}\otimes W_{R}\to V^{\otimes n}$ is given by replacing the interval of sites and contracting

[TABLE]

The map $\alpha$ and its adjoint $\alpha^{\ast}$ are described by, and have properties proved by, pictures completely analogous to those detailed for single-site DMRG in Section 5. The effective problem is also solved the same way. What is not the same is that the vector in $\pazocal{H}_{t+1}=\alpha(W_{L}\otimes V^{\otimes r}\otimes W_{R})$ which solves the effective problem is outside of the model class $\pazocal{M}$ and so one performs a model repair step $\widetilde{\psi}_{t+1}\leadsto\psi_{t+1}$, pictured graphically in ${\pazocal{H}_{\rm eff}}$ by:

[TABLE]

One way to perform the model repair is to choose

[TABLE]

but the flexibility of the model repair step allows for other possibilities. One can use the model repair to implement a dynamic tradeoff between proximity to $\widetilde{\psi}_{t+1}$ and other constraints of interest, such as bond dimension. Many of these implementations have good algorithms arising from singular value decompositions manageable in the effective Hilbert space. Let use denote such a model repair choice as $\psi_{t+1}^{\rm SVD}$. Be aware that if $\psi_{t+1}^{\rm SVD}$ is the vector in $\pazocal{M}\cap\pazocal{H}_{t+1}$ nearest to $\widetilde{\psi}_{t+1}$ as in Equation (31), there is no guarantee that $\psi_{t+1}^{\rm SVD}$ will be nearer to $\psi_{\widehat{\pi}}$ than the previous iterate. In fact, we have experimentally observed the sequence obtained by this kind of model repair to move away from $\psi_{\widehat{\pi}}$. See Figure 3 for an illustration of this possibility.

One might hope to improve the model repair step, say by pre-conditioning the singular value decomposition in a way that is knowledgeable about the target $\psi_{\widehat{\pi}}$. For the experiments reported in this paper, single-site DMRG consistently outperformed multi-site DMRG for several choices of model repair step, and we include multi-site DMRG only for pedagogical reasons. The adaptability of the bond dimension afforded by the multi-site DMRG algorithm could provide benefits that outweigh the challenges of good model repair in some situations.