11institutetext: University of Pennsylvania
11email<EMAIL_ADDRESS>

# Extraction of Hierarchical Functional Connectivity Components in human brain
using Adversarial Learning

Dushyant Sahoo and Christos Davatzikos

###### Abstract

The estimation of sparse hierarchical components reflecting patterns of the
brain’s functional connectivity from rsfMRI data can contribute to our
understanding of the brain’s functional organization, and can lead to
biomarkers of diseases. However, inter-scanner variations and other
confounding factors pose a challenge to the robust and reproducible estimation
of functionally-interpretable brain networks, and especially to reproducible
biomarkers. Moreover, the brain is believed to be organized hierarchically,
and hence single-scale decompositions miss this hierarchy. The paper aims to
use current advancements in adversarial learning to estimate interpretable
hierarchical patterns in the human brain using rsfMRI data, which are robust
to “adversarial effects” such as inter-scanner variations. We write the
estimation problem as a minimization problem and solve it using alternating
updates. Extensive experiments on simulation and a real-world dataset show
high reproducibility of the components compared to other well-known methods.

###### Keywords:

fMRI analysis Matrix factorization Adversarial Learning.

##  1 Introduction

There has been a lot of research on estimating interpretable components of
functional connectivity of the brain. However, these components are often
vulnerable to confounding variations, herein referred to as “adversary” using
ML language, such as inter-scanner and inter-protocol variations, and rsfMRI
noise or irrelevant fluctuations, which can considerably reduce these
components’ reproducibility and hence their utility as biomarkers of diseases
that disrupt functional connectivity. To address this limitation, in this
paper, we introduce adversarial learning aiming to estimate components that
are robust to such confounding variations. Moreover, we use an existing
hierarchical decomposition approach [15], which attempts to capture the
brain’s multi-scale functional organization.

Numerous methods have been developed for estimating interpretable components
that capture the complexity of the human brain. Single-layer matrix
factorization approaches, such as ICA [17] and NMF [13], are not sufficient to
capture these complexities, as these components are believed to have
hierarchical organization [2]. In order to construct hierarchical components
that are interpretable, we extend previous work on hierarchical Sparse
Connectivity Patterns (hSCPs) [3, 15]. This approach’s advantages are that it
can extract sparse, overlapping, and hierarchical components, which are some
of the desired properties that can enhance the understanding of the human
brain’s working. It also extracts subject-specific information of the
components that is useful for capturing heterogeneity in the data.

We introduce the Adversarial hSCPs (Adv. hSCPs) method to enhance the sparse
component’s robustness, which can improve its generalization performance. We
formulate the problem as a bilevel matrix factorization problem and solve it
using alternate minimization. Our method is motivated by recent advances in
matrix factorization approaches that have used adversarial training to achieve
the state of the art performances [7, 9]. In a nutshell, it is a minimax game,
where the adversary perturbs model parameters to maximize or deteriorate our
objective function, and in defense, we minimize the objective function. We
discuss more about adversarial learning in the upcoming sections.

Outline: We start by reviewing hSCP and adversarial learning. Then, in Section
3, we present our method Adversarial hSCP. In Section 4 that follows, we
compare our method against existing methods on simulated and real data.
Additional experiments and calculations are in the Appendix.

Contributions: We propose adversarial learning for extracting hierarchical
Sparse Connectivity Components (hSCP). Empirical studies performed on
simulations and real-world datasets demonstrate that our model can generate
more reproducible components than other related methods. We also extract
interpretable components using the HCP dataset.

Notations and Conventions: We follow the same notations as used in [15].
$\mathbb{S}^{P\times P}_{++}$ is the set of symmetric positive definite
matrices of size $P\times P$. $\textbf{A}\geq 0$ denotes all the elements of
the matrix are greater than or equal to $0$. $\textbf{1}_{P}$ denotes $P\times
P$ matrix with all elements equal to one. $\textbf{I}_{P}$ denotes $P\times P$
identity matrix. $\mathbf{A}\circ\mathbf{B}$ is the element wise product
between two matrices $\mathbf{A}$ and $\mathbf{B}$.

##  2 PRELIMINARIES

###  2.1 Introduction to hierarchical Sparse Connectivity Patterns

Let $\mathbf{X}_{t}^{i}=[X_{t1}^{i},X_{t1}^{i},\ldots,X_{tP}^{i}]$ be the fMRI
BOLD responses of $i^{th}$ subject measured on $P$ nodes at time $t$, for
$t=1,\ldots,T$. The approach takes the correlation matrix for each subject
$\mathbf{\Theta}^{i}\in\mathbb{S}^{P\times P}_{++}$ as an input where each
correlation matrix stores pairwise correlation between the nodes. In hSCP,
$\mathbf{\Theta}^{i}$ is decomposed into non-negative linear combination of
sparse hierarchical components as:

$\displaystyle\mathbf{\Theta}^{i}$
$\displaystyle\approx\mathbf{W}_{1}\mathbf{\Lambda}_{1}^{i}\mathbf{W}_{1}^{\top}$
$\displaystyle\vdots$ $\displaystyle\mathbf{\Theta}^{i}$
$\displaystyle\approx\mathbf{W}_{1}\mathbf{W}_{2}\ldots\mathbf{W}_{K}\mathbf{\Lambda}^{i}_{K}\mathbf{W}_{K}^{\top}\mathbf{W}_{K-1}^{\top}\ldots\mathbf{W}_{1}^{\top},$

where $K$ is depth of hierarchy, $P>k_{1}>\ldots>k_{K}$, set
$\mathcal{W}={\\{\mathbf{W}_{r}\mid r=1,\ldots,K\\}}$ stores information about
sparse components shared across all subjects and set
$\mathcal{D}={\\{\mathbf{\Lambda}_{r}^{i}\mid r=1,\ldots,K;i=1\ldots,S\\}}$
stores subject specific diagonal matrix with $\mathbf{\Lambda}^{i}_{r}\geq 0$.
In the above formulation, $k_{r}$ denotes number of components at $r^{th}$
level. To understand the decomposition, we look at a two-layer hierarchy,
where $\mathbf{W}_{1}\in\mathbb{R}^{P\times k_{1}}$ stores $k_{1}$ components
at the first layer, which is the bottom-most level. If we look at second
level, which is the layer above the first, the components are comprised of
linear combination of components at first level i.e.
$\mathbf{W}_{1}\times\mathbf{W}_{2}\in\mathbb{R}^{P\times k_{2}}$, where
$k_{2}<k_{1}$. This repeated multiplication by low dimensional matrix will
result in the coarser components which are at the top of the hierarchy.
$L_{1}$ constraint on $\mathbf{W}_{s}$ is imposed to make $\mathbf{W}_{1}$
sparse resulting in sparse components and constraining rest of $\mathbf{W}$s
to force the components at next layer to comprise of a sparse non-negative
linear combination of the previous one. hSCP solves the below optimization to
estimate hierarchical patterns:

$\displaystyle\underset{\mathcal{W},\mathcal{D}}{\text{minimize}}$
$\displaystyle
H(\mathcal{W},\mathcal{D},\mathcal{C})=\sum_{i=1}^{S}\sum_{r=1}^{K}\|\mathbf{\Theta}^{i}-(\prod_{j=1}^{r}\mathbf{W}_{j})\mathbf{\Lambda}_{r}^{i}(\prod_{n=1}^{r}\mathbf{W}_{n})^{\top}\|_{F}^{2}$
(1) subject to
$\displaystyle\|\mathbf{w}^{r}_{l}\|_{1}<\lambda_{r},\;\|\mathbf{w}^{r}_{l}\|_{\infty}\leq
1,\;\mathop{\rm trace}\nolimits(\mathbf{\Lambda}_{r}^{i})=1,$
$\displaystyle\mathbf{\Lambda}_{r}^{i}\geq 0,\;\forall i,r,l;\hskip
20.00003pt\mathbf{W}_{j}\geq 0,\;j=2,\ldots,K,$

where $\mathcal{C}={\\{\mathbf{\Theta}^{i}\mid i=1,\ldots,S\\}}$,
$l=1,\ldots,k_{r},\;i=1,\ldots,S\;\text{and}\;r=1,\ldots,K$. $L_{1}$,
$L_{\infty}$ and $\mathop{\rm trace}\nolimits$ constraints help the problem to
identify a decomposition which can provide reproducible components. More
details can be found in the original paper [15].

###  2.2 Adversarial Training

A discriminative model is usually trained by minimizing the empirical expected
loss over a function class
$\mathcal{F}=\\{f_{\mathbf{v}}:\mathbf{v}\in\mathcal{V}\\}$:

$\displaystyle\min_{\mathbf{v}\in\mathcal{V}}\frac{1}{n}\sum_{i=1}^{n}l(f_{\mathbf{v}}(\mathbf{x}_{i}),y_{i}),$

where $l(\cdot,\cdot)$ is a loss function, $f_{\mathbf{v}}$ is the output
function, $\mathbf{x}_{i}$ is the feature vector and $y_{i}$ is the label.
Several recent papers [6, 11] have revealed that adding adversarial noise to
the sample defined for a sample $(\mathbf{x},y)$ as:

$\displaystyle\delta_{\mathbf{v}}^{adv}(\mathbf{x}):=\operatorname*{arg\,max}_{\|\delta\|\leq\epsilon}l(f_{\mathbf{v}}(\mathbf{x}+\delta),y),$
(2)

where $\epsilon>0$ is the adversarial noise power can drastically decrease
model’s performance. Adversarial training [10] was introduced to provide
robustness against the adversaries defined above. The training involves
empirical risk minimization over the perturbed samples by solving

$\displaystyle\min_{\mathbf{v}\in\mathcal{V}}\frac{1}{n}\sum_{i=1}^{n}l(f_{\mathbf{v}}(x_{i}+\delta_{\mathbf{v}}^{adv}(\mathbf{x})),y_{i}).$

The above formulation has a drawback that the accuracy drops drastically. To
overcome this problem, Mix-minibatch adversarial training (MAT) [21] is
performed by solving

$\displaystyle\min_{\mathbf{v}\in\mathcal{V}}\frac{1}{n}\sum_{i=1}^{n}l(f_{\mathbf{v}}(x_{i}+\delta_{\mathbf{v}}^{adv}(\mathbf{x})),y_{i})+l(f_{\mathbf{v}}(x_{i}),y_{i}),$
(3)

which balances between accuracy on the clean examples and robustness on the
adversarial examples. Motivated by the above methodology, we build adversarial
training for learning sparse hierarchical connectivity components. As the
above training regime was supervised, we use a different formulation with the
same idea for the unsupervised hSCP model.

##  3 Adversarial hierarchical Sparse Connectivity Patterns

In this section, we demonstrate how to incorporate adversarial learning at one
level for the hSCP method, which then can be extended to multiple levels. The
idea is to perturb input data $\mathbf{\Theta}^{i}$ and learn stable
components $\mathbf{W}_{1}$ which are robust to adversaries, such as inter-
scanner differences and unwanted rsfMRI noise. There are two parts of the
complete learning procedure-

Attack\- We first manually perturb the input data to get perturbed data
$\mathbf{\Gamma}^{i}=\mathbf{\Theta}^{i}+0.1\sigma\mathbf{1}_{P}$ where
$\sigma$ is the standard deviation of the data and use it to learn a perturbed
weight matrix $\mathbf{\tilde{W}}_{1}$ by minimizing the below cost function:

$\displaystyle
A(\mathbf{\hat{W}}_{1})=\alpha\|\mathbf{\tilde{W}}_{1}-\mathbf{W}_{1}\|_{F}^{2}+\sum_{i=1}^{S}\|\mathbf{\Gamma}^{i}-\mathbf{\tilde{W}}_{1}\mathbf{\Lambda}^{i}\mathbf{\tilde{W}}^{\top}_{1}\|_{F}^{2}.$
(4)

In the above equation, the first part is used to estimate
$\mathbf{\tilde{W}}_{1}$, which is close to $\mathbf{{W}}_{1}$ in Frobenius
norm, thus mimicking the actual components but is learned from the noise-
induced data. The second term is used for learning $\mathbf{\tilde{W}}_{1}$
using a perturbed data matrix. The main goal of the attacker is to learn
$\mathbf{\tilde{W}}_{1}$ for a given $\mathbf{\Gamma}^{i}$ and fool the model
by forcing the model to learn $\mathbf{\Lambda}^{i}$ from the perturbed data.
Our framework does not depend on the perturbations’ assumptions; the type of
perturbation can be varied depending on different types of practical noises
such as site-induced or scanner-induced noise. Finding optimal perturbation is
left for future work.

Defense\- Aim of the learner is to estimate $\mathbf{\Lambda}^{i}$ and
$\mathbf{W}_{1}$ by minimizing the below cost function:

$\displaystyle
D(\mathbf{W},\mathbf{\Lambda})=\sum_{i=1}^{S}\|\mathbf{\Theta}^{i}-\mathbf{\tilde{W}}_{1}\mathbf{\Lambda}^{i}\mathbf{\tilde{W}}^{\top}_{1}\|_{F}^{2}+\beta\sum_{i=1}^{S}\|\mathbf{\Theta}^{i}-\mathbf{W}_{1}\mathbf{\Lambda}^{i}\mathbf{{W}}^{\top}_{1}\|_{F}^{2},$
(5)

for a fixed $\mathbf{\tilde{W}}_{1}$. Learner first estimates subject specific
information $\mathbf{\Lambda}^{i}$ using perturbed weight matrix and then use
it to learn $\mathbf{{W}}_{1}$. We can now define the optimization problem for
the complete adversarial learning at single level using equation 4 and 5 as:

$\displaystyle\underset{\mathbf{W}_{1},\mathbf{\Lambda}^{i}\;\forall
i}{\text{minimize}}$
$\displaystyle\sum_{i=1}^{S}\|\mathbf{\Theta}^{i}-\mathbf{\tilde{W}}_{1}\mathbf{\Lambda}^{i}\mathbf{\tilde{W}}^{\top}_{1}\|_{F}^{2}+\beta\sum_{i=1}^{S}\|\mathbf{\Theta}^{i}-\mathbf{W}_{1}\mathbf{\Lambda}^{i}\mathbf{{W}}^{\top}_{1}\|_{F}^{2}$
(6) subject to
$\displaystyle\mathbf{\tilde{W}}_{1}=\operatorname*{arg\,min}_{\mathbf{\hat{W}}_{1}}\alpha\|\mathbf{\hat{W}}_{1}-\mathbf{{W}}_{1}\|_{F}^{2}+\sum_{i=1}^{S}\|\mathbf{\Gamma}^{i}-\mathbf{\tilde{W}}_{1}\mathbf{\Lambda}^{i}\mathbf{\tilde{W}}^{\top}_{1}\|_{F}^{2}$

The above equation is analogous to discriminate adversarial learning problem
defined in equation 3. Multi level formulation of the above problem can be
written as:

$\displaystyle\underset{\mathcal{W},\mathcal{D}}{\text{minimize}}$
$\displaystyle
J(\mathcal{\tilde{W}},\mathcal{W},\mathcal{D},\mathcal{C})=H(\mathcal{\tilde{W}},\mathcal{D},\mathcal{C})+\beta
H(\mathcal{W},\mathcal{D},\mathcal{C})$ (7) subject to
$\displaystyle\mathbf{\tilde{W}}_{r}=\operatorname*{arg\,min}_{\mathbf{\hat{W}}_{r}}\alpha\|\mathbf{\hat{W}}_{r}-\mathbf{{W}}_{r}\|_{F}^{2}+H(\mathcal{\tilde{W}},\mathcal{D},\mathcal{P})\quad
r=1,\ldots,K$

where $\mathcal{\tilde{W}}={\\{\mathbf{\tilde{W}}_{r}\mid r=1,\ldots,K\\}}$
and $\mathcal{P}={\\{\mathbf{\Gamma}^{i}\mid i=1,\ldots,S\\}}$.

Algorithm 1 Adv. hSCP

1:Input: Data $\mathcal{C}$, perturbed data $\mathcal{P}$; $\mathcal{W}$ and
$\mathcal{D}$ = hSCP($\mathcal{C}$)

2:repeat

3: for $r=1$ to $K$ do

4: Update adversarial perturbations

5: $\mathbf{\hat{W}}_{r}\leftarrow\mathop{\rm
descent}\nolimits(\mathbf{\hat{W}}_{r})$

6: Update model parameters

7: $\mathbf{W}_{r}\leftarrow\mathop{\rm descent}\nolimits(\mathbf{W}_{r})$

8: if $r==1$ then

9: $\mathbf{W}_{r}\leftarrow\mathop{\rm proj}\nolimits_{1}(\mathbf{W}_{r})$

10: else

11: $\mathbf{W}_{r}\leftarrow\mathop{\rm proj}\nolimits_{2}(\mathbf{W}_{r})$

12: $\mathbf{\Lambda}_{r}^{i}\leftarrow\mathop{\rm
descent}\nolimits(\mathbf{\Lambda}_{r}^{i})$;
$\mathbf{\Lambda}_{r}^{i}\leftarrow\mathop{\rm
proj}\nolimits_{2}(\mathbf{\Lambda}_{r}^{i})\quad i=1,\ldots,S$

13:until Stopping criterion is reached

14:Output: $\mathbf{W}$ and $\mathbf{\Lambda}$

###  3.1 Optimization

The complete algorithm to solve the above optimization problem is described in
Algorithm 1 (Adv. hSCP). First, the adversarial perturbations are generated by
performing gradient descent on $\mathbf{\hat{W}}_{r}$, and then the model
parameters are updated using gradient descent. This process is repeated until
the convergence criteria is reached. $\mathop{\rm descent}\nolimits$ is the
update rules defined by AMSgrad [14] for performing gradient descent. The
update rule and gradients are defined in Appendix. $\mathop{\rm
proj}\nolimits_{1}(\mathbf{A})$ [12] function projects each column of
$\mathbf{A}$ into intersection of $L_{1}$ and $L_{\infty}$ ball defined in
equation 1 and $\mathop{\rm proj}\nolimits_{2}(\mathbf{A})$ function makes all
the negative elements of $\mathbf{A}$ equal to zero. The model parameters are
initialized by first optimizing the hSCP model using Algorithm 1 (hSCP)
defined in [15], rather than randomly initialized. This makes algorithm
deterministic, and the algorithm can start from an optimal point on which
adversarial learning can improve if there is overfitting.

##  4 Experiments

###  4.1 Comparison with existing methods

####  4.1.1 Simulated dataset

Table 1: Accuracy on simulated dataset Method | $k_{1}=6$ | $k_{1}=8$ | $k_{1}=10$ | $k_{1}=12$
---|---|---|---|---
hSCP |  $0.801$ |  $0.829$ |  $0.818$ |  $0.814$
Adv. hSCP |  $0.832$ |  $0.847$ |  $0.867$ |  $0.864$
ICA |  $0.656\pm 0.004$ |  $0.696\pm 0.022$ |  $0.734\pm 0.025$ |  $0.748\pm 0.011$
NMF |  $0.650\pm 0.104$ |  $0.701\pm 0.071$ |  $0.708\pm 0.118$ |  $0.712\pm 0.085$
Adv. NMF |  $0.695\pm 0.047$ |  $0.718\pm 0.069$ |  $0.720\pm 0.091$ |  $0.723\pm 0.114$
Table 2: Accuracy on simulated data with two level hierarchy |  | $k_{1}=6$ | $k_{1}=8$ | $k_{1}=10$ | $k_{1}=12$
---|---|---|---|---|---
$k_{2}=4$ |  hSCP |  $0.821$ |  $0.837$ |  $0.827$ |  $0.819$
Adv. hSCP |  $0.859$ |  $0.864$ |  $0.846$ |  $0.823$
$k_{2}=6$ |  hSCP |  $0.816$ |  $0.819$ |  $0.813$ |  $0.805$
Adv hSCP |  $0.848$ |  $0.849$ |  $0.826$ |  $0.814$

We first use a simulated dataset to evaluate the performance of our model
against SCP [3], NMF [13], adv. NMF [9] and ICA [17]. The values of $\alpha$
and $\beta$ are set to be $10^{-3}$ and $0.5$ respectively throughout the
paper. We generate sparse components $\mathbf{S}_{1}\in\mathbb{R}^{P\times
k_{1}}$ with $P=50$ and $k_{1}=8$ and generate network structure from it. This
network is then used as input to NetSim [18] with TR equal to $3$ seconds to
generate time-series data of $100$ subjects, each having $300$ time-points.
NetSim also adds Gaussian noise to the time series of each node. We also add
Poisson noise with a mean equal to $0.4$ to check how different methods
perform in a high noise scenario.

We compare components/factors derived from all the models with
$k_{1}\in\\{6,8,10,12\\}$. Accuracy is used as a performance measure defined
as a normalized inner product between ground truth components and estimated
factors derived from various algorithms. The optimal sparsity parameter
$\lambda_{1}$ in the hSCP is selected from $P\times 10^{[-2:1]}$ having the
highest average split-sample reproducibility in $20$ runs. Split-sample
reproducibility of components is computed by randomly dividing the data into
two equal parts and then calculating normalized inner product between
components extracted from each sample. A high reproducibility value implies
that the same component can be extracted from multiple samples. Table 1 shows
the comparison of the accuracy of different methods averaged over $20$ runs.
As the hSCP method is deterministic, the output remains the same in every run.
From the table, it can be seen that adversarial training can significantly
improve the accuracy of hSCP. An important thing to note here is that
adversarial training has also improved the accuracy of NMF, but it remains
less than that of hSCP. The result for Poisson noise case is presented in
Appendix.

We next generated a two-level hierarchy using the components defined above as
the first layer. We used linear operator for projection to lower dimensional
space to get coarse components with $P=50$ and $k_{2}=4$. Visualization of the
components is in Appendix. Time-series data were then generated under the same
settings presented above. Table 2 shows the components’ average accuracy at
two-level over $20$ runs for hSCP and Adv. hSCP with Adv. hSCP method giving
better results. We did not use ICA and NMF as they can only generate
components at only one level.

####  4.1.2 Resting state fMRI data

We used $100$ unrelated subjects released within the $900$ subjects data
release from the publicly available Human Connectome Project (HCP) [20]
dataset for comparing different methods. ICA+FIX pipeline [5] is used to
process the complete data. Each subject has $4$ scans, with each scan
comprising $1001$ time points and $360$ nodes.

We compare components/factors derived from all the models with
$k_{1}\in\\{5,10,15,20\\}$. As the ground truth is not known, we use split-
sample reproducibility as a performance measure. We first find the optimal
value of $\lambda_{1}$ from $P\times 10^{[-2:1]}$. Training, validation, and
test data are generated by dividing data equally into three parts, and we then
select $\lambda_{1}$ with the highest mean reproducibility over $20$ runs on
training and test data. Training and test are used for final reproducibility
comparison. Table 3 shows that the hSCP method can extract components with
high reproducibility. We have similar results presented in Table 4 for a two-
level decomposition.

Table 3: Reproducibility on HCP dataset Method |  $10$ |  $15$ |  $20$ |  $25$
---|---|---|---|---
hSCP |  $0.749\pm 0.045$ |  $0.750\pm 0.046$ |  $0.712\pm 0.026$ |  $0.701\pm 0.019$
Adv. hSCP |  $0.787\pm 0.052$ |  $0.765\pm 0.059$ |  $0.716\pm 0.020$ |  $0.721\pm 0.016$
ICA |  $0.695\pm 0.067$ |  $0.638\pm 0.046$ |  $0.581\pm 0.039$ |  $0.523\pm 0.027$
NMF |  $0.689\pm 0.038$ |  $0.657\pm 0.067$ |  $0.635\pm 0.053$ |  $0.629\pm 0.020$
Adv. NMF |  $0.709\pm 0.073$ |  $0.659\pm 0.043$ |  $0.653\pm 0.026$ |  $0.633\pm 0.032$
Table 4: Reproducibility on HCP dataset with two level hierarchy |  | $k_{1}=10$ | $k_{1}=15$ | $k_{1}=20$ | $k_{1}=25$
---|---|---|---|---|---
$k_{2}=4$ |  hSCP |  $0.872\pm 0.044$ |  $0.853\pm 0.064$ |  $0.831\pm 0.075$ |  $0.826\pm 0.091$
Adv. hSCP |  $0.895\pm 0.030$ |  $0.866\pm 0.029$ |  $0.848\pm 0.056$ |  $0.830\pm 0.061$
$k_{2}=6$ |  hSCP |  $0.856\pm 0.070$ |  $0.842\pm 0.062$ |  $0.828\pm 0.031$ |  $0.824\pm 0.035$
Adv hSCP |  $0.877\pm 0.076$ |  $0.864\pm 0.067$ |  $0.843\pm 0.045$ |  $0.834\pm 0.048$

###  4.2 Results from rsfMRI data

We extract $10$ components at level $1$, and $4$ components at level $2$ using
Adv. hSCP learning from the HCP dataset. Figure 1 shows two hierarchical
components. Component $1$ stores anti-correlation information between Default
Mode Network and Dorsal Attention Network previously studied using seed-based
correlation method [4]. Component $2$ stores anti-correlation between Default
Mode Network and extrastriate visual areas, which is another well-known
finding [19]. A more thorough discussion is needed for examining the
differences and similarities between the components derived from hSCP and Adv.
hSCP, which we have left for future work.

0-1=dj0-1= 1
---
0-1=
---
0-1=
---

0-1=dj0-1= 2
---
0-1=
---
0-1=
---

Figure 1:  Hierarchical components estimated using Adv. hSCP. Red and blue
color are used for showing negative correlations between regions in a
component.

##  5 Conclusion

In this study, we used adversarial learning to enhance the hSCP method by
increasing the hierarchical components’ reproducibility. We formulate the
problem as a bilevel optimization problem and used adaptive gradient descent
to solve it. Experimental results based on simulated data show that Adv hSCP
can extract components accurately compared to other methods. Results using
real-world rsfMRI data demonstrate the adversarial learning can improve the
reproducibility of the components. We also discuss the interpretability of the
components extracted from the HCP dataset.

There are several applications of this work. Improved reproducibility of the
components can increase accuracy and confidence when applied to clinical
applications such as age prediction, disease diagnosis, etc. Adversarial
learning can be extended to other matrix factorization approaches used for the
analysis of fMRI data, such as dynamic sparse connectivity patterns [1],
sparse granger causality patterns [16], deep non-negative matrix factorization
[8], etc. It would be interesting to assess the impact of the method in
characterizing activity in terms of task-induced activations.

## References

  * [1] Cai, B., Zille, P., Stephen, J.M., Wilson, T.W., Calhoun, V.D., Wang, Y.P.: Estimation of dynamic sparse connectivity patterns from resting state fmri. IEEE transactions on medical imaging 37(5), 1224–1234 (2017)
  * [2] Doucet, G., Naveau, M., Petit, L., Delcroix, N., Zago, L., Crivello, F., Jobard, G., Tzourio-Mazoyer, N., Mazoyer, B., Mellet, E., et al.: Brain activity at rest: a multiscale hierarchical functional organization. Journal of neurophysiology 105(6), 2753–2763 (2011)
  * [3] Eavani, H., Satterthwaite, T.D., Filipovych, R., Gur, R.E., Gur, R.C., Davatzikos, C.: Identifying sparse connectivity patterns in the brain using resting-state fmri. Neuroimage 105, 286–299 (2015)
  * [4] Fox, M.D., Snyder, A.Z., Vincent, J.L., Corbetta, M., Van Essen, D.C., Raichle, M.E.: The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proceedings of the National Academy of Sciences 102(27), 9673–9678 (2005)
  * [5] Glasser, M.F., Sotiropoulos, S.N., Wilson, J.A., Coalson, T.S., Fischl, B., Andersson, J.L., Xu, J., Jbabdi, S., Webster, M., Polimeni, J.R., et al.: The minimal preprocessing pipelines for the human connectome project. Neuroimage 80, 105–124 (2013)
  * [6] Goodfellow, I.J., Shlens, J., Szegedy, C.: Explaining and harnessing adversarial examples. arXiv preprint arXiv:1412.6572 (2014)
  * [7] He, X., He, Z., Du, X., Chua, T.S.: Adversarial personalized ranking for recommendation. In: The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval. pp. 355–364 (2018)
  * [8] Li, H., Zhu, X., Fan, Y.: Identification of multi-scale hierarchical brain functional networks using deep matrix factorization. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 223–231. Springer (2018)
  * [9] Luo, L., Zhang, Y., Huang, H.: Adversarial nonnegative matrix factorization. In: International Conference on Machine Learning. pp. 6479–6488. PMLR (2020)
  * [10] Madry, A., Makelov, A., Schmidt, L., Tsipras, D., Vladu, A.: Towards deep learning models resistant to adversarial attacks. In: International Conference on Learning Representations (2018)
  * [11] Moosavi-Dezfooli, S.M., Fawzi, A., Fawzi, O., Frossard, P.: Universal adversarial perturbations. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 1765–1773 (2017)
  * [12] Podosinnikova, A., Hein, M., Gemulla, R.: Robust Principal Component Analysis as a Nonlinear Eigenproblem. Ph.D. thesis, Saarland University (2013)
  * [13] Potluru, V.K., Calhoun, V.D.: Group learning using contrast nmf: Application to functional and structural mri of schizophrenia. In: 2008 IEEE International Symposium on Circuits and Systems. pp. 1336–1339. IEEE (2008)
  * [14] Reddi, S.J., Kale, S., Kumar, S.: On the convergence of adam and beyond. arXiv preprint arXiv:1904.09237 (2019)
  * [15] Sahoo, D., Satterthwaite, T.D., Davatzikos, C.: Hierarchical extraction of functional connectivity components in human brain using resting-state fmri. IEEE Transactions on Medical Imaging pp. 1–1 (2020). https://doi.org/10.1109/TMI.2020.3042873
  * [16] Sahoo, D., Honnorat, N., Davatzikos, C.: Gpu accelerated extraction of sparse granger causality patterns. In: 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018). pp. 604–607. IEEE (2018)
  * [17] Smith, S.M., Fox, P.T., Miller, K.L., Glahn, D.C., Fox, P.M., Mackay, C.E., Filippini, N., Watkins, K.E., Toro, R., Laird, A.R., et al.: Correspondence of the brain’s functional architecture during activation and rest. Proceedings of the national academy of sciences 106(31), 13040–13045 (2009)
  * [18] Smith, S.M., Miller, K.L., Salimi-Khorshidi, G., Webster, M., Beckmann, C.F., Nichols, T.E., Ramsey, J.D., Woolrich, M.W.: Network modelling methods for fmri. Neuroimage 54(2), 875–891 (2011)
  * [19] Uddin, L.Q., Clare Kelly, A., Biswal, B.B., Xavier Castellanos, F., Milham, M.P.: Functional connectivity of default mode network components: correlation, anticorrelation, and causality. Human brain mapping 30(2), 625–637 (2009)
  * [20] Van Essen, D.C., Ugurbil, K., Auerbach, E., Barch, D., Behrens, T., Bucholz, R., Chang, A., Chen, L., Corbetta, M., Curtiss, S.W., et al.: The human connectome project: a data acquisition perspective. Neuroimage 62(4), 2222–2231 (2012)
  * [21] Wong, E., Kolter, Z.: Provable defenses against adversarial examples via the convex outer adversarial polytope. In: International Conference on Machine Learning. pp. 5286–5295. PMLR (2018)

##  6 Appendix

(a) $k=10$

(b) $k=20$

Figure 2: Convergence of Adv. hSCP algorithm using HCP dataset for different
values of $k$. In the figure, it can be seen that initially, the function
value is low because the initial value is an optimal value of $\mathbf{W}$ and
$\mathbf{\Lambda}$ returned using the hSCP algorithm. As the adversarial
attack begins, the objective function value starts to fluctuate because of the
minimax game, where the adversarial perturbation tries to deviate the result
from the optimal value. In defense, we try to minimize the objective function.
The algorithm converges when the optimal value becomes robust to
perturbations.

###  6.1 Gradients

In this section, we define gradients used during adaptive gradient descent.
Let

$\displaystyle\mathbf{\tilde{W}}_{0}$
$\displaystyle=\mathbf{W}_{0}=\mathbf{I}_{P},\qquad\mathbf{Y}_{r}=\prod_{j=0}^{r}\mathbf{W}_{j},\qquad\mathbf{\tilde{Y}}_{r}=\prod_{j=0}^{r}\mathbf{\tilde{W}}_{j},$
$\displaystyle\mathbf{T}_{n,i}^{r}$
$\displaystyle=(\prod_{j=1}^{n-r}\mathbf{W}_{j})\mathbf{\Lambda}^{i}_{n-r}(\prod_{j=1}^{n-r}\mathbf{W}_{j})^{\top},\qquad\mathbf{\tilde{T}}_{n,i}^{r}=(\prod_{j=1}^{n-r}\mathbf{\tilde{W}}_{j})\mathbf{\Lambda}^{i}_{n-r}(\prod_{j=1}^{n-r}\mathbf{\tilde{W}}_{j})^{\top}.$

We first define gradient for updating adversarial perturbations
$\mathbf{\tilde{W}_{r}}$. The objective function is
$F=\alpha\|\mathbf{\hat{W}}_{r}-\mathbf{{W}}_{r}\|_{F}^{2}+H(\mathcal{\tilde{W}},\mathcal{D},\mathcal{P})$
and gradient with respect to $\mathbf{\tilde{W}_{r}}$ will be

$\displaystyle\frac{F}{\partial\mathbf{\tilde{W}_{r}}}$
$\displaystyle=2\alpha(\mathbf{\hat{W}}_{r}-\mathbf{{W}}_{r})+\frac{\partial
H(\mathcal{\tilde{W}},\mathcal{D},\mathcal{C})}{\partial\mathbf{\tilde{W}}_{r}}=2\alpha(\mathbf{\hat{W}}_{r}-\mathbf{{W}}_{r})$
$\displaystyle+\sum_{i=1}^{S}\sum_{j=r}^{K}-4\mathbf{\tilde{Y}}_{r-1}^{\top}\mathbf{\Gamma}_{i}\mathbf{\tilde{Y}}_{r-1}\mathbf{\tilde{W}}_{r}\mathbf{\tilde{T}}_{j,i}^{r}+4\mathbf{\tilde{Y}}_{r-1}^{\top}\mathbf{\tilde{Y}}_{r-1}\mathbf{\tilde{W}}_{r}\mathbf{\tilde{T}}_{j,i}^{r}\mathbf{\tilde{W}}_{r}^{\top}\mathbf{\tilde{Y}}_{r-1}^{\top}\mathbf{\tilde{Y}}_{r-1}\mathbf{\tilde{W}}_{r}\mathbf{\tilde{T}}_{j,i}^{r}.$

We now define gradients for updating model parameters. The gradient of
objective function $J$ with respect to $\mathbf{\Lambda}_{r}^{i}$ is:

$\displaystyle\frac{\partial
J}{\partial\mathbf{\Lambda}_{r}^{i}}=\frac{\partial
H(\mathcal{\tilde{W}},\mathcal{D},\mathcal{C})}{\partial\mathbf{\Lambda}_{r}^{i}}+\frac{\partial
H(\mathcal{{W}},\mathcal{D},\mathcal{C})}{\partial\mathbf{\Lambda}_{r}^{i}}$
$\displaystyle=[\beta(-2\mathbf{\tilde{Y}}_{r}^{\top}\mathbf{\Theta}_{r}^{i}\mathbf{\tilde{Y}}_{r}+2\mathbf{\tilde{Y}}_{r}^{\top}\mathbf{\tilde{Y}}_{r}\mathbf{\Lambda}_{r}^{i}\mathbf{\tilde{Y}}_{r}^{\top}\mathbf{\tilde{Y}}_{r})$
$\displaystyle+\beta(-2\mathbf{Y}_{r}^{\top}\mathbf{\Theta}_{r}^{i}\mathbf{Y}_{r}+2\mathbf{Y}_{r}^{\top}\mathbf{Y}_{r}\mathbf{\Lambda}_{r}^{i}\mathbf{Y}_{r}^{\top}\mathbf{Y}_{r})]\circ\mathbf{I}_{k_{r}}.$

The gradient of $J$ with respect to $\mathbf{W}_{r}$ is:

$\displaystyle\frac{\partial J}{\partial\mathbf{W}_{r}}=\frac{\partial
H(\mathcal{{W}},\mathcal{D},\mathcal{C})}{\partial\mathbf{W}_{r}}$
$\displaystyle=\sum_{i=1}^{S}\sum_{j=r}^{K}-4\mathbf{Y}_{r-1}^{\top}\mathbf{\Theta}_{i}\mathbf{Y}_{r-1}\mathbf{W}_{r}\mathbf{T}_{j,i}^{r}$
$\displaystyle+4\mathbf{Y}_{r-1}^{\top}\mathbf{Y}_{r-1}\mathbf{W}_{r}\mathbf{T}_{j,i}^{r}\mathbf{W}_{r}^{\top}\mathbf{Y}_{r-1}^{\top}\mathbf{Y}_{r-1}\mathbf{W}_{r}\mathbf{T}_{j,i}^{r}.$

###  6.2 AMSgrad update rule

Let $g_{i}$ be the partial derivative of the objective function with respect
to the parameter $w_{i}$ at $i^{th}$ iteration. Let $m_{i}$ and $v_{i}$ denote
the decaying averages of past and past squared gradients, then the update rule
for AMSgrad is defined as:

$\displaystyle m_{i}$
$\displaystyle=\beta_{1}m_{i-1}+(1-\beta_{1})g_{i},\qquad
v_{i}=\beta_{2}v_{i-1}+(1-\beta_{2})g_{i}^{2}$ $\displaystyle\hat{m}_{i}$
$\displaystyle=\frac{m_{i}}{1-\beta_{1}^{i}}\qquad\hat{v}_{i}=\text{max}(\hat{v}_{i-1},v_{i})\qquad
w_{i+1}=w_{i}-\frac{\eta}{\sqrt{\hat{v}_{i}}+\epsilon}\hat{m}_{i},$

where $\beta_{1}=0.9$, $\beta_{2}=0.99$, $\epsilon=10^{-8}$ and $\eta=0.1$.
$\beta_{1}$ and $\beta_{2}$ are the hyperparameters in the update rules
described above. These are typical values for the practical applications [14].

(a)

$\times$ =

(b)

(c)

Figure 3: (a) Visualization of ground truth components at level $1$. (b) Weight matrix used to generate components at level 2. (c) Visualization of ground truth components at level $2$.  Table 5: Accuracy on simulated dataset with Pois(0.4) noise added Method | $k_{1}=6$ | $k_{1}=8$ | $k_{1}=10$ | $k_{1}=12$
---|---|---|---|---
hSCP |  $0.798$ |  $0.779$ |  $0.739$ |  $0.724$
Adv. hSCP |  $0.804$ |  $0.771$ |  $0.818$ |  $0.843$
ICA |  $0.637\pm 0.015$ |  $0.671\pm 0.034$ |  $0.715\pm 0.027$ |  $0.738\pm 0.012$
NMF |  $0.640\pm 0.101$ |  $0.655\pm 0.109$ |  $0.703\pm 0.079$ |  $0.704\pm 0.128$
Adv. NMF |  $0.690\pm 0.079$ |  $0.681\pm 0.071$ |  $0.694\pm 0.080$ |  $0.682\pm 0.088$


On the stability of the infinite Projected Entangled Pair Operator ansatz for
driven-dissipative 2D lattices

D. Kilda1*, A. Biella2, M. Schiró2,3, R. Fazio4,5, J. Keeling6

1 Division of Chemistry and Chemical Engineering, California Institute of
Technology, Pasadena, CA, USA

2 JEIP, USR 3573 CNRS, Collège de France, PSL Research University, 11 Place
Marcelin Berthelot, 75321 Paris Cedex 05, France

3 Institut de Physique Théorique, Université Paris Saclay, CNRS, CEA, 91191
Gif-sur-Yvette, France

4 Abdus Salam ICTP, Strada Costiera 11, I-34151 Trieste, Italy

5 Dipartimento di Fisica, Università di Napoli ’Federico II’, Monte S. Angelo,
I-80126 Napoli, Italy

6 SUPA, School of Physics and Astronomy, University of St Andrews, St Andrews,
KY16 9SS, United Kingdom

*<EMAIL_ADDRESS>
## Abstract

We present calculations of the time-evolution of the driven-dissipative XYZ
model using the infinite Projected Entangled Pair Operator (iPEPO) method,
introduced by [A. Kshetrimayum, H. Weimer and R. Orús, Nat. Commun. 8, 1291
(2017)]. We explore the conditions under which this approach reaches a steady
state. In particular, we study the conditions where apparently converged
calculations may become unstable with increasing bond dimension of the tensor-
network ansatz. We discuss how more reliable results could be obtained.

###### Contents

  1. 1 Introduction
  2. 2 Application to the dissipative XYZ model
    1. 2.1 Effects of simulation protocol
    2. 2.2 Effects of bond dimension
  3. 3 Conclusion
  4. A Implementation of iPEPO algorithm
    1. A.1 Summary of the iPEPS algorithm
      1. A.1.1 Time evolution: simple update
      2. A.1.2 Contraction: corner transfer matrix
    2. A.2 Extension to iPEPO
    3. A.3 Benchmarking with ground state calculations

##  1 Introduction

Tensor network approaches have provided a route to efficient numerical
simulations across a wide range of physical problems [1, 2, 3, 4]. In one
dimension, matrix product states (MPS) were originally introduced in the
context of one-dimensional quantum ground states [5, 6]. They have
subsequently been extended to finite temperatures [7, 8], and to open quantum
systems and density-matrix evolution [9, 10]. Such methods have been very
fruitful in exploring the nonequilibrium steady states (NESS) of driven-
dissipative one-dimensional systems, using matrix product operators (MPO) [11,
12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. These methods can also be extended
beyond one dimension, either by mapping a finite two-dimensional lattice onto
a one-dimensional chain [22]—see Ref. [21] for a driven-dissipative
implementation—or via the projected entangled pair state (PEPS) algorithm [23,
24, 25, 26]. The PEPS approach represents the two-dimensional lattice directly
as a tensor network [3, 4], and allows a direct simulation of an infinite
(translationally invariant) lattice (iPEPS).

In a significant development, Kshetrimayum et al. [27] presented results
adapting the iPEPS algorithm to simulate open quantum systems on infinite 2D
lattices. By using an infinite projected entangled pair operator (iPEPO)
algorithm, they calculated the NESS of the dissipative XYZ and transverse
field Ising models — i.e. finding the steady state of a many-body Lindblad
master equation. The ability to routinely apply such methods to two-
dimensional open quantum systems is potentially very powerful. While one-
dimensional systems have shown a rich variety of collective behaviour,
symmetry-breaking phase transitions generally do not occur in open one-
dimensional systems, while they can in two dimensions. Several alternative
approaches to approximately simulate two-dimensional open systems have been
proposed, including cluster mean field theory [18], corner space
renormalization [28, 29], and neural network states [30, 31, 32, 33, 34].
However, so far, these methods have generally been restricted to small systems
(or small clusters), making it challenging to extract critical behavior. As
such, the ability to routinely use iPEPO could be extremely powerful to
numerically explore phase transitions and critical behavior in driven
dissipative systems.

Here, we explore in detail the stability of the iPEPO algorithm, introduced by
Kshetrimayum et al [27]. We find that while at short times the algorithm shows
reasonable time evolution, the behaviour at long times varies. In particular,
we find that the algorithm only reaches a steady state in some parameter
regimes, and close to dissipative critical points [35] it can fail to reach a
steady state. The regimes where we fail to find a steady state correspond
closely to the regimes where Ref. [27] found a larger value of their parameter
$\Delta$, which measures how close the state they find is to a steady state.
Moreover, we find that for some parameters, increasing bond dimension of the
iPEPO representation does not systematically improve the accuracy of the
results. On the contrary, it can in some cases destabilize a fixed point
obtained at a lower bond dimension. We also suggest some possible alternatives
to the simple-update iPEPO algorithm, which could help to alleviate the
problem.

Our paper is organised as follows. In Sec. 2 we apply the iPEPO algorithm to
calculate NESS of the dissipative XYZ model in 2D, and analyse whether a
steady state can be found. Section. 3 concludes with some comments on
alternative tensor network approaches for computing NESS in 2D. We also
provide an extended appendix, Sec. A, which summarises our implementation of
the iPEPO algorithm; this implementation can be found at [36]. Since the core
of the iPEPO and iPEPS algorithms is similar, we also present results
benchmarking our implementation against the prototypical applications of
iPEPS: the ground states of the transverse field Ising model and the hardcore
Bose-Hubbard model.

##  2 Application to the dissipative XYZ model

In this section, using the iPEPO implementation described and benchmarked in
Appendix A, we discuss finding the NESS of the dissipative spin-$1/2$ XYZ
model on an infinite square lattice. The specific density-matrix equation of
motion that we consider is:

$\displaystyle\partial_{t}\rho$
$\displaystyle=-i\left[H_{XYZ},\rho\right]+\frac{\kappa}{2}\sum_{j}\left(2\sigma_{j}^{-}\,\rho\sigma_{j}^{+}-\sigma_{j}^{+}\sigma_{j}^{-}\,\rho-\rho\,\sigma_{j}^{+}\sigma_{j}^{-}\right),$
(1) $\displaystyle H_{XYZ}$ $\displaystyle=\sum_{\left\langle
i,j\right\rangle}\left(J_{x}\sigma^{x}_{i}\sigma^{x}_{j}+J_{y}\sigma^{y}_{i}\sigma^{y}_{j}+J_{z}\sigma^{z}_{i}\sigma^{z}_{j}\right),$
(2)

where $\sigma^{x,y,z}_{j}$ are Pauli matrices at lattice site $j$,
$\sigma^{\pm}=\frac{1}{2}\left(\sigma^{x}\pm i\sigma^{y}\right)$, and
$J_{x,y,z}$ are spin-spin coupling constants.

Figure 1: The magnetization order parameter $M_{x}$ of the dissipative XYZ
model as a function of coupling strength $J_{y}$, for $J_{x}=0.5$, $J_{z}=1$
and $D=4$. Energies given in units of $\kappa=1$. The red highlighted areas
indicate parameter regimes where the iPEPO algorithm fails to reach a steady
state. (a) Results computed using timesteps $\delta t=10^{-1},10^{-2}$ run
until a steady state is found. (In the red regions, the run is stopped after
$N=1000$ steps with $\delta t=10^{-1}$ followed by $N=2000$ steps with $\delta
t=10^{-2}$.) (b) Results calculated using a large timestep $\delta t=10^{-1}$
and stopping after $N=1000$ steps.

We begin by computing the time evolution of the dissipative XYZ model for the
same parameters considered by Ref. [27]. Figure 1(a) shows magnetization
averaged over the two sites $i=A,B$ of iPEPO unit cell,
$M_{x}=\frac{1}{2}(|\expectationvalue{\sigma^{x}_{i=A}}|+|\expectationvalue{\sigma^{x}_{i=B}}|)$,
as a function of coupling strength $J_{y}$ using iPEPO bond dimension $D=4$.
We find that iPEPO algorithm only reaches a steady state for some values of
$J_{y}$, while in the red highlighted areas no steady state is found—the
results continue to change with time. Where a steady state is found, our
results closely match Ref. [27]. Figure 1(b) shows that if we use a large
timestep, $\delta t=10^{-1}$, and deliberately stop the simulation early —
i.e. after $N=1000$ steps — one can reproduce results similar to those
presented by Kshetrimayum et al [27]. However, these results do not correspond
to an actual steady state, and also inevitably contain significant Trotter
errors due to the large timestep size. As we discuss further below, the
failure to reach a steady state that we observe occurs specifically in the SU
time evolution. That is, it is completely unaffected by the CTM contractions
needed to compute observables.

We have encountered similar issues of failing to reach a steady state in other
parameter regimes of the dissipative XYZ model, as well as for other systems
such as the dissipative transverse field Ising model in 2D. This raises
important question about the practicality of the iPEPO algorithm as a tool to
find the NESS of open quantum systems. In the following, to keep our
discussion concise, we will restrict our attention to the dissipative XYZ
model. Our goal below will be to understand when the iPEPO algorithm does and
does not reach a steady state, focusing entirely on the SU time evolution, by
studying the relative change of singular values $\epsilon_{\Lambda}$ defined
in Eq. (5).

###  2.1 Effects of simulation protocol

Figure 2: The evolution of $\epsilon_{\Lambda}$ (from $\Lambda^{[U]}$) at
$J_{y}=1.2$ (where no steady state is found) and $J_{y}=1.5$ (where a steady
state is found), using timesteps (a) $\delta t=10^{-1}$, (b) $\delta
t=10^{-2}$, (c) $\delta t=10^{-3}$. Energies given in units of $\kappa=1$.

In Figs. 2(a-c) we plot $\epsilon_{\Lambda}$ against simulation step, for
timestep sizes $\delta t=10^{-1},10^{-2},10^{-3}$ respectively. We show both
$J_{y}=1.2$, in the parameter regime where iPEPO fails to reach a steady state
(blue line), and $J_{y}=1.5$, where a steady state is found (green line). At
$J_{y}=1.5$, we observe clearly that $\epsilon_{\Lambda}$ quickly decreases,
indicating that we approach a steady state. However, at $J_{y}=1.2$ we see
$\epsilon_{\Lambda}$ undergoes noisy oscillations throughout the time
evolution, for all timestep sizes, never approaching zero.

Figure 3: The evolution of $\epsilon_{\Lambda}$ at $J_{y}=1.2$ using different
initial conditions (IC): a random number state (“IC-1”, brown line), an empty
state with all spins pointing ‘down’ (“IC-2”, red line), a full state with all
spins pointing ‘up’ (“IC-3”, orange line), and a state where each spin has
components $\expectationvalue{\sigma^{x}}=1$,
$\expectationvalue{\sigma^{y}}=1$, $\expectationvalue{\sigma^{z}}=-1$ (“IC-4”,
green line), for timesteps (a) $\delta t=10^{-1}$, (b) $\delta t=10^{-2}$, (c)
$\delta t=10^{-3}$. All other parameters are the same as in Fig. 1. Energies
given in units of $\kappa=1$.

We next explore if using different initial conditions affects whether a steady
state is found. Figures 3(a-c) show that the evolution of $\epsilon_{\Lambda}$
at $J_{y}=1.2$ remains noisy for various initial conditions: a random number
state (brown line), a state with all spins pointing ‘down’ (red line), a state
with all spins pointing ‘up’ (orange line), and a state where each spin has
components $\expectationvalue{\sigma^{x}}=1$,
$\expectationvalue{\sigma^{y}}=1$, $\expectationvalue{\sigma^{z}}=-1$ (green
line). Other initial conditions that we have tested (not shown here) produced
a similar behaviour as in Fig. 3.

Figure 4: The evolution of $\epsilon_{\Lambda}$ at selected values of $J_{y}$
during the adiabatic parameter sweep from $J_{y}=1.4$ (where a steady state
exists) to $J_{y}=1.2$ in steps of $\Delta J_{y}=0.01$ and using timesteps (a)
$\delta t=10^{-2}$, (b) $\delta t=10^{-3}$. All other parameters are the same
as in Fig. 1. Energies given in units of $\kappa=1$.

Another possible way to choose initial conditions for the problematic
parameter regime is an adiabatic parameter sweep. We first calculate the NESS
for a value of $J_{y}$ where iPEPO does reach a steady state, then change
$J_{y}$ in small steps, using the steady state at each value of $J_{y}$ as the
initial state for the next value. This strategy can bypass highly entangled
intermediate states where a very high $D$ may be needed. In our case, we start
at $J_{y}=1.4$, and gradually reduce $J_{y}$ in steps of $\Delta J_{y}=0.01$
to $J_{y}=1.2$. Figures 4(a,b) show the evolution of $\epsilon_{\Lambda}$ at
selected values of $J_{y}$ during the parameter sweep, with timesteps $\delta
t=10^{-2},10^{-3}$ respectively. We observe that for $J_{y}\geq 1.33$,
$\epsilon_{\Lambda}$ shows a decreasing trend, indicating that iPEPO finds a
steady state, while for $J_{y}\leq 1.32$ we again find noisy oscillations.
Smaller timesteps $\delta t$ and smaller sweeping steps $\Delta J_{y}$ (not
shown) lead to the same conclusion.

Figure 5: The evolution of $\epsilon_{\Lambda}$ at $J_{y}=1.2$ for selected
values of $\kappa$ during the adiabatic parameter sweep from a strong
dissipation regime with $\kappa=8$ to a weak dissipation regime with
$\kappa=1$ in steps of $\Delta\kappa=0.1$ and using timesteps (a) $\delta
t=10^{-2}$, (b) $\delta t=10^{-3}$. All other parameters are the same as in
Fig. 1. Energies given in units of $\kappa=1$.

A similar strategy is to start from a strong dissipation regime (i.e. large
$\kappa$) where we know the steady state is approximately factorizable, and
then perform an adiabatic sweep to lower $\kappa$. Figures 5(a,b) show the
time evolution of $\epsilon_{\Lambda}$ at $J_{y}=1.2$ for selected values of
$\kappa$, in a sweep starting from $\kappa=8$ reducing $\kappa$ in steps of
$\Delta\kappa=0.1$, with timesteps $\delta t=10^{-2},10^{-3}$ respectively.
Similarly to the $J_{y}$ sweep, we find a steady state exists for $\kappa\geq
5.2$, but beyond this we again find noisy oscillations. To summarise this
section, for those points where a steady state is not found, this result
appears to be robust to a variety of initial states and simulation protocols.

###  2.2 Effects of bond dimension

Figure 6: Evolution of $\epsilon_{\Lambda}$ with bond dimension at $J_{y}=1.2$
(a-c,g-i) and $J_{y}=1.5$ (d-f). Left column (panels a,d,g) are for timestep
$\delta t=10^{-1}$, middle (b,e,h) $\delta t=10^{-2}$, and right (c,f,i)
$\delta t=10^{-3}$. Panels (a-f) show $D=3,4,5,6$, while panels (g-i) show
$D=10,12,14,15$. All other parameters are the same as in Fig. 1. Energies
given in units of $\kappa=1$.

Figure 6 presents the effect of changing iPEPO bond dimensions D. Panels (a-c)
show the time evolution of $\epsilon_{\Lambda}$ at $J_{y}=1.2$ for timesteps
$\delta t=10^{-1},10^{-2},10^{-3}$ respectively. Each panel shows simulations
performed using different bond dimensions $3\leq D\leq 6$; no steady state is
found for any of these values of $D$. Panels (d-f) show the same quantities
but for $J_{y}=1.5$, where a steady state is known to occur at $D=4$. In this
case, notably, while iPEPO reaches a steady state for $D=3,4$, no steady state
is found for $D=5,6$. To explore this further, Panels (g-i) show the behavior
at $J_{y}=1.2$ for larger bond dimensions $10\leq D\leq 15$. We observe that
for $D=12$ a steady state is found for timesteps $\delta t=10^{-2},10^{-3}$.
However, increasing the bond dimension further to $D=14,15$ leads again to
noisy oscillations. These results suggest that while larger bond dimension may
eventually yield a meaningful steady state, spurious steady states can arise
at small bond dimension which then change as the bond dimension increases
further. In addition, we note that while we can run the SU time evolution for
$D=15$, the CTM calculations for this bond dimensions would require over $128$
GB of RAM, making it challenging to find observables without a support of high
performance distributed-memory calculations and quantum symmetries.

##  3 Conclusion

From the results of the last section, we may conclude that the SU iPEPO
algorithm at low bond dimensions is not always stable, reaching a steady state
only in some parameter regimes, typically away from dissipative critical
points. In other regimes, the algorithm failed to reach a steady state for all
bond dimensions $D$ that we could access. Moreover, in some cases, even when a
steady state is found for a given value of $D$, this may change as the bond
dimension is increased, switching instead to noisy time-dependent dynamics.
While we believe that there exists a value of $D$ allowing for a faithful
representation of the steady-state density matrix (when spatial correlations
decay exponentially), this study is unable to conclude what typical value of
the bond dimension is needed for a prototypical driven-dissipative lattice
model. Overall, the results shown here suggest that a significant caution is
required when extending the SU iPEPS algorithm in 2D to Liouvillian evolution.
Below we discuss a number of alternative approaches which may be employed
instead.

One alternative approach is to adapt the Full Update (FU) iPEPS algorithm [24,
3, 26, 37] to the Lindblad time evolution of mixed states. For closed systems,
the FU scheme [24, 38, 3] achieves an optimal truncation by using a
variational update scheme that computes the full environment at every step.
Since the Liouvillian evolution involves non-Hermitian operators, an issue
here is to find a reliable non-Hermitian algorithm that could substitute the
alternating least-squares scheme used in the two-site variational minimization
in the standard FU algorithm. There are, however, recent works on time
evolution in closed systems [39, 40, 41, 42] for which FU presents problems
with stability, meaning SU can be more accurate. A closely related idea is to
consider a global variational search algorithm that targets the null
eigenstate $\ket{\rho}$ of either Liouvillian $\mathcal{L}$ or a Hermitian
positive semidefinite object $\mathcal{L}^{\dagger}\mathcal{L}$; such
approaches were successfully used in one dimension [15, 16], and extension of
these methods to iPEPO has been discussed in Ref. [43]. Solving the
variational problem with $\mathcal{L}^{\dagger}\mathcal{L}$ is particularly
appealing since it allows reusing the standard and robust Hermitian
optimization algorithms. However, even when
$\mathcal{L}=\sum_{l}\mathcal{L}_{l,l+1}$ contains only nearest-neighbour
terms, the product
$\mathcal{L}^{\dagger}\mathcal{L}=\sum_{l,r}\mathcal{L}^{\dagger}_{l,l+1}\mathcal{L}_{l+r,l+r+1}$
will introduce highly nonlocal terms. While this is manageable in 1D [15], for
2D the nonlocal couplings may easily lead to unfeasibly large bond dimensions.
This may perhaps be adressed by truncating the range of these nonlocal terms
as has been discussed in 1D [44]. In addition, variational optimization iPEPO
approaches would require computationally expensive tensor contractions
involving both the iPEPS representing $\ket{\rho}$ and the iPEPO representing
either $\mathcal{L}$ or $\mathcal{L}^{\dagger}\mathcal{L}$.

A more promising approach, also discussed in Ref. [43] may be to extend novel
variational iPEPS techniques for ground state calculations in 2D introduced in
Refs. [45, 46], optimizing iPEPS tensors using tangent space methods or by
solving a local generalized eigenvalue problem. Notably, both approaches avoid
the need to construct a full PEPO for the Hamiltonian. Adapting these
algorithms to either $\mathcal{L}$ or $\mathcal{L}^{\dagger}\mathcal{L}$ could
dramatically reduce the computational costs that limit the practical use of
variational iPEPS methods. The global variational optimization could also
offer a potentially much more robust way of finding the NESS of $\mathcal{L}$
than one could hope to achieve with the standard iPEPS algorithm relying on
two-body updates.

## Acknowledgements

We acknowledge helpful discussions with Marzena Szymańska, Conor McKeever, and
Andrew Daley. We are grateful for comments from Augustine Kshetrimayum on an
earlier version of this manuscript.

##### Author contributions

The calculations presented here were performed by DK, using code developed by
DK with contributions from AB. The project was initially conceived by JK and
RF. All authors contributed to the writing of the manuscript.

##### Funding information

DK acknowledges support from the EPSRC Condensed Matter Centre for Doctoral
Training (EP/L015110/1). JK, RF and AB acknowledge the Kavli Institute for
Theoretical Physics, University of California, Santa Barbara (USA) for the
hospitality and support during the early stages of this work; this research
was supported in part by the National Science Foundation under Grant No. NSF
PHY11-25915.

##  Appendix A Implementation of iPEPO algorithm

As noted above, iPEPO is a simple extension of the iPEPS algorithm to
nonequilibrium steady states of Lindblad superoperators. As with other
extensions of tensor network methods [10] the main idea is to frame the
density matrix evolution through superoperators applied to vectorized many
body density matrices, which we denote $\ket{\rho}_{\sharp}$. In this appendix
we first briefly summarize the iPEPS algorithm and then discuss its extension
to open quantum systems. While these algorithms are standard and have been
described in full by Orús in [3], we summarise them here to provide a self-
contained description of the method that we apply to the open quantum systems.

###  A.1 Summary of the iPEPS algorithm

The basic idea behind PEPS is to parameterize the quantum state tensor
$\Psi_{k_{1},k_{2},\ldots,k_{N}}$ by a two-dimensional array of interconnected
rank-$5$ tensors (see Fig. 7). Each individual tensor represents a single site
of the quantum many-body system, with one vertical leg corresponding to the
local Hilbert space of dimension $d$, and four in-plane legs corresponding to
the bonds between different lattice sites. We denote the bond dimension of
PEPS by $D$, which limits the amount of entanglement that can be captured by
PEPS.

For translationally invariant systems, one may use the infinite PEPS (iPEPS)
ansatz [24] working directly in the thermodynamic limit. We can construct an
iPEPS by choosing a unit cell and representing its sites with tensors. We will
consider problems with a square unit cell. Since we use a Trotterised time
evolution that propagates pairs of sites, we will need only two on-site
tensors $A$ and $B$ to define iPEPS. We next describe the two main ingredients
of the iPEPS approach: the imaginary time propagation of iPEPS and the
calculation of the environment needed to extract observables.

Figure 7: The iPEPS time evolution using Eq. (4) involves propagating four
different bonds ‘U’, ‘D’, ‘R’, and ‘L’, indicated by different colours. The SU
algorithm uses Vidal form with $\Gamma^{[A,B]}$ site tensors and
$\Lambda^{[U,R,D,L]}$ diagonal bond matrices to represent iPEPS with a two-
site unit cell.

####  A.1.1 Time evolution: simple update

Time evolution can be performed by the Simple Update (SU) method, which
follows essentially the same main steps as the imaginary time infinite Time
Evolving Block Decimation (iTEBD) algorithm [47, 48, 49, 50]. In two
dimensions we perform Trotter decomposition by splitting our Hamiltonian into
four terms $H_{U}$, $H_{D}$, $H_{R}$ and $H_{L}$, describing respectively the
‘U’ (up), ‘D’ (down), ‘R’ (right), and ‘L’ (left) bonds of the lattice:

$H=H_{U}+H_{D}+H_{R}+H_{L}$ (3)

The first order Trotter decomposition of the time evolution operator
$U(\delta\tau)=e^{-H\delta\tau}$ then reads:

$U(\delta\tau)=e^{-\delta\tau H_{U}}e^{-\delta\tau H_{R}}e^{-\delta\tau
H_{D}}e^{-\delta\tau H_{L}}+O(\delta\tau^{2})$ (4)

where $\delta\tau$ is the imaginary timestep. Similarly to iTEBD in one
dimension, in SU we represent iPEPS using Vidal form: i.e. the iPEPS with a
two-site unit cell is fully specified by two $\Gamma^{[A,B]}$ site tensors and
four $\Lambda^{[U,R,D,L]}$ diagonal matrices that store the singular values of
iPEPS bonds, as seen in Fig. 7. We denote the local Hilbert space dimension by
$d$, and the bond dimension by $D$. The SU then consists of the following
steps:

  1. 1.

Absorb $\Lambda^{[R,D,L]}$ tensors on the external bonds into $\Gamma^{[A,B]}$
to obtain $Q^{[A,B]}$.

  2. 2.

Decompose each of $Q^{[A,B]}$ into subtensors $\mathrm{v}_{A,B}$ and $X_{A}$,
$Y_{B}$ using an exact SVD or QR/LQ decompositions. The original rank-5
tensors $\Gamma^{[A,B]}$ had the dimensionality of $dD^{4}$ giving rise to a
large computational cost $O(d^{3}D^{9})$ of the update procedure. However, an
update performed using the new rank-3 subtensors has a substantially reduced
cost of $O(d^{6}D^{3})$ since the dimensions of $\mathrm{v}_{A,B}$ are
considerably smaller and equal to $d\times dD\times D$ and $d\times D\times
dD$ respectively.

  3. 3.

Contract the two-body propagator $e^{-\delta\tau H_{U}}$ with
$\mathrm{v}_{A,B}$ and $\Lambda^{[U]}$ to form $\theta$ tensor.

  4. 4.

Decompose $\theta$ tensor into $\tilde{\mathrm{v}}_{A,B}$ and
$\tilde{\Lambda}^{[U]}$ tensors using SVD. To prevent the bond dimension of
our tensors from growing indefinitely, we must truncate
$\tilde{\mathrm{v}}_{A,B}$ and $\tilde{\Lambda}^{[U]}$ by retaining $D$
largest singular values and discarding the rest.

  5. 5.

Recover the updated rank-5 tensors $\tilde{Q}^{[A,B]}$ by contracting the
rank-3 subtensors $\tilde{\mathrm{v}}_{A,B}$ with $X_{A}$, $Y_{B}$
respectively.

  6. 6.

To restore $\Lambda^{[R,D,L]}$ on the external bonds, we divide each of
$Q^{[A,B]}$ by $\Lambda^{[R]}$, $\Lambda^{D}$, and $\Lambda^{[L]}$. This
procedure brings iPEPS back to its original Vidal form, with updated tensors
$\tilde{\Gamma}^{[A,B]}$ and $\tilde{\Lambda}^{[U]}$ for each ‘U’ bond on the
lattice.

All other steps are the same as in the one dimensional iTEBD algorithm. To
find the ground state we propagate iPEPS for $N$ imaginary timesteps
$\delta\tau$ until a steady state is achieved with respect to the spectrum of
singular values in $\Lambda^{[U,R,D,L]}$.

The SU algorithm is both simple and very efficient, with computational cost
$O(D^{3}d^{6})$ of the time evolution. However, SU is suboptimal since it
employs local truncations without taking into account the full environment of
a unit cell. For MPS, this issue can be resolved relatively easily by
transforming the tensor network into a canonical form, which orthonormalizes
other bonds surrounding the bond being truncated. This solution is not
possible in 2D, since there is no known canonical form for PEPS. To achieve an
optimal truncation, one must use a variational update scheme that computes the
full environment at every step. This procedure, known as the Full Update (FU)
[24, 3], is considerably more expensive and bears the computational cost of
$O(N\chi^{3}D^{6}+N\chi^{2}D^{8})$ where $N$ is the number of steps of the
imaginary time evolution. In practice, SU has been applied extensively to
various models and yields sufficiently accurate results for systems with large
gaps and sufficiently short correlation lengths [51]. Due to its simplicity
and efficiency, it allows shorter computation times and significantly higher
bond dimensions than FU, and thus remains popular. However, the suboptimal
truncation becomes an issue near quantum critical points when correlation
lengths become long, and in these cases FU should be used instead.

####  A.1.2 Contraction: corner transfer matrix

For the iPEPS representation to be of practical use, we must be able to
extract expectation values from it. Unlike MPS where we could evaluate
overlaps exactly at a polynomial cost, the exact contraction of two PEPS is an
exponentially hard problem that scales as $O(e^{L})$ with PEPS size $L$ [3].
Fortunately, there exist various computational algorithms that can perform
this contraction approximately with high precision. For infinite systems,
these methods typically proceed by computing the approximate environment of an
iPEPS unit cell. This effective environment consists of a small set of tensors
that represent the infinite tensor network surrounding the unit cell. Possibly
the most successful technique for computing iPEPS environments is the corner
transfer matrix (CTM) method [25, 52, 53], which will be the method we use. In
this section we will explain the details of the CTM algorithm.

Since observables for quantum states involve the overlap of two copies of the
state, the starting point for CTM is the contraction of two iPEPS, which
produces an infinite 2D network made of reduced tensors $a$. Each reduced
tensor $a$ results from the contraction of $\left[M^{A}\right]^{{\dagger}}$
and $M^{A}$ iPEPS tensors by their physical indices, except at the sites where
any operator is applied, leading to a different reduced tensor $a_{O}$.
Supposing the bond indices of $M^{A}$ had dimension $D$, the reduced tensors
$a$ now have bonds of dimension $D^{2}$.

For simplicity of exposition, we will start by considering a one-site unit
cell. However, methods based on Trotter decomposition into even and odd bonds
modify the translational invariance from one-site to two-site. Therefore, in
practice we always use the two-site version of the CTM algorithm. Let us now
subdivide this network into a $1\times 1$ unit cell made of $a$ tensors, and
its environment that contains the remaining infinite tensor network in which
the unit cell in embedded (see Fig. 8). The key idea is to represent the
environment by a set of four corner matrices $\\{C_{1,2,3,4}\\}$, and four
transfer tensors $\\{T_{1,2,3,4}\\}$ – these tensors are connected by new
virtual bond indices of size $\chi$. Similarly to the bond dimension of MPS
and PEPS, the environmental bond dimension $\chi$ is the parameter that
controls the accuracy of the CTM approximation of environment. The goal of CTM
algorithm is to obtain the environmental tensors by performing a series of
coarse graining moves:

Figure 8: (a) The main steps (3-5) of the left move of the CTM algorithm for a
one-site unit cell containing tensor $a$. (b) The renormalization step (5) is
done by inserting isometries $ZZ^{+}=I$ into the left edge of CTM network. (c)
The insertion step of the left move of the CTM algorithm for a two-site unit
cell containing tensors $a$, $b$. (d) The renormalization step for a two-site
unit cell is done by inserting two types of isometries, $ZZ^{+}=I$ and
$WW^{+}=I$.

  1. 1.

Initialize the CTM tensors, e.g. using a random-number initialization, or a
mean field environment with $\chi=1$.

  2. 2.

Perform four coarse graining moves in the left, right, up, and down
directions. The left move involves the following steps, illustrated
graphically in Fig. 8(a):

  3. 3.

Insertion: insert an extra column into the CTM network that contains the unit
cell tensor $a$, and the transfer tensors $T_{1,3}$.

  4. 4.

Absorption: absorb the new column into the left side of the CTM environment by
contracting their respective tensors. This increases the environmental bond
dimension by $\chi\to D^{2}\chi$: to prevent the bonds from growing
indefinitely we must implement an appropriate truncation scheme.

  5. 5.

Renormalization: truncate the environmental tensors by inserting appropriate
isometries $ZZ^{+}=I$ that reduce the bond dimensions $D^{2}\chi\to\chi$ by
projecting onto a relevant subspace, as shown in Fig. 8(b).

  6. 6.

Repeat steps (2-5) to let the CTM environment grow in all four directions
until it converges. Convergence is typically achieved when the eigenspectrum
of each corner matrix $\\{C_{1,2,3,4}\\}$ reaches the fixed point.

The CTM algorithm for a two-site unit cell, containing two $a$ and $b$
tensors, follows the same steps as the one-site algorithm outlined above. The
environment is now specified by a set of four corner matrices
$\\{C_{1,2,3,4}\\}$, and eight transfer tensors
$\\{T^{a}_{1,2,3,4},\;T^{b}_{1,2,3,4}\\}$. The main difference is that in the
’Insertion’ step we now insert two new columns instead of just one, as shown
in Fig. 8(c). There are now two ’Absorption’ and ’Renormalization’ steps in
the algorithm: the absorption of each column is followed by renormalization to
reduce the bond dimension $D^{2}\chi\to\chi$. The two-site algorithm also
needs two types of isometries $ZZ^{+}=I$ and $WW^{+}=I$ in the ’Renormalize’
step, to obtain renormalized transfer tensors $\tilde{T}^{a,b}_{4}$ and corner
tensors $\tilde{C}_{1,4}$ in Fig. 8(d).

Figure 9: (a,b) The steps to compute $W,W^{+}$ isometries for the bonds split
by the ’cut-1’. (c) To compute $Z,Z^{+}$ isometries for the bonds split by the
’cut-2’, we perform a translationally-invariant shift by inserting two extra
rows of tensors in the green box, and repeat the steps in (a).

Clearly, the crucial step of CTM algorithm is calculating the isometries.
Several different methods exist, for instance the ones described in Refs. [25,
54, 37], which we have implemented and tested in the process of developing our
iPEPS code. The prescription we have found to work best was the one in Ref.
[37], which achieved a smoother convergence and a more efficient
representation of the environment than its predecessors. Figure 9 illustrates
graphically the calculation of $Z,W$ isometries to be used in the
’Renormalize’ step of the left move. In Fig. 9(a,b) we compute $W,W^{+}$
isometries to be inserted into the bonds split by the ’cut-1’. In the first
stage in Fig. 9(a), we contract the lower and upper parts of the network,
producing the tensors $Q_{A}$ and $Q_{B}$ respectively. In the second stage,
also in Fig. 9(a), we decompose $Q_{A}$ and $Q_{B}$ using an exact SVD to
obtain the $R_{A,B}$ and $V^{\dagger}_{A,B}$ tensors. In the third stage in
Fig. 9(b), we form the product
$I=R_{A}\left[R^{-1}_{A}R^{-1}_{B}\right]R_{B}$, and decompose
$\left[R^{-1}_{A}R^{-1}_{B}\right]=U\Lambda V^{\dagger}$ using SVD, this time
truncating to the $\chi$ dominant singular values. The $R_{A,B}$ matrices and
the SVD matrices are then combined in the symmetrized fashion
$I\approx\left[R_{A}U\Lambda^{1/2}\right]\left[\Lambda^{1/2}V^{\dagger}R_{B}\right]=WW^{+}$
to construct the isometries $W=\left[R_{A}U\Lambda^{1/2}\right]$ and
$W^{+}=\left[\Lambda^{1/2}V^{\dagger}R_{B}\right]$. To obtain the $Z,Z^{+}$
isometries for the bonds split by the ’cut-2’, we perform a translationally-
invariant shift by inserting two extra rows of tensors in the green box, as
shown in Fig. 9(c). Contracting the lower and upper parts of the network then
gives tensors $Q_{A,B}$ for the cut-2. We can now compute $Z,Z^{+}$ repeating
exactly the same steps as for the $W,W^{+}$ before. Once all isometries
$Z,Z^{+}$ and $W,W^{+}$ are available, one can finally use them to carry out
the renormalization in Fig. 8(d).

The CTM algorithm described above has a computational complexity of
$O(\chi^{3}D^{6}\\\ +\chi^{2}D^{8})$. Once we have found a converged CTM
environment, it can be used to compute various observables.

###  A.2 Extension to iPEPO

To extend the above method to an open quantum system, as noted above, one may
first represent the density matrix $\rho$ as a Projected Entangled Pair
Operator (PEPO), and then reshape (vectorize) it into a PEPS
$\ket{\rho}_{\sharp}$ by combining both physical indices at each site. The
problem of computing NESS for an infinite 2D lattice thus becomes equivalent
to the problem of finding the ground states of Hamiltonians using the SU iPEPS
algorithm, by replacing imaginary time Hamiltonian propagation with the real
time Liouvillian propagation. To distinguish between the iPEPS representing
wavefunctions and vectorized density operators, we will refer to the density
matrix version as the iPEPO algorithm.

The main steps of the iPEPO algorithm are the same as in Sec. A.1, except for
two differences that we discuss next. The first difference is the propagator.
The imaginary time two-body propagators $U_{\alpha}(\delta\tau)=e^{-\delta\tau
H_{\alpha}}$ with Hamiltonian $H_{\alpha}$ for a bond $\alpha\in\\{U,R,D,L\\}$
are replaced by the real time two-body propagators $U_{\alpha}(\delta
t)=e^{-\delta t\mathcal{L}_{\alpha}}$, where $\mathcal{L}_{\alpha}$ is the
two-body Liouvillian for a bond $\alpha\in\left\\{U,R,D,L\right\\}$ and $t$ is
the real time. The second difference is that observables are calculated using
$\expectationvalue{O}=\Tr\left[O\rho\right]$, instead of
$\expectationvalue{O}=\expectationvalue{O}{\Psi}$. Similarly, the correct
normalization $\Tr\left[\rho\right]=1$ in contrast to
$\innerproduct{\Psi}{\Psi}=1$. As such, when extracting observables from
iPEPO, the reduced tensors that are contracted to find the environment come
from tracing out local indices instead of computing inner products. We may
then apply the CTM method from Sec. A.1.2 to observables.

We have implemented our iPEPO code in Fortran [36], including the CTM
algorithm, the functionality required for computing local observables and two-
point correlators, the SU procedure for both NESS and ground state
calculations, as well as the FU procedure for ground state calculations (for
benchmarking purposes only). In our calculations we gradually decrease the
time step $\delta t$ during the simulation to reduce the effects of Trotter
error while keeping the computational cost low. To determine when the
calculation has reached a steady state we require that the spectrum of
singular values contained in each diagonal bond matrix
$\Lambda\in\left\\{\Lambda^{[U,D,R,L]}\right\\}$ in Fig. 7 stops changing
within some accuracy $\epsilon$. More specifically, we take the largest
difference between singular values in diagonal matrices $\Lambda_{n}$ and
$\Lambda_{n-1}$, at timesteps $n$ and $n-1$ respectively, rescaled by the
largest singular value $|\Lambda_{n}|_{\text{max}}$ and by timestep size
$\delta t$:

$\epsilon_{\Lambda}=\frac{|\Lambda_{n}-\Lambda_{n-1}|_{\text{max}}}{\delta
t\,|\Lambda_{n}|_{\text{max}}}.$ (5)

For a steady state we expect $\epsilon_{\Lambda}$ to approach zero (or more
precisely, a value depending on machine precision), as the eigenvalue spectrum
should cease changing. To determine numerically when to stop, we define a
steady state as being reached when $\epsilon_{\Lambda}<\epsilon$ for each
$\Lambda\in\left\\{\Lambda^{[U,D,R,L]}\right\\}$.

###  A.3 Benchmarking with ground state calculations

Due to the similarity between the iPEPS and iPEPO algorithms, we may benchmark
our implementation against the known numerical results from ground state
calculations of two models: the transverse field Ising model, and the hardcore
Bose-Hubbard model.

##### Transverse-field Ising model

Our first test problem is the transverse field Ising model on an infinite
square lattice,

$H=-J\sum_{\left\langle
i,j\right\rangle}\sigma^{z}_{i}\sigma^{z}_{j}-g\sum_{i}\sigma^{x}_{i}.$ (6)

Here $\sigma^{x,z}_{i}$ are Pauli matrices at site $i$, $J$ is the nearest-
neighbour coupling between spins, and $g$ is the transverse magnetic field
along the $x$ axis. The ground state of this model exhibits a second order
phase transition between a paramagnetic phase at large $g$, and a
ferromagnetic phase at small $g$; the order parameter of this transition is
the longitudinal magnetization
$M_{z}=\expectationvalue{\sigma^{z}}{\Psi_{\text{GS}}}$. We show results in
units where $J=1$.

Figure 10(a,b) shows the longitudinal magnetization $M_{z}$ and transverse
magnetization $M_{x}$ as a function of transverse field $g$ in the vicinity of
phase transition, for different values of iPEPS bond dimension $D$. Our
implementation of iPEPS reproduces accurately both the SU and FU results
reported in Refs. [24, 25, 26]. As expected, the SU and FU calculations match
well far from the critical point where correlations are short-ranged, and the
results converge fast with increasing values of $D$. Near the critical point
FU becomes considerably more accurate than SU at a given bond dimension. That
is, due to the diverging correlation length, a much higher value of $D$ is
needed for SU to achieve the same level of accuracy as FU with $D=2,3$. As
seen from Fig. 10(a), the FU calculation with $D=3$ predicts the critical
point around $g=3.05$, in good agreement with previous iPEPS results in Refs.
[24, 25, 26] and Quantum Monte Carlo results in Ref. [55].

Figure 10: Ground state results benchmarking our iPEPS implementation. Panels
(a,b) show (a) the longitudinal magnetization $M_{z}$ and (b) transverse
magnetization $M_{x}$ as a function of transverse field $g$ (in units of $J$)
in the vicinity of phase transition of the transverse field Ising model. These
are computed for different iPEPS bond dimensions $D$ using both SU and FU as
indicated. Panel (c) shows number density of bosons $n$ per lattice site and
the condensate fraction $n_{0}$ for the hardcore Bose-Hubbard model as a
function of chemical potential $\mu$ (in units of $J$), computed for different
iPEPS bond dimensions $D$ using SU.

##### Hardcore Bose-Hubbard model

Our second test case is the Bose-Hubbard model (BHM) on an infinite square
lattice,

$H=-J\sum_{\left\langle
i,j\right\rangle}\left(a^{\dagger}_{i}a_{j}+\text{H.c.}\right)-\mu\sum_{i}a^{\dagger}_{i}a_{i},$
(7)

where the occupations are restricted to $0$ or $1$ bosons on each lattice
site. Here, $a^{\dagger}_{i}$, $a_{i}$ are hardcore bosonic creation and
annihilation operators at site $i$, satisfying the commutation relation
$\left[a^{\dagger}_{i},a_{j}\right]=(1-2a^{\dagger}_{i}a_{i})\delta_{ij}$. $J$
is the hopping rate between adjacent sites, and $\mu$ is the onsite chemical
potential. This model undergoes a second order phase transition between the
superfluid and the Mott insulator phases at the critical value of $\mu/J$. As
before, we present results in units where $J=1$. Figure 10 shows the number
density of bosons $n=\expectationvalue{a^{\dagger}a}{\Psi_{\text{GS}}}$ and
the condensate fraction $n_{0}=|\expectationvalue{a}{\Psi_{\text{GS}}}|^{2}$
as a function of chemical potential $\mu$, for different bond dimensions $D$
of SU iPEPS. Again, our iPEPS calculations reproduce accurately the ground
state results of Refs. [56, 26, 57].

## References

  * [1] F. Verstraete, V. Murg and J. I. Cirac,  _Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems_ ,  Adv. Phys. 57, 143 (2008),  10.1080/14789940801912366.
  * [2] U. Schollwöck,  _The density-matrix renormalization group in the age of matrix product states_ ,  Ann. Phys. 326, 96 (2011),  /10.1016/j.aop.2010.09.012.
  * [3] R. Orús,  _A practical introduction to tensor networks: Matrix product states and projected entangled pair states_ ,  Ann. Phys. 349, 117 (2014),  10.1016/j.aop.2014.06.013.
  * [4] I. Cirac, D. Perez-Garcia, N. Schuch and F. Verstraete,  _Matrix product states and projected entangled pair states: Concepts, symmetries, and theorems_ ,  URL http://arxiv.org/abs/2011.12127,  Preprint (2020).
  * [5] S. R. White,  _Density matrix formulation for quantum renormalization groups_ ,  Phys. Rev. Lett. 69, 2863 (1992),  10.1103/PhysRevLett.69.2863.
  * [6] S. R. White,  _Density-matrix algorithms for quantum renormalization groups_ ,  Phys. Rev. B 48, 10345 (1993),  /10.1103/PhysRevB.48.10345.
  * [7] A. E. Feiguin and S. R. White,  _Finite-temperature density matrix renormalization using an enlarged Hilbert space_ ,  Phys. Rev. B 72, 220401 (2005),  10.1103/PhysRevB.72.220401.
  * [8] S. R. White,  _Minimally entangled typical quantum states at finite temperature_ ,  Phys. Rev. Lett. 102, 190601 (2009),  10.1103/PhysRevLett.102.190601.
  * [9] A. J. Daley, C. Kollath, U. Schollwöck and G. Vidal,  _Time-dependent density-matrix renormalization-group using adaptive effective Hilbert spaces_ ,  J. Stat. Mech. Theory Exp. 2004, P04005 (2004),  10.1088/1742-5468/2004/04/P04005.
  * [10] M. Zwolak and G. Vidal,  _Mixed-state dynamics in one-dimensional quantum lattice systems: a time-dependent superoperator renormalization algorithm_ ,  Phys. Rev. Lett. 93, 207205 (2004),  10.1103/PhysRevLett.93.207205.
  * [11] M. J. Hartmann,  _Polariton crystallization in driven arrays of lossy nonlinear resonators_ ,  Phys. Rev. Lett. 104, 113601 (2010),  10.1103/PhysRevLett.104.113601.
  * [12] C. Joshi, F. Nissen and J. Keeling,  _Quantum correlations in the one-dimensional driven dissipative XY model_ ,  Phys. Rev. A 88, 063835 (2013),  10.1103/PhysRevA.88.063835.
  * [13] L. Bonnes, D. Charrier and A. M. Läuchli,  _Dynamical and steady-state properties of a bose-hubbard chain with bond dissipation: A study based on matrix product operators_ ,  Phys. Rev. A 90, 033612 (2014),  10.1103/PhysRevA.90.033612.
  * [14] M. Schiró, C. Joshi, M. Bordyuh, R. Fazio, J. Keeling and H. Türeci,  _Exotic attractors of the nonequilibrium rabi-hubbard model_ ,  Phys. Rev. Lett. 116, 143603 (2016),  10.1103/PhysRevLett.116.143603.
  * [15] J. Cui, J. I. Cirac and M. C. Bañuls,  _Variational matrix product operators for the steady state of dissipative quantum systems_ ,  Phys. Rev. Lett. 114, 220601 (2015),  10.1103/PhysRevLett.114.220601.
  * [16] E. Mascarenhas, H. Flayac and V. Savona,  _Matrix-product-operator approach to the nonequilibrium steady state of driven-dissipative quantum arrays_ ,  Phys. Rev.A 92, 022116 (2015),  10.1103/PhysRevA.92.022116.
  * [17] A. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert and S. Montangero,  _Positive tensor network approach for simulating open quantum many-body systems_ ,  Phys. Rev. Lett. 116, 237201 (2016),  10.1103/PhysRevLett.116.237201.
  * [18] J. Jin, A. Biella, O. Viyuela, L. Mazza, J. Keeling, R. Fazio and D. Rossini,  _Cluster mean-field approach to the steady-state phase diagram of dissipative spin systems_ ,  Phys. Rev. X 6, 031011 (2016),  10.1103/PhysRevX.6.031011.
  * [19] D. Kilda and J. Keeling,  _Fluorescence spectrum and thermalization in a driven coupled cavity array_ ,  Phys. Rev. Lett. 122, 043602 (2019),  10.1103/PhysRevLett.122.043602.
  * [20] E. Gillman, F. Carollo and I. Lesanovsky,  _Numerical simulation of critical dissipative non-equilibrium quantum systems with an absorbing state_ ,  New Journal of Physics 21(9), 093064 (2019),  10.1088/1367-2630/ab43b0.
  * [21] H. Landa, M. Schiró and G. Misguich,  _Multistability of driven-dissipative quantum spins_ ,  Phys. Rev. Lett. 124, 043601 (2020),  10.1103/PhysRevLett.124.043601.
  * [22] E. M. Stoudenmire and S. R. White,  _Studying two-dimensional systems with the density matrix renormalization group_ ,  Annu. Rev. Condens. Matter Phys. 3, 111 (2012),  10.1146/annurev-conmatphys-020911-125018.
  * [23] F. Verstraete, M. M. Wolf, D. Perez-Garcia and J. I. Cirac,  _Criticality, the area law, and the computational power of projected entangled pair states_ ,  Phys. Rev. Lett. 96, 220601 (2006).
  * [24] J. Jordan, R. Orús, G. Vidal, F. Verstraete and J. I. Cirac,  _Classical simulation of infinite-size quantum lattice systems in two spatial dimensions_ ,  Phys. Rev. Lett. 101, 250602 (2008),  10.1103/PhysRevLett.101.250602.
  * [25] R. Orús and G. Vidal,  _Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction_ ,  Phys. Rev. B 80, 094403 (2009),  10.1103/PhysRevB.80.094403.
  * [26] J. Jordan,  _Studies of Infinite Two-Dimensional Quantum Lattice Systems with Projected Entangled Pair States_ ,  Ph.D. thesis, The University of Queensland, Australia (2011).
  * [27] A. Kshetrimayum, H. Weimer and R. Orús,  _A simple tensor network algorithm for two-dimensional steady states_ ,  Nat. Commun. 8, 1291 (2017),  10.1038/s41467-017-01511-6.
  * [28] S. Finazzi, A. Le Boité, F. Storme, A. Baksic and C. Ciuti,  _Corner-space renormalization method for driven-dissipative two-dimensional correlated systems_ ,  Phys. Rev. Lett. 115, 080604 (2015),  10.1103/PhysRevLett.115.080604.
  * [29] R. Rota, F. Minganti, C. Ciuti and V. Savona,  _Quantum critical regime in a quadratically driven nonlinear photonic lattice_ ,  Phys. Rev. Lett. 122, 110405 (2019),  10.1103/PhysRevLett.122.110405.
  * [30] G. Carleo and M. Troyer,  _Solving the quantum many-body problem with artificial neural networks_ ,  Science 355, 602 (2017),  10.1126/science.aag2302.
  * [31] N. Yoshioka and R. Hamazaki,  _Constructing neural stationary states for open quantum many-body systems_ ,  Phys. Rev. B 99, 214306 (2019),  10.1103/PhysRevB.99.214306.
  * [32] A. Nagy and V. Savona,  _Variational quantum monte carlo method with a neural-network ansatz for open quantum systems_ ,  Phys. Rev. Lett. 122, 250501 (2019),  10.1103/PhysRevLett.122.250501.
  * [33] M. J. Hartmann and G. Carleo,  _Neural-network approach to dissipative quantum many-body dynamics_ ,  Phys. Rev. Lett. 122, 250502 (2019),  10.1103/PhysRevLett.122.250502.
  * [34] F. Vicentini, A. Biella, N. Regnault and C. Ciuti,  _Variational neural-network ansatz for steady states in open quantum systems_ ,  Phys. Rev. Lett. 122, 250503 (2019),  10.1103/PhysRevLett.122.250503.
  * [35] F. Minganti, A. Biella, N. Bartolo and C. Ciuti,  _Spectral theory of Liouvillians for dissipative phase transitions_ ,  Phys. Rev. A 98, 042118 (2018),  10.1103/PhysRevA.98.042118.
  * [36] D. Kilda, A. Biella, F. Nissen and J. Keeling,  _iPEPO_ ,  https://github.com/The-iPEPO-Project/iPEPO (2019).
  * [37] B. Bruognolo,  _Tensor network techniques for strongly correlated systems_ ,  Ph.D. thesis, Ludwig-Maximilians-Universität München (2017).
  * [38] M. Lubasch, J. I. Cirac and M.-C. Banuls,  _Algorithms for finite projected entangled pair states_ ,  Phys. Rev. B 90, 064425 (2014),  10.1103/PhysRevB.90.064425.
  * [39] C. Hubig and J. I. Cirac,  _Time-dependent study of disordered models with infinite projected entangled pair states_ ,  SciPost Physics 6(3) (2019),  10.21468/SciPostPhys.6.3.031.
  * [40] A. Kshetrimayum, M. Goihl and J. Eisert,  _Time evolution of many-body localized systems in two spatial dimensions_ ,  URL http://arxiv.org/abs/1910.11359,  Preprint (2019).
  * [41] C. Hubig, A. Bohrdt, M. Knap, F. Grusdt and J. I. Cirac,  _Evaluation of time-dependent correlators after a local quench in iPEPS: hole motion in the t-J model_ ,  SciPost Phys. 8, 21 (2020),  10.21468/SciPostPhys.8.2.021.
  * [42] A. Kshetrimayum, M. Goihl, D. M. Kennes and J. Eisert,  _Quantum time crystals with programmable disorder in higher dimensions_ ,  URL http://arxiv.org/abs/2004.07267,  Preprint (2020).
  * [43] H. Weimer, A. Kshetrimayum and R. Orús,  _Simulation methods for open quantum many-body systems_ ,  URL http://arxiv.org/abs/1907.07079,  In press (2020).
  * [44] A. A. Gangat, I. Te and Y.-J. Kao,  _Steady states of infinite-size dissipative quantum chains via imaginary time evolution_ ,  Phys. Rev. Lett. 119, 010501 (2017),  10.1103/PhysRevLett.119.010501.
  * [45] L. Vanderstraeten, J. Haegeman, P. Corboz and F. Verstraete,  _Gradient methods for variational optimization of projected entangled-pair states_ ,  Phys. Rev. B 94, 155123 (2016),  10.1103/PhysRevB.94.155123.
  * [46] P. Corboz,  _Variational optimization with infinite projected entangled-pair states_ ,  Phys. Rev. B 94, 035133 (2016),  10.1103/PhysRevB.94.035133.
  * [47] G. Vidal,  _Efficient classical simulation of slightly entangled quantum computations_ ,  Phys. Rev. Lett. 91, 147902 (2003),  10.1103/PhysRevLett.91.147902.
  * [48] G. Vidal,  _Efficient simulation of one-dimensional quantum many-body systems_ ,  Phys. Rev. Lett. 93, 040502 (2004),  10.1103/PhysRevLett.93.040502.
  * [49] G. Vidal,  _Classical simulation of infinite-size quantum lattice systems in one spatial dimension_ ,  Phys. Rev. Lett. 98, 070201 (2007),  10.1103/PhysRevLett.98.070201.
  * [50] R. Orus and G. Vidal,  _Infinite time-evolving block decimation algorithm beyond unitary evolution_ ,  Phys. Rev. B 78, 155117 (2008),  10.1103/PhysRevB.78.155117.
  * [51] P. Corboz, R. Orús, B. Bauer and G. Vidal,  _Simulation of strongly correlated fermions in two spatial dimensions with fermionic projected entangled-pair states_ ,  Phys. Rev. B 81, 165104 (2010),  10.1103/PhysRevB.81.165104.
  * [52] T. Nishino and K. Okunishi,  _Corner transfer matrix renormalization group method_ ,  J. Phys. Soc. Jpn. 65, 891 (1996),  10.1143/JPSJ.65.891.
  * [53] T. Nishino and K. Okunishi,  _Corner transfer matrix algorithm for classical renormalization group_ ,  J. Phys. Soc. Jpn. 66, 3040 (1997),  10.1143/JPSJ.66.3040.
  * [54] P. Corboz, J. Jordan and G. Vidal,  _Simulation of fermionic lattice models in two dimensions with projected entangled-pair states: Next-nearest neighbor hamiltonians_ ,  Phys. Rev. B 82, 245119 (2010),  10.1103/PhysRevB.82.245119.
  * [55] H. W. Blöte and Y. Deng,  _Cluster Monte Carlo simulation of the transverse Ising model_ ,  Phys. Rev. E 66, 066110 (2002),  10.1103/PhysRevE.66.066110.
  * [56] A. Kshetrimayum, M. Rizzi, J. Eisert and R. Orús,  _Tensor network annealing algorithm for two-dimensional thermal states_ ,  Phys. Rev. Lett. 122, 070502 (2019),  10.1103/PhysRevLett.122.070502.
  * [57] J. Jordan, R. Orús and G. Vidal,  _Numerical study of the hard-core Bose-Hubbard model on an infinite square lattice_ ,  Phys. Rev. B 79, 174515 (2009),  10.1103/PhysRevB.79.174515.