CBMM Memo No. 047 April 12, 2016

Bridging the Gaps Between Residual Learning, Recurrent Neural Networks and
Visual Cortex
by
Qianli Liao and Tomaso Poggio
Center for Brains, Minds and Machines, McGovern Institute, MIT

Abstract: We discuss relations between Residual Networks (ResNet), Recurrent
Neural Networks (RNNs) and the primate visual cortex. We begin with the
observation that a shallow RNN is exactly equivalent to a very deep ResNet
with weight sharing among the layers. A direct implementation of such a RNN,
although having orders of magnitude fewer parameters, leads to a performance
similar to the corresponding ResNet. We propose 1) a generalization of both
RNN and ResNet architectures and 2) the conjecture that a class of moderately
deep RNNs is a biologically-plausible model of the ventral stream in visual
cortex. We demonstrate the effectiveness of the architectures by testing them
on the CIFAR-10 dataset.

This work was supported by the Center for Brains, Minds and Machines (CBMM),
funded by NSF STC award CCF - 1231216.

##  1 Introduction

Residual learning [8], a novel deep learning scheme characterized by ultra-
deep architectures has recently achieved state-of-the-art performance on
several popular vision benchmarks. The most recent incarnation of this idea
[10] with hundreds of layers demonstrate consistent performance improvement
over shallower networks. The 3.57% top-5 error achieved by residual networks
on the ImageNet test set arguably rivals human performance.

Because of recent claims [33] that networks of the AlexNet[16] type
successfully predict properties of neurons in visual cortex, one natural
question arises: how similar is an ultra-deep residual network to the primate
cortex? A notable difference is the depth. While a residual network has as
many as 1202 layers[8], biological systems seem to have two orders of
magnitude less, if we make the customary assumption that a layer in the NN
architecture corresponds to a cortical area. In fact, there are about half a
dozen areas in the ventral stream of visual cortex from the retina to the
Inferior Temporal cortex. Notice that it takes in the order of 10ms for neural
activity to propagate from one area to another one (remember that spiking
activity of cortical neurons is usually well below 100 Hz). The evolutionary
advantage of having fewer layers is apparent: it supports rapid (100msec from
image onset to meaningful information in IT neural population) visual
recognition, which is a key ability of human and non-human primates [31, 28].

It is intriguingly possible to account for this discrepancy by taking into
account recurrent connections within each visual area. Areas in visual cortex
comprise six different layers with lateral and feedback connections [17],
which are believed to mediate some attentional effects [2, 17, 14, 25, 12] and
even learning (such as backpropagation [20]). “Unrolling” in time the
recurrent computations carried out by the visual cortex provides an equivalent
“ultra-deep” feedforward network, which might represent a more appropriate
comparison with the state-of-the-art computer vision models.

In addition, we conjecture that the effectiveness of recent “ultra-deep”
neural networks primarily come from the fact they can efficiently model the
recurrent computations that are required by the recognition task. We show
compelling evidences for this conjecture by demonstrating that 1. a deep
residual network is formally equivalent to a shallow RNN; 2. such a RNN with
weight sharing, thus with orders of magnitude less parameters (depending on
the unrolling depth), can retain most of the performance of the corresponding
deep residual network.

Furthermore, we generalize such a RNN into a class of models that are more
biologically-plausible models of cortex and show their effectiveness on
CIFAR-10.

##  2 Equivalence of ResNet and RNN

###  2.1 Intuition

We discuss here a very simple observation: a Residual Network (ResNet)
approximates a specific, standard Recurrent Neural Network (RNN) implementing
the discrete dynamical system described by

$h_{t}=K\circ(h_{t-1})+h_{t-1}$ (1)

where $h_{t}$ is the activity of the neural layer at time $t$ and $K$ is a
nonlinear operator. Such a dynamical systems corresponds to the feedback
system of Figure 1 (B). Figure 1 (A) shows that unrolling in (discrete) time
the feedback system gives a deep residual network with the same (that is,
shared) weights among the layers. The number of layers in the unrolled network
corresponds to the discrete time iterations of the dynamical system. The
identity shortcut mapping that characterizes residual learning appears in the
figure.

Figure 1: A formal equivalence of a ResNet (A) with weight sharing and a RNN
(B). I is the identity operator. K is an operator denoting the nonlinear
transformation called $f$ in the main text. $x_{t}$ is the value of the input
at time $t$. $\delta_{t}$ is a Kronecker delta function.

Thus, ResNets with shared weights can be reformulated into the form of a
recurrent system. In section 5.2) we show experimentally that a ResNet with
shared weights retains most of its performance (on CIFAR-10).

A comparison of a plain RNN and ResNet with shared weights is in the Appendix
Figure 11.

###  2.2 Formulation in terms of Dynamical Systems (and Feedback)

We frame recurrent and residual neural networks in the language of dynamical
systems. We consider here dynamical systems in discrete time, though most of
the definitions carry over to continuous time. A neural network (that we
assume for simplicity to have a single layer with $n$ neurons) can be a
dynamical system with a dynamics defined as

$h_{t+1}=f(h_{t};w_{t})+x_{t}$ (2)

where $h_{t}\in\mathbb{R}^{n}$ is the activity of the $n$ neurons in the layer
at time $t$ and $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a continuous, bounded
function parametrized by the vector of weights $w_{t}$. In a typical neural
network, $f$ is synthesized by the following relation between the activity
$y_{t}$ of a single neuron and its inputs $x_{t-1}$:

$y_{t}=\sigma(\left\langle{w},{x_{t-1}}\right\rangle+b),$ (3)

where $\sigma$ is a nonlinear function such as the linear rectifier
$\sigma(\cdot)=|\cdot|_{+}$.

A standard classification of dynamical systems defines the system as

  1. 1.

homogeneous if $x_{t}=0,\,\,\forall t>0$ (alternatively the equation reads as
$h_{t+1}=f(h_{t};w_{t})$ with the inital condition $h_{0}=x_{0}$)

  2. 2.

time invariant if $w_{t}=w$.

Residual networks with weight sharing thus correspond to homogeneous, time-
invariant systems which in turn correspond to a feedback system (see Figure 1)
with an input which is non-zero only at time $t=0$
($x_{t=0}=x_{0},x_{t}=0\,\,\forall t>0$) and with $f(z)=(K+I)\circ z$:

$h_{n}=f(h_{t};w_{t})=(K+I)^{n}\circ x_{0}$ (4)

“Normal” residual networks correspond to homogeneous, time-variant systems. An
analysis of the corresponding inhomogeneous, time-invariant system is provided
in the Appendix.

##  3 A Generalized RNN for Multi-stage Fully Recurrent Processing

As shown in the previous section, the recurrent form of a ResNet is actually
shallow (if we ignore the possible depth of the operator $K$). In this
section, we generalize it into a moderately deep RNN that reflects the multi-
stage processing in the primate visual cortex.

###  3.1 Multi-state Graph

We propose a general formulation that can capture the computations performed
by a multi-stage processing hierarchy with full recurrent connections. Such a
hierarchy can be characterized by a directed (cyclic) graph G with vertices V
and edges E.

$G=\\{V,E\\}$ (5)

where vertices V is a set contains all the processing stages (i.e., we also
call them states). Take the ventral stream of visual cortex for example,
$V=\\{LGN,V1,V2,V4,IT\\}$. Note that retina is not listed since there is no
known feedback from primate cortex to the retina. The edges $E$ are a set that
contains all the connections (i.e., transition functions) between all
vertices/states, e.g., V1-V2, V1-V4, V2-IT, etc. One example of such a graph
is in Figure 2 (A).

Figure 2: Modeling the ventral stream of visual cortex using a multi-state
fully recurrent neural network

###  3.2 Pre-net and Post-net

The multi-state fully recurrent system does not have to receive raw inputs.
Rather, a (deep) neural network can serve as a preprocesser. We call the
preprocesser a “pre-net” as shown in Figure 2 (B). On the other hand, one also
needs a “post-net” as a postprocessor and provide supervisory signals to the
recurrent system and the pre-net. The pre-net, recurrent system and post-net
are trained in an end-to-end fashion with backpropagation.

For most models in this paper, unless stated otherwise, the pre-net is a
simple 3x3 convolutional layer and the post-net is a pipeline of a batch
normalization, a ReLU, a global average pooling and a fully connected layer
(or a 1x1 convolution, we use these terms interchangeably).

Take primate visual system for instance, the retina is a part of the “pre-
net”. It does not receive any feedback from the cortex and thus can be
separated from the recurrent system for simplicity. In Section 5.3.3, we also
tried 3 layers of 3x3 convolutions as an pre-net, which might be more similar
to a retina, and observed slightly better performance.

###  3.3 Transition Matrix

The set of edges $E$ can be represented as a 2-D matrix where each element
(i,j) represents the transition function from state $i$ to state $j$.

One can also extend the representation of $E$ to a 3-D matrix, where the third
dimension is time and each element (i,j,t) represents the transition function
from state $i$ to state $j$ at time $t$. In this formulation, the transition
functions can vary over time (e.g., being blocked from time $t_{1}$ to time
$t_{2}$, etc.). The increased expressive power of this formulation allows us
to design a system where multiple locally recurrent systems are connected
sequentially: a downstream recurrent system only receives inputs when its
upstream recurrent system finishes, similar to recurrent convolutional neural
networks (e.g., [19]). This system with non-shared weights can also represent
exactly the state-of-the-art ResNet (see Figure 3).

Nevertheless, time-dependent dynamical systems, that is recurrent networks of
real neurons and synapses, offer interesting questions about the flexibility
in controlling their time dependent parameters.

Example transition matrices used in this paper are shown in Figure 4.

When there are multiple transition functions to a state, their outputs are
summed together.

Figure 3: We show two types of ResNet and corresponding RNNs. (A) and (C):
Single state ResNet. The spacial and featural sizes are fixed. (B) and (D)
3-state ResNet. The spacial size is reduced by a factor of 2 and the featural
size is doubled at time n and 2n, where n is a meta parameter. This type
corresponds to the ones proposed by He et. al. [8]. Figure 4: The transition
matrices used in the paper. “BN” denotes Batch Normalization and “Conv”
denotes convolution. Deconvolution layer (denoted by “Deconv”) is [34] used as
a transition function from a spacially small state to a spacially large one.
BRCx2/BRDx2 denotes a BN-ReLU-Conv/Deconv-BN-ReLU-Conv/Deconv pipeline
(similar to a residual module [10]). There is always a 2x2
subsampling/upsampling between nearby states (e.g., V1/h1: 32x32, V2/h2:
16x16, V4/h3:8x8, IT:4x4). Stride 2 (convolution) or upsampling 2
(deconvolution) is used in transition functions to match the spacial sizes of
input and output states. The intermediate feature sizes of transition function
BRCx2/BRDx2 or BRCx3/BRDx3 are chosen to be the average feature size of input
and output states. “+I” denotes a identity shortcut mapping. The design of
transition functions could be an interesting topic for future research.

###  3.4 Shared vs. Non-shared Weights

Weight sharing is described at the level of an unrolled network. Thus, it is
possible to have unshared weights with a 2D transition matrix — even if the
transitions are stable over time, their weights could be time-variant.

Given an unrolled network, a weight sharing configuration can be described as
a set $S$, whose element is a set of tied weights
$s=\\{W_{i_{1},j_{1},t_{1}},...,W_{i_{m},j_{m},t_{m}}\\}$, where
$W_{i_{m},j_{m},t_{m}}$ denotes the weight of the transition functions from
state $i_{m}$ to $j_{m}$ at time $t_{m}$. This requires: 1. all weights
$W_{i_{m},j_{m},t_{m}}\in s$ have the same initial values. 2. the actual
gradients used for updating each element of $s$ is the sum of the gradients of
all elements in $s$:

$\forall W\in s,\,\,(\frac{\partial E}{\partial W})_{used}=\sum_{W^{\prime}\in
s}(\frac{\partial E}{\partial W^{\prime}})_{original}$ (6)

where $E$ is the training objective.

For RNNs, weights are usually shared across time, but one could unshare the
weights, share across states or perform more complicated sharing using this
framework.

###  3.5 Notations: Unrolling Depth vs. Readout Time

The meaning of “unrolling depth” may vary in different RNN models since
“unrolling” a cyclic graph is not well defined. In this paper, we adopt a
biologically-plausible definition: we simulate the time after the onset of the
visual stimuli assuming each transition function takes constant time 1. We use
the term “readout time” to refer to the time the post-net reads the data from
the last state.

This definition in principle allows one to have quantitive comparisons with
biological systems. e.g., for a model with readout time $t$ in this paper, the
wall clock time can be estimated to be $20t$ to $50t$ ms, considering the
latency of a single layer of biological neurons.

Regarding the initial values, at $t=0$ all states are empty except that the
first state has some data received from the pre-net. We only start simulate a
transition function when its input state is populated.

###  3.6 Sequential vs. Static Inputs/Outputs

As an RNN, our model supports sequential data processing and in principle all
other tasks supported by traditional RNNs. See Figure 5 for illustrations.
However, if there is a batch normalization in the model, we have to use “time-
specific normalization” described in Section 3.7, which might not be feasible
for some tasks.

Figure 5: Our model supports sequential inputs/outputs, which includes one-to-
one, many-to-one, one-to-many and many-to-many input-output mappings.

###  3.7 Batch Normalizations for RNNs

As an additional observation, we found that it generally hurts performance
when the normalization statistics (e.g., average, standard deviation,
learnable scaling and shifting parameters) in batch normalization are shared
across time. This may be consistent with the observations from [18].

However, good performance is restored if we apply a procedure we call a “time-
specific normalization”: mean and standard deviation are calculated
independently for every $t$ (using training set). The learnable scaling and
shifting parameters should be time-specific. But in most models we do not use
the learnable parameters of BN since they tend not to affect the performance
much.

We expect this procedure to benefit other RNNs trained with batch
normalization. However, to use this procedure, one needs to have a initial
$t=0$ and enumerate all possible $t$s. This is feasible for visual processing
but needs modifications for other tasks.

##  4 Related Work

Deep Recurrent Neural Networks: Our final model is deep and similar to a
stacked RNN [27, 5, 7] with several main differences: 1. our model has
feedback transitions between hidden layers and self-transition from each
hidden layer to itself. 2. our model has identity shortcut mappings inspired
by residual learning. 3. our transition functions are deep and convolutional.

As suggested by [23], the term depth in RNN could also refer to input-to-
hidden, hidden-to-hidden or hidden-to-output connections. Our model is deep in
all of these senses. See Section 3.2.

Recursive Neural Networks and Convolutional Recurrent Neural Networks:  When
unfolding RNN into a feedforward network, the weights of many layers are tied.
This is reminiscent of Recursive Neural Networks (Recursive NN), first
proposed by [29]. Recursive NN are characterized by applying same operations
recursively on a structure. The convolutional version was first studied by
[4]. Subsequent related work includes [24] and [19]. One characteristic
distinguishes our model and residual learning from Recursive NN and
convolutional recurrent NN is whether there are identity shortcut mappings.
This discrepancy seems to account for the superior performance of residual
learning and of our model over the latters.

A recent report [3] we became aware of after we finished this work discusses
the idea of imitating cortical feedback by introducing loops into neural
networks.

A Highway Network [30] is a feedforward network inspired by Long Short Term
Memory [11] featuring more general shortcut mappings (instead of hardwired
identity mappings used by ResNet).

##  5 Experiments

###  5.1 Dataset and training details

We test all models on the standard CIFAR-10 [15] dataset. All images are 32x32
pixels with color. Data augmentation is performed in the same way as [8].

Momentum was used with hyperparameter 0.9. Experiments were run for 60 epochs
with batchsize 64 unless stated otherwise. The learning rates are 0.01 for the
first 40 epochs, 0.001 for epoch 41 to 50 and 0.0001 for the last 10 epochs.
All experiments used the cross-entropy loss function and softmax for
classification. Batch Normalization (BN) [13] is used for all experiments. But
the learnable scaling and shifting parameters are not used (except for the
last BN layer in the post-net). Network weights were initialized with the
method described in [9]. Although we do not expect the initialization to
matter as long as batch normalization is used. The implementations are based
on MatConvNet[32].

###  5.2 Experiment A: ResNet with shared weights

####  5.2.1 Sharing Weights Across Time

We conjecture that the effectiveness of ResNet mainly comes from the fact that
it efficiently models the recurrent computations required by the recognition
task. If this is the case, one should be able to reinterpret ResNet as a RNN
with weight sharing and achieve comparable performance to the original
version. We demonstrate various incarnations of this idea and show it is
indeed the case.

We tested the 1-state and 3-state ResNets described in Figure 3 and 4. The
results are shown in Figure 6.

Figure 6: All models are robust to sharing weights across time. This supports
our conjecture that deep networks can be well approximated by
shallow/moderately-deep RNNs. The transition matrices of all models are shown
in Figure 4. “#Param” denotes the number of parameters. The 1-state ResNet has
a single state of size 32x32x64 (height $\times$ width $\times$ #features). It
was trained and tested with readout time t=10. The 3-state ResNet has 3 states
of size 32x32x16, 16x16x32 and 8x8x64 – there is a transition (via a simple
convolution) at time 4 and 8 — each state has a self-transition unrolled 3
times. The 2-state fully recurrent NN has 2 states of the same size: 32x32x64.
It was trained and tested with readout time t=10 (same as 1-state ResNet). It
is a generalization of and directly comparable with 1-state ResNet, showing
the benefit of having more states. The 2-state fully recurrent NN with shared
weights and fewer parameters outperforms 1-state and 3-state ResNet with non-
shared weights.

####  5.2.2 Sharing Weights Across All Convolutional Layers (Less
Biologically-plausible)

Out of pure engineering interests, one could further push the limit of weight
sharing by not only sharing across time but also across states. Here we show
two 3-state ResNets that use a single set of convolutional weights across all
convolutional layers and achieve reasonable performance with very few
parameters (Figure 7).

Figure 7: A single set of convolutional weights is shared across all
convolutional layers in a 3-state ResNet. The transition (at time t= 4 and 8)
between nearby states is a 2x2 max pooling with stride 2. This means that each
state has a self-transition unrolled 3 times. “Feat.” denotes the number of
feature maps, which is the same across all 3 states. The learning rates were
the same as other experiments except that more epochs are used (i.e., 100, 20
and 20).

###  5.3 Experiment B: Multi-state Fully/Densely Recurrent Neural Networks

####  5.3.1 Shared vs. Non-shared Weights

Although an RNN is usually implemented with shared weights across time, it is
however possible to unshare the weights and use an independent set of weights
at every time $t$. For practical applications, whenever one can have a initial
$t=0$ and enumerate all possible $t$s, an RNN with non-shared weights should
be feasible, similar to the time-specific batch normalization described in
3.7. The results of 2-state fully recurrent neural networks with shared and
non-shared weights are shown in Figure 6.

####  5.3.2 The Effect of Readout Time

In visual cortex, useful information increases as time proceeds from the onset
of the visual stimuli. This suggests that recurrent system might have better
representational power as more time is allowed. We tried training and testing
the 2-state fully recurrent network with various readout time (i.e., unrolling
depth, see Section 3.5) and observe similar effects. See Figure 8.

Figure 8: A 2-state fully recurrent network with readout time t=2, 3, 5 or 10
(See Section 3.5 for the definition of readout time). There is consistent
performance improvement as $t$ increases. The number of parameters changes
since at some $t$, some recurrent connections have not been contributing to
the output and thus their number of parameters are subtracted from the total.

####  5.3.3 Larger Models With More States

We have shown the effectiveness of 2-state fully recurrent network above by
comparing it with 1-State ResNet. Now we discuss several observations
regarding 3-state and 4-state networks.

First, 3-state models seem to generally outperform 2-state ones. This is
expected since more parameters are introduced. With a 3-state models with
minimum engineering, we were able to get 7.47% validation error on CIFAR-10.

Next, for computational efficiency, we tried only allowing each state to have
transitions to adjacent states and to itself by disabling bypass connections
(e.g., V1-V3, V2-IT, etc.). In this case, the number of transitions scales
linearly as the number of states increases, instead of quadratically. This
setting performs well with 3-state networks and slightly less well with
4-state networks (perhaps as a result of small feature/parameter sizes). With
only adjacent connections, the models are no longer fully recurrent.

Finally, for 4-state fully recurrent networks, the models tend to become
overly computationally heavy if we train it with large $t$ or large number of
feature maps. With small $t$ and feature maps, we have not achieved better
performance than 3-state networks. Reducing the computational cost of training
multi-state densely recurrent networks would be an important future work.

For experiments in this subsection, we choose a moderately deep pre-net of
three 3x3 convolutional layers to model the layers between retina and V1:
Conv-BN-ReLU-Conv-BN-ReLU-Conv. This is not essential but outperforms shallow
pre-net slightly (within 1% validation error).

The results are shown in Figure 9.

Figure 9: The performance of 4-state and 3-state models. The state sizes of
the 4-state model are: 32x32x8, 16x16x16, 8x8x32, 4x4x64. The state sizes of
the 3-state model are: 32x32x64, 16x16x128, 8x8x256. 4-state models are small
since they are computationally heavy. The readout time is t=5 for both models.
All models are time-invariant systems (i.e., weights are shared across time).

####  5.3.4 Generalization Across Readout Time

As an RNN, our model supports training and testing with different readout
time. Based on our theoretical analyses in Section 2.2, the representation is
usually not guaranteed to converge when running a model with time
$t\to\infty$. Nevertheless, the model exhibits good generalization over time.
Results are shown in Figure 10. As a minor detail, the model in this
experiment has only adjacent connections and does not have any self-
transition, but we do not expect this to affect the conclusion.

Figure 10: Training and testing with different readout time. A 3-state
recurrent network with adjacent connections is trained with readout time t=10
and test with t=6,7,8,9,10,11,12 and 15.

##  6 Discussion

The dark secret of Deep Networks: trying to imitate Recurrent Shallow
Networks?

A radical conjecture would be: the effectiveness of most of the deep
feedforward neural networks, including but not limited to ResNet, can be
attributed to their ability to approximate recurrent computations that are
prevalent in most tasks with larger $t$ than shallow feedforward networks.
This may offer a new perspective on the theoretical pursuit of the long-
standing question “why is deep better than shallow” [22, 21].

Equivalence between Recurrent Networks and Turing Machines

Dynamical systems (in particular discrete time systems, that is difference
equations) are Turing universal (the game “Life" is a cellular automata that
has been demonstrated to be Turing universal). Thus dynamical systems such as
the feedback systems we discussed can be equivalent to Turing machine. This
offers the possibility of representing a computation more complex than a
single (for instance boolean) function with the same number of learnable
parameters. Consider for instance the powercase of learning a mapping $F$
between an input vector $x$ and an output vector $y=F(x)$ that belong to the
same n-dimensional space. The output can be thought as the asymptotic states
of the discrete dynamical system obtained iterating some map $f$. We expect
that in many cases the dynamical system that asymptotically performs the
mapping may have a much simpler structure than the direct mapping $F$. In
other words, we expect that the mapping $f$ such that $f^{(n)}(x)=F(x)$ for
appropriate, possibly very large $n$ can be much simpler than the mapping $F$
(here $f^{(n)}$ means the n-th iterate of the map $f$).

Empirical Finding: Recurrent Network or Residual Networks with weight sharing
work well

Our key finding is that multi-state time-invariant recurrent networks seem to
perform as well as very deep residual networks (each state corresponds to a
cortical area) without shared weights. On one hand this is surprising because
the number of parameters is much reduced. On the other hand a recurrent
network with fixed parameters can be equivalent to a Turing machine and
maximally powerful.

Conjecture about Cortex and Recurrent Computations in Cortical Areas

Most of the models of cortex that led to the Deep Convolutional architectures
and followed them – such as the Neocognitron [6], HMAX [26] and more recent
models [1] – have neglected the layering in each cortical area and the
feedforward and recurrent connections within each area and between them. They
also neglected the time evolution of selectivity and invariance in each of the
areas. The conjecture we propose in this paper changes this picture quite
drastically and makes several interesting predictions. Each area corresponds
to a recurrent network and thus to a system with a temporal dynamics even for
flashed inputs; with increasing time one expects asymptotically better
performance; masking with a mask an input image flashed briefly should disrupt
recurrent computations in each area; performance should increase with time
even without a mask for briefly flashed images. Finally, we remark that our
proposal, unlike relatively shallow feedforward models, implies that cortex,
and in fact its component areas are computationally as powerful a universal
Turing machines.

#### Acknowledgments

This work was supported by the Center for Brains, Minds and Machines (CBMM),
funded by NSF STC award CCF – 1231216.

## References

## References

  * [1] Using goal-driven deep learning models to understand sensory cortex.  Nature Neuroscience, 19,3:356–365, 2016.
  * [2] Christian Büchel and KJ Friston.  Modulation of connectivity in visual pathways by attention: cortical interactions evaluated with structural equation modelling and fmri.  Cerebral cortex, 7(8):768–778, 1997.
  * [3] Isaac Caswell, Chuanqi Shen, and Lisa Wang.  Loopy neural nets: Imitating feedback loops in the human brain.  CS231n Report, Stanford, http://cs231n.stanford.edu/reports2016/110_Report.pdf. Google Scholar time stamp: March 25th, 2016.
  * [4] David Eigen, Jason Rolfe, Rob Fergus, and Yann LeCun.  Understanding deep architectures using a recursive convolutional network.  arXiv preprint arXiv:1312.1847, 2013.
  * [5] Salah El Hihi and Yoshua Bengio.  Hierarchical recurrent neural networks for long-term dependencies.  Citeseer.
  * [6] Kunihiko Fukushima.  Neocognitron: A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position.  Biological Cybernetics, 36(4):193–202, April 1980.
  * [7] Alex Graves.  Generating sequences with recurrent neural networks.  arXiv preprint arXiv:1308.0850, 2013.
  * [8] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  arXiv preprint arXiv:1512.03385, 2015.
  * [9] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE International Conference on Computer Vision, pages 1026–1034, 2015.
  * [10] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Identity mappings in deep residual networks.  arXiv preprint arXiv:1603.05027, 2016.
  * [11] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [12] JM Hupe, AC James, BR Payne, SG Lomber, P Girard, and J Bullier.  Cortical feedback improves discrimination between figure and background by v1, v2 and v3 neurons.  Nature, 394(6695):784–787, 1998.
  * [13] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  arXiv preprint arXiv:1502.03167, 2015.
  * [14] Minami Ito and Charles D Gilbert.  Attention modulates contextual influences in the primary visual cortex of alert monkeys.  Neuron, 22(3):593–604, 1999.
  * [15] Alex Krizhevsky.  Learning multiple layers of features from tiny images, 2009.
  * [16] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in neural information processing systems, pages 1097–1105, 2012.
  * [17] Victor AF Lamme, Hans Super, and Henk Spekreijse.  Feedforward, horizontal, and feedback processing in the visual cortex.  Current opinion in neurobiology, 8(4):529–535, 1998.
  * [18] César Laurent, Gabriel Pereyra, Philémon Brakel, Ying Zhang, and Yoshua Bengio.  Batch normalized recurrent neural networks.  arXiv preprint arXiv:1510.01378, 2015.
  * [19] Ming Liang and Xiaolin Hu.  Recurrent convolutional neural network for object recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3367–3375, 2015.
  * [20] Qianli Liao, Joel Z Leibo, and Tomaso Poggio.  How important is weight symmetry in backpropagation?  arXiv preprint arXiv:1510.05067, 2015.
  * [21] Hrushikesh Mhaskar, Qianli Liao, and Tomaso Poggio.  Learning real and boolean functions: When is deep better than shallow.  arXiv preprint arXiv:1603.00988, 2016.
  * [22] Guido F Montufar, Razvan Pascanu, Kyunghyun Cho, and Yoshua Bengio.  On the number of linear regions of deep neural networks.  In Advances in neural information processing systems, pages 2924–2932, 2014.
  * [23] Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, and Yoshua Bengio.  How to construct deep recurrent neural networks.  arXiv preprint arXiv:1312.6026, 2013.
  * [24] Pedro HO Pinheiro and Ronan Collobert.  Recurrent convolutional neural networks for scene parsing.  arXiv preprint arXiv:1306.2795, 2013.
  * [25] Rajesh PN Rao and Dana H Ballard.  Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects.  Nature neuroscience, 2(1):79–87, 1999.
  * [26] M Riesenhuber and T Poggio.  Hierarchical models of object recognition in cortex.  Nature Neuroscience, 2(11):1019–1025, November 1999.
  * [27] Jürgen Schmidhuber.  Learning complex, extended sequences using the principle of history compression.  Neural Computation, 4(2):234–242, 1992.
  * [28] Thomas Serre, Aude Oliva, and Tomaso Poggio.  A feedforward architecture accounts for rapid categorization.  Proceedings of the National Academy of Sciences of the United States of America, 104(15):6424–6429, 2007.
  * [29] Richard Socher, Cliff C Lin, Chris Manning, and Andrew Y Ng.  Parsing natural scenes and natural language with recursive neural networks.  In Proceedings of the 28th international conference on machine learning (ICML-11), pages 129–136, 2011.
  * [30] Rupesh Kumar Srivastava, Klaus Greff, and Jürgen Schmidhuber.  Highway networks.  arXiv preprint arXiv:1505.00387, 2015.
  * [31] Simon Thorpe, Denis Fize, Catherine Marlot, et al.  Speed of processing in the human visual system.  nature, 381(6582):520–522, 1996.
  * [32] Andrea Vedaldi and Karel Lenc.  Matconvnet: Convolutional neural networks for matlab.  In Proceedings of the 23rd Annual ACM Conference on Multimedia Conference, pages 689–692. ACM, 2015.
  * [33] D.L.K. Yamins and J.D. Dicarlo.  Using goal-driven deep learning models to understand sensory cortex, 2016\.
  * [34] Matthew D Zeiler, Dilip Krishnan, Graham W Taylor, and Rob Fergus.  Deconvolutional networks.  In Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on, pages 2528–2535. IEEE, 2010.

##  Appendix A An Illustrative Comparison of a Plain RNN and a ResNet

Figure 11: A ResNet can be reformulated into a recurrent form that is almost
identical to a conventional RNN.

##  Appendix B Inhomogeneous, Time-invariant ResNet

The inhomogeneous, time-invariant version of ResNet is shown in Figure 12.

Figure 12: Inhomogeneous, Time-invariant ResNet

Let $K^{\prime}=K+I$, asymptotically we have:

$h=Ix+K^{\prime}h$ (7) $\displaystyle(I-K^{\prime})h=x$ (8) $\displaystyle
h=(I-K^{\prime})^{-1}x$ (9)

The power series expansion of above equation is:

$\displaystyle h=(I-K^{\prime})^{-1}x=(I+K^{\prime}+K^{\prime}\circ
K^{\prime}+K^{\prime}\circ K^{\prime}\circ K^{\prime}+...)x$ (10)

Inhomogeneous, time-invariant version of ResNet corresponds to the standard
ResNet with shared weights and shortcut connections from input to every layer.
If the model has only one state, it is experimentally observed that these
shortcuts undesirably add raw inputs to the final representations and degrade
the performance. However, if the model has multiple states (like the visual
cortex), it might be biologically-plausible for the first state (V1) to
receive constant inputs from the pre-net (retina and LGN). Figure 13 shows the
performance of an inhomogeneous 3-state recurrent network in comparison with
homogeneous ones.

Figure 13: Inhomogeneous 3-state models. All settings are the same as Figure
9. All models are time-invariant.


# Phase Transitions in the Classical Simulability of Open Quantum Systems

F. Azad  London Centre for Nanotechnology, University College London, Gordon
St., London, WC1H 0AH, United Kingdom  A. Hallam  School of Physics and
Astronomy, University of Leeds, Leeds LS2 9JT United Kingdom  J. Morley
London Centre for Nanotechnology, University College London, Gordon St.,
London, WC1H 0AH, United Kingdom  A. G. Green  London Centre for
Nanotechnology, University College London, Gordon St., London, WC1H 0AH,
United Kingdom

###### Abstract

We study the evolution of an open quantum system using a Langevin unravelling
of the density matrix evolution over matrix product states. As the strength of
coupling to and temperature of the environment is increased, we find a
transition where the entanglement of the individual trajectories saturates,
permitting a classical simulation of the system for all times. This is the
Hamiltonian open system counterpart of the saturation in entanglement found in
random circuits with projective or weak measurements. If a system is open,
there is a limit to the advantage in simulating its behaviour on a quantum
computer, even when that evolution harbours important quantum effects.

The term classical is applied to quantum systems in at least two different
ways. On the one hand, if a closed quantum system is in a weakly-entangled
state, it may be considered classical as long as the entanglement remains low.
In this limit, the equations of motion of the system are termed semi-
classical. On the other hand, an open system behaves classically once the
coupling to the environment has caused dephasing of the off-diagonal elements
of the density matrix. These two limits occur on very different timescales —
the semi-classical limit at early times and the dephasing limit at late times.

Can these views be reconciled and a classical description developed that works
from the earliest to the latest times? Recent insights have made steps towards
such an understanding for open systems. The transition in entanglement growth
in random circuits as a function of the rate of projective or weak measurement
allows a classical, weakly-entangled description of the system for all
timesnahum2017quantum ; li2018quantum ; nahum2018operator ; von2018operator ;
chan2019unitary ; li2019measurement ; szyniszewski2019entanglement .

Both are cousins to the quantum Zeno effect, by which frequent measurement in
a channel impedes transitions in that channelmisra1977zeno . The nature of the
many-body transition has been studied extensivelyszyniszewski2020universality
; li2020conformal ; chen2020emergent ; tang2020measurement ;
zhang2020nonuniversal ; zabalo2020critical ; gullans2020scalable ;
gullans2020dynamical ; jian2020measurement ; choi2020quantum ; bao2020theory
and similar analyses extended to measurement-induced transitions in open
Hamiltonian systemsgoto2020measurement ; alba2021spreading ;
fuji2020measurement ; PhysRevLett.126.170602 . These later cases are closely
related to a transition in classical describability as a function of coupling
to the environment, which we consider.

Our approach is to consider an unravelling of the density matrix evolution for
an open system over trajectories described using matrix product states. The
equations of motion of each trajectory can be considered a Langevin extension
of the time-dependent variational principle (TDVP) for matrix product states.
For closed systems, the TDVP equations constitute a semi-classical limit; they
correspond to classical Hamilton equations of motion on the variational
manifoldHaegemanTDVP ; hallam2019lyapunov . As the entanglement grows during
the Hamiltonian evolution, the TDVP equations break down as a larger and
larger variational manifold (measured by an exponentially growing bond order)
is required in order to capture the state and its dynamics. In this sense, the
semi-classical description is confined to early times. In our stochastic TDVP
Langevin equation, we find thresholds in the dynamics of individual
trajectories as a function of coupling to and temperature of the environment,
whereby the entanglement saturates at a low value, and a low-bond order
description gives high fidelity for all time. This quantum Zeno
phaseli2018quantum  constitutes a transition in the classical describability
of the open quantum system: the low bond order TDVP Langevin equation is an
effective semi-classical description that works for all time.

##  I A TDVP Langevin Equation

Figure 1:  Basic Properties of MPS Langevin Evolution: The TDVP Langevin
equation, Eq.(1) describes the evolution of the density matrix through an
ensemble of stochastic pure-state trajectories. Here, the von Neumann
entanglement as a function of time for a typical trajectory is shown
calculated is shown calculated for different coupling strengths. Here the bond
dimension is $D=128$ and we have kept temperature fixed, $T=0.2$, while we
vary the noise, $\gamma T$, i.e. coupling $\gamma$. Such plots show three
different regimes of behaviour. An initial transient, followed by an
approximately linear in $t$ or logarithmic in $t$ growth and finally a
saturation at the longest times. This saturation can either be determined by
the variational approximation or intrinsically by the interplay of the
Hamiltonian and the dissipative bath. It signals whether a lower bond-
dimension (hence less computationally intensive) simulation suffices in
underpinning the entropy dynamics.

Langevin equations describe the motion of a system coupled to an environment
(or alternatively the motion of slow collective degrees of freedom in an
effective bath described by the faster degrees of freedomzwanzig1960ensemble )
by adding noise and friction terms to the basic equations of motion of the
system. If the environmental degrees of freedom are in thermal equilibrium,
the friction and noise satisfy a fluctuation-dissipation relation. Applied to
quantum systems, the Schrödinger equation provides the basic equations of
motion. The ensemble of the resulting stochastic Schrödinger trajectories
recovers the density matrix evolution and is said to be an unravelling of it.

The Langevin equation over matrix product states studied here, can be written
in its Markovian limit as:

$\displaystyle\langle\partial_{i}\psi|\partial_{j}\psi\rangle\dot{X}_{j}$
$\displaystyle=$
$\displaystyle-i\langle\partial_{i}\psi|\hat{H}|\psi\rangle-i\sum_{n}\langle\partial_{i}\psi|\hat{F}_{n}|\psi\rangle\eta(t)$
(1)
$\displaystyle-i\sum_{n}\gamma\frac{\langle\psi|\hat{F}_{n}|\psi\rangle}{dt}\langle\partial_{i}\psi|\hat{F}_{n}|\psi\rangle.$

The terms on the left-hand side and the first term on the right constitute the
conventional time-dependent variational principle (TDVP) equationsHaegemanTDVP
. The second and third terms on the right are, respectively, the noise and
friction due to coupling to the environment. $\hat{F}_{n}$ are the operators
by which the system is coupled to the bath displacement operators. We
generally assume these to be spatially local. For spin-half chains they are
given by the $x,y$ and $z-$components of the spin operators on each site, each
of which couples to a separate bath. The bath is described as a collection of
harmonic oscillators and the noise-correlator is determined by the spectrum of
oscillators and the temperature of the bath;
$\langle\langle\eta(t)\eta(t^{\prime})\rangle\rangle=2\gamma
T\delta(t-t^{\prime})$ in the combined classical and Markovian limits.

The TDVP Langevin Equation is a new approach to unravelling the density matrix
evolution of an open system. It combines the study of quantum many-body
dynamics on a variational manifold — here we use matrix product states
(MPS)Orus:2014zl ; schollwock2011density ; perez2006matrix  — with the
Langevin limit of Keldysh field theorieskamenev2011field ; sieberer2016keldysh
. We combine these approaches by constructing a Keldysh path integral over MPS
statesgreen2016feynman . To date, MPS techniques have been employed in the
study of open systems largely in three ways: by starting with the Lindblad
master equation and either describing the density matrix directly as a matrix
product operatorverstraete2004matrix ; cui2015variational ;
weimer2019simulation  or else unravelling its evolution over MPS
representations of quantum trajectoriesdaley2009atomic ; daley2014quantum ;
bonnes2014superoperators ; alternatively, tensor network techniques can be
employed to describe the bath directly within the state
ansatzstrathearn2018efficient . Of these, the latter is numerically costly and
the Lindblad-based approaches formally apply in a different limitFergusThesis
. Questions of applicability aside, there are some features of our approach
that make it particularly attractive. It naturally leads to thermal
equilibrium, and it can treat the non-Markovian limit. There also exists a
natural hydrodynamic limit and Fokker-Plank description. The scaling of the
algorithm for each trajectory is the same as the usual TDVP for matrix product
states – moreover, there appears to be a degree of self averaging so that
certain quantities are well-approximated by relatively fewer trajectories.

A derivation of Eq.(1) is given in Appendix A. We develop the Langevin
equation from the Keldysh path integral for the time-evolution of the density
matrix. The method follows that of Ref.kamenev2011field  with the modification
that the Keldysh path integral is constructed over matrix product
statesgreen2016feynman . The result adds noise and friction to the time-
dependent variational principle constructed over matrix product
statesHaegemanTDVP . A similar construction for any variational class would
lead to a similar stochastic equation of motion, which we dub the TDVP-
Langevin equation. An alternative, heuristic derivation involves solving the
Schrödinger equation for the system and bath and substituting the solution of
the latter into the equation of motion of the former. Alternative unravellings
of the Lindblad equation for the density matrix evolution over
MPSdaley2009atomic  apply in different circumstances of relative time and
energy scales of the bath and system.

We construct Eq.(1) over matrix product states using conventional
methodsHaegemanTDVP . Integration of this equation is complicated by the
friction term. Naively, this requires inversion of a matrix that is
proportional to the system size and dimension of the variational manifold.
However, recognising that it consists of an outer product of vectors allows an
efficient inversion and integration of the equations of motion. Details are
given in Appendix B and our code is available at
https://github.com/AndrewHallam/Langevin.

Basic properties of the TDVP Langevin equation are summarised in Fig.1. A
thermal distribution over the variational manifold is given by a Boltzmann-
weighted Haar average over the variational manifold. In the case the of MPS of
bond dimension $D$, this average can be performed as a Haar integral over the
group $SU(dD)$ (with $d$ the local Hilbert space dimension). Such thermal
distributions are fixed points of the Langevin evolution (see Appendix A).
Fig. 1 captures the dynamics of the von Neumann entanglement entropy typical
for trajectories over the parameters we consider. In the case shown,
temperature $T=0.2$ is kept fixed and friction is increased. Beyond a critical
value, $\gamma T\approx 0.1$, entanglement growth is suppressed. The
entanglement undergoes a transition from being determined by the variational
approximation to becoming intrinsic to the interplay between the Hamiltonian
and the dissipative bath. Similar transitions occur when other parameters are
kept constant (friction $\gamma$, or noise $\gamma T$). We detail these
results in the following section.

##  II Open Evolution of a Rapidly Entangling System

In the absence of coupling to a bath, TDVP equations eventually fail as the
entanglement grows beyond that which can be represented on the variational
manifold111TDVP equations for the thermofield purification of the density
matrix may escape this fatehallam2019lyapunov  at least as far as local
observations are concerned. However, just as observed in projective
measurements of random circuits, the effects of the environment may restrict
the growth of entanglement. In extremis this might limit entanglement of
individual trajectories so that they can be represented on low dimensional
variational manifolds. The TDVP Langevin equation will then give a good
account of the dynamics at all times, signifying a transition in its classical
representability. This is our interpretation of the sequence of results
presented in this section.

The Hamiltonian that we consider is the tilted field Ising model

$\hat{H}=-\sum_{i}\left[J\sigma_{i}^{z}\sigma_{i+1}^{z}+h\sigma^{z}_{i}+g\sigma^{x}_{i}\right],$
(2)

with $J=1$, $g=-1.05$ and $h=0.5$. With these parameters, the Hamiltonian is a
far from any integrable point and rapidly thermalisingbanuls2011strong ;
leviatan2017quantum .

Infinite temperature and vanishing friction: Fig.2 shows the variation in von
Neumann entanglement across the central bond as a function of time for
simulations with a range of bond orders and noise strengths. The broad result
of these simulations is that the entanglement saturates at long times at a
value determined by the bond order of the simulation. This is consistent with
an infinite temperature final state with the maximum entanglement supported by
the variational manifold. The most interesting aspect of these results is the
decreasing rate of early-time entanglement growth with increasing noise
strength. Crucially we do not find evidence of an intrinsic saturation of
entanglement – only that dictated by the limitations of the variational
approximation.

Figure 2:  Evolution of Entanglement at Infinite Temperature: Here we consider
the evolution under Eq.(1) with Hamiltonian Eq.(2) with $J=1$, $g=-1.05$ and
$h=0.5$. The friction coefficient $\gamma=0$ with noise $\gamma T$ finite and
initial state $Z$. a) von Neumann Entropy versus time at a fixed bond order
$D=160$ for different values of $\gamma T$. Without noise, and after an
initial transient, the entanglement grows linearly with time before saturating
at a constant value less than the maximum determined by the MPS manifold. With
non-zero noise, the entanglement growth is linear in time with a growth rate
that reduces with increasing $\gamma T$ due to the dephasing effects of the
bath. Ultimately, all curves saturate (not shown) to the same value
characteristic of an infinite temperature state on the variational manifold.
b) The growth rate of entanglement $\dot{S}$, extracted beyond the initial
transient ($t\approx 4$), versus noise.

Finite temperature and friction: Including both noise and friction, we do see
such an intrinsic saturation. This is demonstrated in two ways; by considering
the saturation of entanglement at long times and by a high fidelity between
low- and high-bond order simulations at long times.
i. Saturating entanglement In order to demonstrate this, we first show in Fig.
3 the long-time average of the von-Neumann entropy. A graph showing the
typical time-dependence from which such saturation values are computed is
shown in Fig.1b). For low noise and friction, the saturation is determined by
the limitations of the variational manifold. Panels a), b) and c) show that as
a function of $\gamma T$ at fixed $\gamma$, then $T$, and $\gamma T$ at fixed
$\gamma$, respectively. A threshold is reached for each bond order where it
adequately captures the saturation entanglement, thus indicating a transition
to increasingly classically simulatable dynamics. The transition can be seen
from the point where the trajectories obtained at different bond orders give
the same saturation entanglement. From this we can extract a critical $\gamma$
or $\gamma T$ as a function of bond order that we show in each corresponding
inset figure.
ii. High fidelity as  $t\rightarrow\infty$: We can identify an analogous
transition in the fidelity of each trajectory at different bond orders versus
a reference trajectory with bond order $D=128$. In this case, we find that
beyond a critical combination of $\gamma$ or $\gamma T$, the fidelity of the
state at low bond dimension remains close to $1$ for long times. We expand
upon this result in Fig. 4, where we identify a divergent classical simulation
time. We note that the fidelity is more sensitive to the time-step as friction
is increased – an issue typical of numerical integration of systems of
stochastic differential equations. This makes accessing the critical point of
the transition numerically intensive for the parameters and Hamiltonian we
consider. The entropy is less sensitive to this.

Figure 3:  Evolution of Saturation Entanglement at Finite Temperature: Here we
consider the evolution under Eq.(1) with Hamiltonian Eq.(2) where $J=1$,
$g=-1.05$ and $h=0.5$ and finite $\gamma$ and $T$. In the main figures we show
the dependence of the von Neumann entropy as a function of noise and friction:
a) versus $\gamma T$ at fixed $\gamma$, b) versus $\gamma T$ at fixed $T$, c)
versus $\gamma$ at fixed $\gamma T$. In each case, at low values of noise and
friction, the saturation entanglement $\bar{S}$, is determined by the choice
of variational parametrization through the bond order. As the noise and
friction are increased, there is a cross-over where the saturation
entanglement decreases from this value. Each bond order captures the
saturation entanglement for a sufficiently large noise and friction. This is
indicated when the entanglement begins to follow the entanglement given by the
highest bond order simulation. Where the saturation entanglement of the
individual trajectories (distinguished by bond order) dramatically fall to a
single curve, so that the lower bond dimension curves coincide with higher
bond dimensions, illustrates this transition to situations where the
saturation entanglement is intrinsically determined. In each corresponding
inset figure, we have extracted critical dissipation strengths where these
transitions occur as a function of bond order.

Figure 4:  Divergent Classical Simulation Time: Simulations carried out at a
bond order $D$ give a good account of the system evolution up to a time
$t^{*}(D)$. We extract these values versus a reference $D=128$ simulation,
which serves as the good account of the system. We do this in two ways, by
comparing the difference in von Neumann entanglement entropy between these
states and the fidelity with this state. $t^{*}(D)$ is the time when the
simulation with varying bond dimension deviates appreciably from the reference
trajectory. The row with panels a), b), c), shows $t^{*}(D)$ extracted from
the entropy, while panels d), e), f) demonstrate this for the fidelity. The
fixed variables are split across the columns – a), d), shows varying $\gamma
T$ at fixed $\gamma=0.2$, b), e) $\gamma T$ at fixed $T=0.2$, and c), f)
$\gamma$ at fixed $\gamma T=0.25$. Panels g) and h) show typical evolution of
entanglement and fidelity with time. This is data from the near-critical point
for $T=0.2$, which is for the trajectory with noise $\gamma T=0.1$. We say
that a simulation has failed to provide a good account of the system when the
trajectory deviates beyond $\epsilon=0.05$, and the time at which this occurs
is $t^{*}(D)$. In g), this is the point where $\Delta
S/S_{D=128}=|S_{D=128}-S_{D}|/S_{D=128}>\epsilon$. Analogously in h),
$t^{*}(D)$ is the time when the fidelity is appreciably different to $1$, i.e.
$|\langle\psi_{D=128}(t)|\psi_{D}(t)\rangle|<1-\epsilon$. A divergent
$t^{*}(D)$, within either method of extraction, indicates a transition in the
classical simulability of the open quantum system.

##  III Discussion

This work introduces a new method to investigate the dynamics of open many-
body quantum systems, the TDVP-Langevin equation. We derive this by
considering an appropriate limit of the Keldysh path integral constructed over
the MPS manifold. Our investigations reveal a phase transition in the
applicability of this approach as a function of coupling to the environment –
when the bath temperature and induced friction are sufficiently high,
entanglement growth in individual trajectories is suppressed, and a low bond
order description works for all time. This is a transition in the classical
simulabiity of the open quantum system.

We believe that this transition is related to several other transitions in
quantum dynamics that have been observed as a function of coupling the the
environment or measurement, including the restriction of entanglement growth
in random circuits with projective measurement, the quantum Zeno effect (and
perhaps the KT transition in the spin-boson model leggett1987dynamics ;
Florens2010quantum ; barratt2020dissipative ).

The implications of this result may be far-reaching. In the context of using
the TDVP-Langevin equation to simulate open quantum systems, an efficient
description for long times is possible for systems in the many-body quantum
Zeno phase. Indeed, when a target system is in such a phase, there is no
(asymptotic) advantage in using a quantum computer to simulate it. Since many
chemical reactions of potential interest for quantum computation occur
embedded in a dissipative aqueous environment, this is certainly a point
worthy of consideration.

Moreover, viewed from the perspective of a description of the quantum
computational device, the transition into the many-body Zeno phase might
indicate transitions in the ability to solve quantum problems. While
thresholds of noise for quantum error correction have been identified in the
case of gate-based quantum computation, no such thresholds currently exist for
adiabatic computation. It is intriguing to speculate that determining whether
a putative adiabatic computational device is in its Zeno phase or not might
provide similar bounds on performancecrowley2014quantum ;
barratt2020dissipative .

We envisage a number of ways in which this work might be developed. Extending
the approach to local observables in closed quantum systems presents some
exciting possibilities. In this case the bath would refer to other elements of
the system itself and its properties self-consistently determined through the
evolutionzwanzig1960ensemble ; zwanzig1961memory , Such a description has the
promise of connecting early-time semi-classical descriptions to late-time
hydrodynamics and thermalisation. Exploring the Fokker-Planck formulation of
the TDVP-Langevin equation would bring a complementary perspective to our
analysisberta2017thermal .

The accurate description of a quantum system from early to late times is
generally not possible because of growing entanglement. However, coupling to
the environment can limit this growth and render this achievable. This work
has coordinated physical insights from several different perspectives to
develop such a numerical scheme. We hope both that the algorithm itself will
prove useful and that it will inspire further insights.

##  IV Acknowledgements

We gratefully acknowledge funding from the EPSRC under grants EP/L015242/,
EP/S005021/1 and EP/R020612/1.

## References

  * (1) Nahum, A., Ruhman, J., Vijay, S. & Haah, J.  Quantum entanglement growth under random unitary dynamics.  _Physical Review X_ 7, 031016 (2017).
  * (2) Li, Y., Chen, X. & Fisher, M. P.  Quantum zeno effect and the many-body entanglement transition.  _Physical Review B_ 98, 205136 (2018).
  * (3) Nahum, A., Vijay, S. & Haah, J.  Operator spreading in random unitary circuits.  _Physical Review X_ 8, 021014 (2018).
  * (4) von Keyserlingk, C., Rakovszky, T., Pollmann, F. & Sondhi, S.  Operator hydrodynamics, otocs, and entanglement growth in systems without conservation laws.  _Physical Review X_ 8, 021013 (2018).
  * (5) Chan, A., Nandkishore, R. M., Pretko, M. & Smith, G.  Unitary-projective entanglement dynamics.  _Physical Review B_ 99, 224307 (2019).
  * (6) Li, Y., Chen, X. & Fisher, M. P.  Measurement-driven entanglement transition in hybrid quantum circuits.  _Physical Review B_ 100, 134306 (2019).
  * (7) Szyniszewski, M., Romito, A. & Schomerus, H.  Entanglement transition from variable-strength weak measurements.  _Physical Review B_ 100, 064204 (2019).
  * (8) Misra, B. & Sudarshan, E. G.  The zeno’s paradox in quantum theory.  _Journal of Mathematical Physics_ 18, 756–763 (1977).
  * (9) Szyniszewski, M., Romito, A. & Schomerus, H.  Universality of entanglement transitions from stroboscopic to continuous measurements.  _Physical review letters_ 125, 210602 (2020).
  * (10) Li, Y., Chen, X., Ludwig, A. W. & Fisher, M.  Conformal invariance and quantum non-locality in hybrid quantum circuits.  _arXiv preprint arXiv:2003.12721_ (2020).
  * (11) Chen, X., Li, Y., Fisher, M. P. & Lucas, A.  Emergent conformal symmetry in nonunitary random dynamics of free fermions.  _Physical Review Research_ 2, 033017 (2020).
  * (12) Tang, Q. & Zhu, W.  Measurement-induced phase transition: A case study in the nonintegrable model by density-matrix renormalization group calculations.  _Physical Review Research_ 2, 013022 (2020).
  * (13) Zhang, L. _et al._ Nonuniversal entanglement level statistics in projection-driven quantum circuits.  _Physical Review B_ 101, 235104 (2020).
  * (14) Zabalo, A. _et al._ Critical properties of the measurement-induced transition in random quantum circuits.  _Physical Review B_ 101, 060301 (2020).
  * (15) Gullans, M. J. & Huse, D. A.  Scalable probes of measurement-induced criticality.  _Physical review letters_ 125, 070606 (2020).
  * (16) Gullans, M. J. & Huse, D. A.  Dynamical purification phase transition induced by quantum measurements.  _Physical Review X_ 10, 041020 (2020).
  * (17) Jian, C.-M., You, Y.-Z., Vasseur, R. & Ludwig, A. W.  Measurement-induced criticality in random quantum circuits.  _Physical Review B_ 101, 104302 (2020).
  * (18) Choi, S., Bao, Y., Qi, X.-L. & Altman, E.  Quantum error correction in scrambling dynamics and measurement-induced phase transition.  _Physical Review Letters_ 125, 030505 (2020).
  * (19) Bao, Y., Choi, S. & Altman, E.  Theory of the phase transition in random unitary circuits with measurements.  _Physical Review B_ 101, 104301 (2020).
  * (20) Goto, S. & Danshita, I.  Measurement-induced transitions of the entanglement scaling law in ultracold gases with controllable dissipation.  _Physical Review A_ 102, 033316 (2020).
  * (21) Alba, V. & Carollo, F.  Spreading of correlations in markovian open quantum systems.  _Physical Review B_ 103, L020302 (2021).
  * (22) Fuji, Y. & Ashida, Y.  Measurement-induced quantum criticality under continuous monitoring.  _Physical Review B_ 102, 054302 (2020).
  * (23) Alberton, O., Buchhold, M. & Diehl, S.  Entanglement transition in a monitored free-fermion chain: From extended criticality to area law.  _Phys. Rev. Lett._ 126, 170602 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.126.170602.
  * (24) Haegeman, J. _et al._ Time-dependent variational principle for quantum lattices.  _Physical Review Letters_ 107, 070601 (2011).
  * (25) Hallam, A., Morley, J. & Green, A. G.  The lyapunov spectra of quantum thermalisation.  _Nature communications_ 10, 1–8 (2019).
  * (26) Zwanzig, R.  Ensemble method in the theory of irreversibility.  _The Journal of Chemical Physics_ 33, 1338–1341 (1960).
  * (27) Orús, R.  A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  _Annals of physics_ 349, 117–158 (2014).
  * (28) Schollwöck, U.  The density-matrix renormalization group in the age of matrix product states.  _Annals of Physics_ 326, 96–192 (2011).
  * (29) Perez-Garcia, D., Verstraete, F., Wolf, M. M. & Cirac, J. I.  Matrix product state representations.  _arXiv preprint quant-ph/0608197_ (2006).
  * (30) Kamenev, A.  _Field theory of non-equilibrium systems_ (Cambridge University Press, 2011).
  * (31) Sieberer, L. M., Buchhold, M. & Diehl, S.  Keldysh field theory for driven open quantum systems.  _Reports on Progress in Physics_ 79, 096001 (2016).
  * (32) Green, A., Hooley, C., Keeling, J. & Simon, S.  Feynman path integrals over entangled states.  _arXiv preprint arXiv:1607.01778_ (2016).
  * (33) Verstraete, F., Garcia-Ripoll, J. J. & Cirac, J. I.  Matrix product density operators: Simulation of finite-temperature and dissipative systems.  _Physical review letters_ 93, 207204 (2004).
  * (34) Cui, J., Cirac, J. I. & Banuls, M. C.  Variational matrix product operators for the steady state of dissipative quantum systems.  _Physical review letters_ 114, 220601 (2015).
  * (35) Weimer, H., Kshetrimayum, A. & Orús, R.  Simulation methods for open quantum many-body systems.  _arXiv preprint arXiv:1907.07079_ (2019).
  * (36) Daley, A., Taylor, J., Diehl, S., Baranov, M. & Zoller, P.  Atomic three-body loss as a dynamical three-body interaction.  _Physical review letters_ 102, 040402 (2009).
  * (37) Daley, A. J.  Quantum trajectories and open many-body quantum systems.  _Advances in Physics_ 63, 77–149 (2014).
  * (38) Bonnes, L. & Läuchli, A. M.  Superoperators vs. trajectories for matrix product state simulations of open quantum system: a case study.  _arXiv preprint arXiv:1411.4831_ (2014).
  * (39) Strathearn, A., Kirton, P., Kilda, D., Keeling, J. & Lovett, B. W.  Efficient non-markovian quantum dynamics using time-evolving matrix product operators.  _Nature communications_ 9, 1–9 (2018).
  * (40) Barratt, F.  _The dynamics and control of quantum information out of equilibrium_.  Ph.D. thesis, Kings College London (2021).
  * (41) TDVP equations for the thermofield purification of the density matrix may escape this fatehallam2019lyapunov  at least as far as local observations are concerned.
  * (42) Bañuls, M. C., Cirac, J. I. & Hastings, M. B.  Strong and weak thermalization of infinite nonintegrable quantum systems.  _Physical review letters_ 106, 050405 (2011).
  * (43) Leviatan, E., Pollmann, F., Bardarson, J. H., Huse, D. A. & Altman, E.  Quantum thermalization dynamics with matrix-product states.  _arXiv preprint arXiv:1702.08894_ (2017).
  * (44) Leggett, A. J. _et al._ Dynamics of the dissipative two-state system.  _Reviews of Modern Physics_ 59, 1 (1987).
  * (45) Florens, S., Venturelli, D. & Narayanan, R.  Quantum phase transition in the spin boson model.  In _Quantum Quenching, Annealing and Computation_ , 145–162 (Springer, 2010).
  * (46) Barratt, F. _et al._ Dissipative failure of adiabatic quantum transport as a dynamical phase transition.  _arXiv preprint arXiv:2012.15212_ (2020).
  * (47) Crowley, P., Đurić, T., Vinci, W., Warburton, P. & Green, A.  Quantum and classical dynamics in adiabatic computation.  _Physical Review A_ 90, 042317 (2014).
  * (48) Zwanzig, R.  Memory effects in irreversible thermodynamics.  _Physical Review_ 124, 983 (1961).
  * (49) Berta, M., Brandao, F. G., Haegeman, J., Scholz, V. B. & Verstraete, F.  Thermal states as convex combinations of matrix product states.  _arXiv preprint arXiv:1709.07423_ (2017).
  * (50) Frenkel, J.  _Wave mechanics, advanced general theory_ , vol. 1 (Oxford, 1934).
  * (51) Kamenev, A.  Keldysh and doi-peliti techniques for out-of-equilibrium systems.  In _Strongly Correlated Fermions and Bosons in Low-Dimensional Disordered Systems_ , 313–340 (Springer, 2002).
  * (52) Crowley, P. J. & Green, A.  Anisotropic landau-lifshitz-gilbert models of dissipation in qubits.  _Physical Review A_ 94, 062106 (2016).
  * (53) Haegeman, J., Lubich, C., Oseledets, I., Vandereycken, B. & Verstraete, F.  Unifying time evolution and optimization with matrix product states.  _Physical Review B_ 94, 165116 (2016).

##  Appendix A Haar averaging on the MPS manifold

In the left canonical form, MPS of left bond order $D$ and local Hilbert space
dimension $D$ are given by $SU(dD)$ isometriesOrus:2014zl . In the case of a
finite chain of length $L$, the left bond order at the $n^{th}$ site
$D_{n}=d^{n}$ up to the maximum bond order at site $n=\log_{d}D_{\hbox{max}}$.
The bond order remains $D_{\hbox{max}}$ up to site
$L-1+\log_{d}D_{\hbox{max}}$ beyond which it reduces as $D_{n}=d^{(n-L+1)}$.

The thermal expectation of an operator can be calculated by a Boltzmann-
weighted Haar average over this variational manifold. The average of an
operator $\hat{O}$ is given by

$\displaystyle\langle\langle\hat{O}\rangle\rangle$ $\displaystyle=$
$\displaystyle\frac{\int\prod_{n}DA_{n}\langle\phi|\hat{O}|\phi\rangle\exp\left[-\beta\langle\phi|\hat{H}|\phi\rangle\right]}{\int\prod_{n}DA_{n}\exp\left[-\beta\langle\phi|\hat{H}|\phi\rangle\right]}.$
(3)

The expectation of the Hamiltonian $\langle\phi|\hat{H}|\phi\rangle$ and the
operator $\langle\phi|\hat{O}|\phi\rangle$ are calculated by usual MPS
techniques. In practice, we carry out the integrals by by sampling a Haar
random distribution of isometric MPS tensors; $A^{\sigma}_{ij}\equiv
U_{i\otimes\sigma,j\otimes 1}\in SU(dD_{n})$ . These are obtained by a QR
decomposition of a tensor with elements drawn randomly from a normal
distribution.

Fig. 5 shows this Boltzmann-weighted Haar distribution as a function of energy
at different temperatures. These thermal distributions are fixed-points of the
dynamics described by Eq.(1)

Figure 5:  Density of occupied states at finite-temperature for a length 15
chain at bond order $D_{\hbox{max}}=2$: a) The density of occupied states is
computed for the Hamiltonian Eq.(2) with $J=1$ $g=-1.05$, and $h=0.5$ from a
sample of Haar random distributed isometric MPS initial states. At $T=\infty$
— an unweighted Haar average — the majority of the states are in the middle of
the spectrum and are highly entangled. A Boltzmann weight shifts this
distribution to lower energy. Low-entanglement states are a subset of measure
zero in the thermodynamic limit. b) The Boltzmann-weighted Haar average is a
fixed point of the TDVP-Langevin equation. Here we compare distributions
obtained after evolving with the Hamiltonian for $t=1$, $2$, and $3$, with
$\gamma=0.1$ and $T=0.1$. The distribution for $T=0.1$ is shown in bold in a),
we see the Langevin evolution preserves this distribution.

##  Appendix B Deriving the TDVP Langevin Equation

Here we outline how the MPS TDVP Langevin equation can be obtained. We follow
two separate routes: a heuristic route using a solution of the Schródinger
equation for the bath and system, and constructing a Langevin limit of the
Keldysh path integral constructed over MPS.

Bath Model: We model the bath as a collection of independent non-interacting
harmonic oscillators. These are coupled to local system operators
$\hat{F}_{n}$ at site $n$ by their displacements. The Hamiltonian for the bath
and its coupling to the system are given by

$\displaystyle\hat{H}_{bath}$ $\displaystyle=$
$\displaystyle\sum_{n}\sum_{\alpha}\hbar\omega_{\alpha}\left(\hat{a}^{\dagger}_{n,\alpha}\hat{a}_{n,\alpha}+1/2\right),$
$\displaystyle\hat{H}_{I}$ $\displaystyle=$
$\displaystyle-\sum_{n}\sum_{\alpha}\lambda_{\alpha}\left(\hat{a}^{\dagger}_{n,\alpha}+\hat{a}_{n,\alpha}\right)\hat{F}_{n},$
(4)

respectively. The index $\alpha$ labels the different oscillator modes at the
site $n$. The distribution and temperature of the oscillator modes is assumed
to be the same at each site. Moreover, we assume no back reaction of the
system on the bath, so that the bath distribution remains equilibrium. This is
a subtle assumption - the bath must be non-linear in order to thermalise
energy absorbed from the system, but these non-linearities must operate on
timescales such that the bath’s effect on the system is the same as
independent oscillators. The assumptions are standard - but non-trivial - and
permit the simple manipulations that follow.

A heuristic derivation of the TDVP Langevin equation can be made in the spirit
of the Frenkel principle for deriving the TDVP equations frenkel1934wave . The
state of the system and bath is parametrized as
$|\psi(\bf{z})\rangle\otimes_{n,\alpha}|\phi_{n,\alpha}\rangle$ where ${\bf
z}$ corresponds to some set of variational parameters of the system and
$\phi_{n,\alpha}$ are coherent state parameters of the $\alpha$-oscillator on
site $n$. The time derivative of the wavefunction in the Schrödinger equation
for the system and bath is expanded in a chain rule over the system and bath
parameters:

$\displaystyle\dot{\bf
z}|\partial_{z_{i}}\psi\rangle\otimes_{n,\alpha}|\phi_{n,\alpha}\rangle+|\psi\rangle\sum_{m,\beta}\dot{\phi}_{m,\beta}\partial_{\phi_{m,\beta}}\left(\otimes_{n,\alpha}|\phi_{n,\alpha}\rangle\right)$
$\displaystyle\approx$
$\displaystyle-i\hat{H}|\psi\rangle\otimes_{n,\alpha}|\phi_{\alpha}\rangle,$

where the inequality is because the dynamics might take the state outside of
the variational manifold. Taking an inner product with
$\langle\partial_{z_{i}}\psi|\otimes_{n,\alpha}\langle\phi_{n,\alpha}|$ allows
us to obtain an equation of motion for the system in the presence of the bath,
and an inner product with
$\langle\psi|\;\partial_{\phi_{o,\gamma}}\left(\otimes_{n,\alpha}\langle\phi_{n,\alpha}|\right)$
allows us to obtain an equation of motion for the bath in the presence of the
system:

$\displaystyle
i\langle\partial_{z_{i}}\psi|\partial_{z_{j}}\psi\rangle\dot{z}_{j}=\langle\partial_{z_{i}}\psi|\hat{H}_{S}|\psi\rangle$
$\displaystyle+\sum_{n}\langle\partial_{z_{i}}\psi|\hat{F}_{n}|\psi\rangle\sum_{\alpha}\tilde{\lambda}_{\alpha}(\bar{\phi}_{n,\alpha}(t)+\phi_{n,\alpha}(t))$
$\displaystyle
i\dot{\phi}_{n,\alpha}=\omega_{\alpha}+\lambda_{\alpha}\langle\psi|\hat{F}_{n}|\psi\rangle$
(5)

The equalities are attained since the inner products with the tangent vectors
project the Hamiltonian evolution back onto the variational manifold. These
equations correspond to the usual TDVP equations with addition of coupling
between the system and bath.

The remaining steps involve integrating the equation of motion for the bath
degrees of freedom and substituting back into the equation of motion for the
system. The formal solution of the equation of motion of the bath degrees of
freedom is

$\phi_{n,\alpha}(t)=\phi_{n,\alpha}(0)e^{-i\omega_{\alpha}t}-\int_{0}^{\infty}D^{R}(t^{\prime}-t)\langle\psi|\hat{F}_{n}|\psi\rangle
dt^{\prime}$

where we have identified the retarded bath correlator
$D^{R}(t^{\prime}-t)=-i\Theta(t-t^{\prime})\langle\hat{a}^{\dagger}(t^{\prime})\hat{a}(t)\rangle=ie^{i\omega(t^{\prime}-t)}$.
The assumptions of no back reaction of the system on the bath are taken
account of by treating the terms
$\eta_{n}(t)=\sum_{\alpha}\phi_{n,\alpha}(0)e^{-i\omega_{\alpha}t}$ as a
stochastic random field with variance appropriate to the thermal distribution.
After further identifying $\partial_{t}\Gamma(t)=D^{R}(t)$ the system of
equation of motion can be put in the form

$\displaystyle
i\langle\partial_{A_{i}}\psi|\partial_{A_{j}}\psi\rangle\dot{A}_{j}=\langle\partial_{A_{i}}\psi|\hat{H}_{S}|\psi\rangle+\sum_{n}\langle\partial_{A_{i}}\psi|\hat{F}_{n}|\psi\rangle\left(\int
dt^{\prime}\Gamma(t-t^{\prime})\langle\partial_{A_{j}}\psi|\hat{F}_{n}|\psi\rangle\dot{A}_{j}+\eta(t)\right)$
(6)

the Markovian limit of which recovers Eq.(1).

A formal derivation from the Keldysh path integral can also be made. We follow
the approach described in Ref.kamenev2002keldysh ; kamenev2011field ;
crowley2016anisotropic  constructing a Keldysh path integral for the density
matrix and integrate out the bath in an appropriate limit to obtain a Langevin
equation. Our main modification is to construct the path integral over MPS
following Ref.green2016feynman .

The initial density matrix is assumed to factorise into density matrices for
the system and bath as $\hat{\rho}\otimes\hat{\rho}_{bath}$ \- the bath being
in thermal equilibrium and its distribution assumed to be unchanged in time.
This evolves to
$\bar{T}e^{-i\int_{0}^{t}dt^{\prime}(\hat{H}+\hat{H}_{int}+\hat{H}_{bath})}\hat{\rho}\otimes\hat{\rho}_{bath}Te^{-i\int_{0}^{t}dt^{\prime}(\hat{H}+\hat{H}_{int}+\hat{H}_{bath})}$
at time $t$. Construction of the Keldysh path integral proceeds by dividing up
the time-ordered ($T$) and anti-time-ordered ($\bar{T}$) exponentials into
many infinitessimal evolutions and inserting resolutions of the identity over
coherent state variables for the bath and using a Haar measure over MPS states
for the systemgreen2016feynman :

$\displaystyle{\bm{1}}$ $\displaystyle=$
$\displaystyle\int\prod_{n}d\bar{\phi}_{n}(t)d\phi_{n}(t)e^{-\sum_{n}\bar{\phi}_{n}\phi_{n}}\otimes_{n}|\bar{\phi}_{n}\rangle\langle\phi_{n}|$
$\displaystyle{\bm{1}}$ $\displaystyle=$ $\displaystyle\int
DA(t)|A\rangle\langle A|$

Here $DA$ symbolises an integral over all tensors in the MPS chain using the
Haar measure introduced in Ref.green2016feynman . A Keldysh rotation
transforms from the fields on the time-ordered contour ($+$) and the anti
time-ordered contour ($-$) to the sum and difference between them
$A^{\pm}=A^{cl}\pm A^{q}$, $\phi^{\pm}=\phi^{cl}\pm\phi^{q}$, known as the
classical and quantum components of the quantum field. These manipulations
give the following path integral for the evolution kernal of the density
matrix:

$\displaystyle K(t)$ $\displaystyle=$ $\displaystyle\int
DAD(\bar{\phi},\phi)e^{S[A^{c},A^{q},\phi^{c},\phi^{q}]}$ $\displaystyle S$
$\displaystyle=$ $\displaystyle S[A^{cl}+A^{q}]-S[A^{cl}-A^{q}]+\int
dt\Bigg{[}\sum_{n,\alpha}\lambda_{\alpha}\left(F^{q}_{n},F^{cl}_{n}\right)\left(\begin{array}[]{c}\phi^{cl}_{n,\alpha}+\bar{\phi}^{cl}_{n,\alpha}\\\
\phi^{q}_{n,\alpha}+\bar{\phi}^{q}_{n,\alpha}\end{array}\right)\Bigg{]}$
$\displaystyle+\int
dt\;dt^{\prime}\Bigg{[}\sum_{n,\alpha}\left(\bar{\phi}^{cl}_{n,\alpha}(t),\bar{\phi}^{q}_{n,\alpha}(t)\right)\begin{pmatrix}0&[D_{\alpha}^{A}]^{-1}(t-t^{\prime})\\\
[D_{\alpha}^{R}]^{-1}(t-t^{\prime})&[D_{\alpha}^{-1}]^{K}(t-t^{\prime})\\\
\end{pmatrix}(t-t^{\prime})\left(\begin{array}[]{c}\phi^{cl}_{n,\alpha}(t^{\prime})\\\
\phi^{q}_{n,\alpha}(t^{\prime})\end{array}\right)\Bigg{]}$

where $D_{\alpha}^{A}$, $D_{\alpha}^{R}$ and $D_{\alpha}^{K}$ are the
advanced, retarded and Keldysh components of the bath Greens function.
Consistent with our assumption of a thermal equilibrium bath and no back
reaction, they are related by the fluctuation dissipation relation:
$D^{K}_{\alpha}(\omega)=\coth(\omega/2T)[D^{R}_{\alpha}(\omega)-D^{A}_{\alpha}(\omega)]$
with $D^{R(A)}=1/(\omega\pm i\delta)$. $S[A]$ is the action of the system in
the absence of coupling to the bath.

The simple quadratic form of the bath action follows from our assumptions and
modeling of its affects as independent harmonic oscillators. It enables one to
integrate out the bath and, depending upon timescalesFergus, to construct
either a Lindblad or Langevin limit. We construct the latter limit in three
steps. First the bath degrees of freedom are integrated out. The resulting
dissipative contribution to the action has cross terms between classical and
quantum components of the expectations of $\hat{F}$, and a term quadratic in
the quantum component.,

$S_{diss}=\int
dt\;dt^{\prime}\sum_{n}\left(F^{cl}_{n},F^{q}_{n}\right)(t)\begin{pmatrix}0&D^{A}(t-t^{\prime})\\\
D^{R}(t-t^{\prime})&D^{K}(t-t^{\prime})\\\
\end{pmatrix}\left(\begin{array}[]{c}F^{cl}_{n}\\\
F^{q}_{n,}\end{array}\right)(t^{\prime}),$

where the bath propagators without indices indicate a sum over all modes, for
example $D^{R}=\sum_{\alpha}\lambda^{2}_{\alpha}D^{R}_{\alpha}$. The quadratic
term in $F^{q}$ is decoupled with a Hubbard-Stratonovich field $\eta(t)$ that
ultimately will play the role of the stochastic noise field in Eq.(1). The
final trick to bring this integral to the Langevin form is to Taylor expand
the action to linear order in the in the quantum fields, $A^{q}$. The result
is a path integral over the MPS tensors $A^{q}$ and $A^{cl}$ and the noise
field $\eta$:

$\displaystyle K(t)$ $\displaystyle=$ $\displaystyle\int DA^{q}DA^{cl}D\eta
e^{iS[A^{c},A^{q},\eta]}$ $\displaystyle S$ $\displaystyle=$
$\displaystyle\int dt\sum_{n}A^{q}_{n}(t)\underbrace{\Bigg{[}2\delta
S[A^{cl}]/\delta A_{n}^{cl}(t)+2\sum_{m}\partial F^{cl}_{m}/\partial
A_{n}^{cl}(t)\left(\int
dt^{\prime}\;D^{R}(t-t^{\prime})F^{cl}_{m}(t^{\prime})+\eta_{m}(t)\right)\Bigg{]}}_{\hbox{Eq.(1)}}$
$\displaystyle-\int
dt\;dt^{\prime}\sum_{n}\eta_{n}(t)[D^{K}]^{-1}(t-t^{\prime})\eta_{n}(t^{\prime})\Bigg{]}.$

This is equivalent to the TDVP Langevin of Eq.(6); the quantum field $A^{q}$
plays the role of a Lagrange multiplier that imposes Eq.(6) and the remaining
term gives the bath correlations. The tensor indices have been suppressed for
clarity in this expression. To make the comparison with the TDVP Langevin
equation, note that $\delta S[A^{cl}]/\delta A_{n}^{cl}(t)=0$ recovers the
usual TDVP equations for matrix product states. The additional terms
correspond to the dissipative effects of the bath. These terms are non-local
in the chain indices $n$ and $m$, despite our model of local independent
baths. This is due to the potential long-ranged entanglement of the matrix
product state. The long-range effects of the noise term reflect those already
found in the usual TDVP equations, since the noise term arises from a random
local potential. The non-locality of the friction term is more problematic and
some insight is required to implement it efficiently.

##  Appendix C Implementing the TDVP Langevin Equation

Here we outline how the TDVP Langevin equation, Eq.(1), can be implemented
numerically for MPS. This equation comprises three parts. The first is the
closed-system TDVP equation. This is implemented by standard means. The second
is the random noise induced by the environment. This is essentially a time-
dependent Hamiltonian term. We integrate it in a Stratonovich scheme. The
final part is friction. Even in the Markovian limit, this term is generally
spatially non-local — however, a significant simplification can be acheived by
working with purely local operators.

The TDVP equations for MPS can be written n the form

$\displaystyle\langle\partial_{i}\psi|\partial_{j}\psi\rangle\dot{X}_{j}$
$\displaystyle=$ $\displaystyle-i\langle\partial_{i}\psi|\hat{H}|\psi\rangle$
$\displaystyle\Rightarrow\;\;\;\;\langle\partial_{A_{n}}\psi|\partial_{A_{m}}\psi\rangle\dot{A}_{m}$
$\displaystyle=$
$\displaystyle-i\langle\partial_{A_{n}}\psi|\hat{H}|\psi\rangle$ (9)

where we have suppressed the tensor indices of $A$ for clarity, retaining only
the site index. The solution of this equation is well established. A judicious
choice of gauge fixing for the tangent vectors to the MPS manifold puts the
Gramm matrix $\langle\partial_{A_{n}}\psi|\partial_{A_{m}}\psi\rangle$ in a
diagonal form. Various algorithms for evaluating the TDVP equations for finite
systems exist, in this report we have used a modification of the method
introduced in haegeman2016unifying . A single time-step of the algorithm is
achieved by sweeping through the system from right to left and applying a
unitary rotation to the local variables on each site $A_{n}(t+\delta
t)=e^{iH_{eff}}A_{n}(t)$, followed by repeating this process by sweeping from
left to right.

The Noise  contribution to the TDVP Langevin equation for MPS can be written
in the form

$\displaystyle-i\sum_{m}\langle\partial_{i}\psi|\hat{F}_{m}|\psi\rangle\eta(t)$
(10) $\displaystyle\Rightarrow$
$\displaystyle-i\sum_{m}\langle\partial_{A_{n}}\psi|\hat{F}_{m}|\psi\rangle\eta_{m}(t)$

This evidently takes the same form as the right hand side of Eq.(9) and no
substantial modification to the TDVP algorithm is required. At the beginning
of each timestep, $\eta_{m}(t)$ is chosen by sampling from a normal
distribution with mean zero and variance $2\delta t\gamma T$.

Friction: The friction term can be written in the Markovian limit and in terms
of MPS tensors as follows:

$\displaystyle-i\sum_{m}\gamma\frac{\langle\psi|\hat{F}_{m}|\psi\rangle}{dt}\langle\partial_{i}\psi|\hat{F}_{m}|\psi\rangle,$
(11)

where

$\frac{\langle\psi|\hat{F}_{m}|\psi\rangle}{dt}=\langle\psi|\hat{F}_{m}|\partial_{A_{n}}\psi\rangle\dot{A}_{n}+\langle\partial_{\bar{A}_{n}}\psi|\hat{F}_{m}|\psi\rangle\dot{\bar{A}}_{n}.$
(12)

Eq.(12) can be evaluated by substituting in Eq.(9) introducing the poisson
bracket notation,
$\\{O_{1},O_{2}\\}=i\langle\psi|\hat{O}_{1}|\partial_{A_{n}}\psi\rangle\langle\partial_{\bar{A}_{n}}\psi|\hat{O}_{2}|\psi\rangle-i\langle\psi|\hat{O}_{2}|\partial_{A_{n}}\psi\rangle\langle\partial_{\bar{A}_{n}}\psi|\hat{O}_{1}|\psi\rangle$,

$\begin{split}&(\delta_{m,n}+\sum_{n}\gamma\\{F_{m},F_{n}\\})\frac{\langle\psi|\hat{F}_{n}|\psi\rangle}{dt}=\\{F_{m},H\\}+\sum_{n}\\{F_{m},F_{n}\\}\eta_{n}(t)\\\
&\rightarrow\frac{\langle\psi|\hat{F}_{n}|\psi\rangle}{dt}=(\mathbb{I}+\gamma\mathbb{F})^{-1}\left(\\{F_{m},H\\}+\sum_{n}\\{F_{m},F_{n}\\}\eta_{n}(t)\right)\end{split}$
(13)

where $\mathbb{F}_{ij}=\\{F_{i},F_{j}\\}$. Evaluating Eq.(13) exactly for
arbitrary operators $\hat{F}_{n}$ is quite numerically inefficient, scaling
quadratically in the number of noise operators. Moreover, it is inconsistent
with the site by site sweep algorithm introduced above for the TDVP equations.
Fortunately, Eq.(13) is substantially simplified in the case of single-site,
local noise fields. For local fields $\\{F_{i},F_{j}\\}$ is only non-zero
provided the two operators are located on the same site of the system and so
$\mathbb{F}$ becomes a simple, block diagonal matrix.

The modified TDVP algorithm therefore works as follows: Before each sweep
through the system $\eta_{m}(t)$ is sampled from a normal distribution with
mean $0$ and variance $2\delta t\gamma T$. The Hamiltonian and noise terms are
then used to calculate $\frac{\langle\psi|\hat{F}_{n}|\psi\rangle}{dt}$ for
all noise operators $F_{n}$ using Eq.(13). The noise and friction terms are
then combine with the Hamiltonian using Eq.(9) and the state is evolved using
the standard TDVP algorithm.