# Optimal pointwise sampling for $L^{2}$ approximation

Matthieu Dolbeault and Albert Cohen Sorbonne Université, UPMC Univ Paris 06,
CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, 4 place Jussieu, 75005 Paris,
France (matthieu.dolbeault@sorbonne-universite.fr)Sorbonne Université, UPMC
Univ Paris 06, CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, 4 place
Jussieu, 75005 Paris, France<EMAIL_ADDRESS>
###### Abstract

Given a function $u\in L^{2}=L^{2}(D,\mu)$, where $\mu$ is a measure on a set
$D$, and a linear subspace $V_{n}\subset L^{2}$ of dimension $n$, we show that
near-best approximation of $u$ in $V_{n}$ can be computed from a near-optimal
budget of $Cn$ pointwise evaluations of $u$, with $C>1$ a universal constant.
The sampling points are drawn according to some random distribution, the
approximation is computed by a weighted least-squares method, and the error is
assessed in expected $L^{2}$ norm. This result improves on the results in [6,
8] which require a sampling budget that is sub-optimal by a logarithmic
factor, thanks to a sparsification strategy introduced in [17, 18]. As a
consequence, we obtain for any compact class $\mathcal{K}\subset L^{2}$ that
the sampling number $\rho_{Cn}^{\rm rand}(\mathcal{K})_{L^{2}}$ in the
randomized setting is dominated by the Kolmogorov $n$-width
$d_{n}(\mathcal{K})_{L^{2}}$. While our result shows the existence of a
randomized sampling with such near-optimal properties, we discuss remaining
issues concerning its generation by a computationally efficient algorithm.

_MSC 2020:_ 41A65, 41A81, 93E24, 62E17, 94A20

##  1 Introduction

We study the approximation of a function $u\in L^{2}(D,\mu)$, where $\mu$ is a
measure on a set $D$, by an element $\tilde{u}$ of $V_{n}$, a subspace of
$L^{2}(D,\mu)$ of finite dimension $n$, based on pointwise data of of $u$.
Therefore, to construct $\tilde{u}$, we are allowed to evaluate $u$ on a
sample of $m$ points $X=\\{x^{1},\dots,x^{m}\\}\in D^{m}$. In addition, we
consider randomized sampling and reconstruction, in the sense that $X$ will be
drawn according to a distribution $\sigma$ over $D^{m}$, so the error
$u-\tilde{u}$ should be evaluated in some probabilistic sense. For the sake of
notational simplicity, having fixed $D$ and $\mu$, we write throughout the
paper

$L^{2}:=L^{2}(D,\mu)\quad{\rm
and}\quad\|v\|:=\|v\|_{L^{2}}=\left(\int_{D}|v|^{2}d\mu\right)^{1/2},$ (1)

as well as

$e_{n}(u):=\min_{v\in V_{n}}\|u-v\|.$ (2)

One typical applicative setting is the reconstrution of multivariate
functions, which corresponds to $D$ being a domain in $\mathbb{R}^{d}$.

Our main result is the following:

######  Theorem 1.

For some universal constants $C,K\geqslant 1$, and for any $n$-dimensional
space $V_{n}\subset L^{2}$, there exists a random sampling
$X=\\{x^{1},\dots,x^{m}\\}$ with $m\leqslant Cn$ and a reconstruction map
$R:D^{m}\times\mathbb{C}^{m}\mapsto V_{n}$, such that for any $u\in L^{2}$,

$\mathbb{E}_{X}\left(\|u-\tilde{u}\|^{2}\right)\leqslant Ke_{n}(u)^{2}$ (3)

where $\tilde{u}:=R(x^{1},\dots,x^{m},u(x^{1}),\dots,u(x^{m}))$.

The reconstruction map $R$ is obtained through a weighted least-squares method
introduced in [6], which has already been discussed in several papers, see [1,
5, 8, 9, 7, 15, 16]. The weights involved are given by the expression

$w:x\in D\mapsto n\min_{v\in
V_{n}}\frac{\|v\|^{2}}{|v(x)|^{2}}=\frac{n}{\sum_{j=1}^{n}|L_{j}(x)|^{2}},$
(4)

where the last formula holds for any $L^{2}$-orthonormal basis
$(L_{1},\dots,L_{n})$ of $V_{n}$, which, up to the factor $n$, is the
Christoffel function associated to the space $V_{n}$ and the space
$L^{2}(D,\mu)$. The weighted least-squares solution is then simply defined as

$\tilde{u}:=\underset{v\in
V_{n}}{\operatorname{argmin}}\sum_{i=1}^{m}w(x^{i})|u(x^{i})-v(x^{i})|^{2}.$
(5)

Introducing the discrete $\ell^{2}$ norm

$\|v\|_{X}^{2}:=\frac{1}{m}\sum_{i=1}^{m}w(x^{i})|v(x^{i})|^{2}$ (6)

and its associated scalar product $\langle\cdot,\cdot\rangle_{X}$, we get a
computable formula for $\tilde{u}$:

$\tilde{u}=\underset{v\in
V_{n}}{\operatorname{argmin}}\|u-v\|_{X}^{2}=P^{X}_{V_{n}}u,$ (7)

where $P^{X}_{V_{n}}$ denotes the orthogonal projection on $V_{n}$ with
respect to $\langle\cdot,\cdot\rangle_{X}$. Note that strictly speaking
$\|\cdot\|_{X}$ is not a norm over $L^{2}$, however the existence and
uniqueness of $P^{X}_{V_{n}}$ will be ensured by the second condition in Lemma
1 below, see Remark 2.

Therefore our main achievements lie in the particular choice of the random
sample $X$ for ensuring the near-optimal approximation and sampling budget in
Theorem 1.

Now, the proof of Theorem 1 relies on two conditions: first, the expectation
of $\|\cdot\|_{X}^{2}$ has to be bounded by $\|\cdot\|^{2}$ up to a constant.
Second, an inverse bound should hold almost surely, instead of just in
expectation, for functions $v$ in $V_{n}$. More precisely, one has:

######  Lemma 1.

Assume that $m$ and the law $\sigma$ of $X=\\{x^{1},\dots,x^{m}\\}$ are such
that

$\mathbb{E}(\|v\|_{X}^{2})\leqslant\alpha\|v\|^{2},\qquad v\in L^{2},$ (8)

and

$\|v\|^{2}\leqslant\beta\|v\|_{X}^{2}\;\;\text{ a.s.},\qquad v\in V_{n}.$ (9)

Then

$\mathbb{E}(\|u-\tilde{u}\|^{2})\leqslant(1+\alpha\beta)e_{n}(u)^{2}.$ (10)

###### Proof.

Denote $u_{n}$ the orthogonal projection of $u$ on $V_{n}$ with respect to the
$L^{2}(D,\mu)$ norm. Applying Pythagoras theorem both for $\|\cdot\|$ and
$\|\cdot\|_{X}$, one obtains

$\displaystyle\mathbb{E}(\|u-\tilde{u}\|^{2})$
$\displaystyle=\|u-u_{n}\|^{2}+\mathbb{E}(\|u_{n}-\tilde{u}\|^{2})$
$\displaystyle\leqslant\|u-u_{n}\|^{2}+\beta\mathbb{E}(\|u_{n}-\tilde{u}\|_{X}^{2})$
$\displaystyle=\|u-u_{n}\|^{2}+\beta\mathbb{E}(\|u_{n}-u\|_{X}^{2}-\|u-\tilde{u}\|_{X}^{2})$
$\displaystyle\leqslant\|u-u_{n}\|^{2}+\beta\mathbb{E}(\|u_{n}-u\|_{X}^{2})$
$\displaystyle\leqslant(1+\alpha\beta)\|u-u_{n}\|^{2},$

which proves (10) since $\|u-u_{n}\|=e_{n}(u)$. ∎

In section 2, we recall how both conditions (8) and (9) can be obtained with
$m$ quasi-linear in $n$, that is, of order $n\log n$. We reduce this budget to
$m$ of order $n$ in section 3, by randomly subsampling the set of evaluation
points, based on results from [17, 18]. The proof of Theorem 1 follows. We
compare it to the recent results [13, 14, 19] in section 4, in particular
regarding the domination of sampling numbers by $n$-widths. We conclude in
section 5 by a discussion on the offline computational cost for practically
generating the sample $X$.

##  2 Weighted least-squares

A first approach consists in drawing the $x^{i}$ independently according to
the same distribution $\rho$, that is, taking $\sigma=\rho^{\otimes m}$. The
natural choice for $\rho$ is $d\rho=\frac{1}{w}d\mu$, which is a probability
measure since

$\int_{D}\frac{1}{w}d\mu=\frac{1}{n}\sum_{j=1}^{n}\int_{D}|L_{j}(x)|^{2}d\mu(x)=\frac{1}{n}\sum_{j=1}^{n}\|L_{j}\|^{2}=1.$
(11)

We denote by $Z=\\{x^{1},\dots,x^{m}\\}$ this first random sample and by
$\|\cdot\|_{Z}$ the corresponding discrete $\ell^{2}$ norm. With this sampling
measure,

$\mathbb{E}(\|v\|_{Z}^{2})=\frac{1}{m}\sum_{i=1}^{m}\int_{D}w(x)|v(x)|^{2}d\rho=\int_{D}|v|^{2}d\mu=\|v\|^{2},$
(12)

so condition (8) is ensured for $X=Z$ with $\alpha=1$. To study the second
condition, we introduce the Hermitian positive semi-definite Gram matrix

$G_{Z}:=(\langle L_{j},L_{k}\rangle_{Z})_{j,k=1,\dots,n}$ (13)

and notice that (9) is equivalent to

$|\nu|^{2}=\Big{\|}\sum_{j=1}^{n}\nu_{j}L_{j}\Big{\|}^{2}\leqslant\beta\Big{\|}\sum_{j=1}^{n}\nu_{j}L_{j}\Big{\|}_{Z}^{2}=\beta\nu^{*}G_{Z}\nu,\qquad\nu\in\mathbb{C}^{n},$
(14)

which in turn rewrites as $\lambda_{\min}(G_{Z})\geqslant\beta^{-1}$.

By the central limit theorem, as $m$ tends to infinity, the scalar products
$\langle L_{j},L_{k}\rangle_{Z}$ converge almost surely to $\langle
L_{j},L_{k}\rangle=\delta_{j,k}$, so $G_{Z}$ converges to the identity matrix,
and we expect that $\lambda_{\min}(G_{Z})\geqslant\beta^{-1}$ holds for
$\beta>1$ with high probability as $m$ gets large. A quantitative formulation
can be obtained by studying the concentration of $G_{Z}$ around $I$ in the
matrix spectral norm

$\|M\|_{2}:=\max\\{|Mx|\,:\,|x|=1\\}.$

This is based on the matrix Chernoff bound, see [2, 22] for the original
inequality and [5], Lemma 2.1, for its application to our problem:

######  Lemma 2.

For $m\geqslant 10n\ln(\frac{2n}{\varepsilon})$, if $X\sim\rho^{\otimes m}$,
then

$\mathbb{P}\left(\|G_{Z}-I\|_{2}\leqslant\frac{1}{2}\right)\geqslant
1-\varepsilon.$ (15)

In particular,
$\mathbb{P}\left(\lambda_{\min}(G_{Z})\geqslant\frac{1}{2}\right)\geqslant
1-\varepsilon$.

Thus assumption (9) is satisfied with $\beta=2$, but only with probability
$1-\varepsilon$. As we would like it to hold almost surely, we condition the
sampling to the event

$E:=\left\\{\|G_{Z}-I\|_{2}\leqslant\frac{1}{2}\right\\}$ (16)

which defines a new sample

$Y=Z|E.$ (17)

In practice, $Y$ can be obtained through a rejection method, which consists in
drawing successively sets of points $Z^{1},Z^{2},\dots$ according to
$\rho^{\otimes m}$, and defining $Y=Z^{k}$ for the first value $k$ such that
$E$ holds. We then define $\tilde{u}$ as the weighted least-square estimator
based on this conditioned sample, that is

$\tilde{u}:=P_{V_{n}}^{Y}u.$ (18)

This approach was introduced and analyzed in [8], see in particular Theorem
3.6 therein. A simpler version of this result, sufficient for our purposes, is
the following:

######  Lemma 3.

For $m\geqslant 10n\ln(4n)$, if $Z\sim\rho^{\otimes m}$ and $Y=Z|E$, then

$\|G_{Y}-I\|_{2}\leqslant\frac{1}{2},$ (19)

and

$\mathbb{E}_{Y}(\|u-\tilde{u}\|^{2})\leqslant 5e_{n}(u)^{2}.$ (20)

###### Proof.

The first part immediately results from the definition of $Y$ and $E$, and
implies condition (9) with $\beta=2$. Moreover,
$\mathbb{P}(E)\geqslant\frac{1}{2}$ by Lemma 2 with $\varepsilon=\frac{1}{2}$,
so for any $v\in L^{2}(D,\mu)$,

$\mathbb{E}_{Y}(\|v\|_{Y}^{2})=\mathbb{E}_{Z}(\|v\|_{Z}^{2}|E)=\frac{\mathbb{E}_{Z}(\|v\|_{Z}^{2}\raise
1.29167pt\hbox{\large$\chi$}_{E})}{\mathbb{P}(E)}\leqslant\frac{\mathbb{E}_{Z}(\|v\|_{Z}^{2})}{\mathbb{P}(E)}\leqslant
2\|v\|^{2},$ (21)

so condition (8) holds with $\alpha=2$. The conclusion follows from Lemma 1. ∎

######  Remark 1.

The number of redraws $k$ for reaching $Y$ follows a geometric law of
expectation $\mathbb{E}(k)=\mathbb{P}(E)^{-1}=(1-\varepsilon)^{-1}$, that is
$\mathbb{E}(k)\leqslant 2$ for the particular choice of $m$ in the above
lemma. It should be well noted that $u$ is not evaluated at the intermediately
generated samples $Z^{1},\dots,Z^{k-1}$, which thus enter the offline cost of
the sampling algorithm.

######  Remark 2.

The fact that the Gramian $G_{Y}$ is non-singular implies that for any $u$
with given values $y^{i}$ at the points $x^{i}$, we can uniquely define

$\tilde{u}=P^{Y}_{V_{n}}u=\sum_{j=1}^{n}a_{j}L_{j},$ (22)

since $a=(a_{1},\dots,a_{n})^{*}$ solves the system of normal equations

$G_{Y}a=b,$ (23)

where the right-side vector has coordinates

$b_{j}=\langle
L_{j},u\rangle_{Y}=\frac{1}{m}\sum_{i=1}^{m}w(x^{i})L_{j}(x^{i})\overline{y^{i}}.$
(24)

If $u$ is in $L^{2}$, the $y^{i}$ are only defined up to a representer,
however since two representers $u^{1}$ and $u^{2}$ coincide $\mu$-almost
surely, we find that $P^{Y}_{V_{n}}u$ is well defined almost surely over the
draw of $Y$.

##  3 Random subsampling

With Lemma 3, we already have an error bound similar to that of Theorem 1.
However, the sampling budget is larger than $n$ by a logarithmic factor, which
we seek to remove in this section. To do so, we partition the sample $Y$ into
subsets of size comparable to $n$, and randomly pick one of these subsets to
define the new sample. An appropriate choice of the partitioning is needed to
circumvent the main obstacle, namely the preservation of condition (9). It
relies on the following lemma, taken from Corollary B of [18], itself a
consequence of Corollary 1.5 in [17]. The relevance of these two results to
sampling problems were exploited in [19] and noticed in [10], respectively.

######  Lemma 4.

Let $a_{1},\dots,a_{m}\in\mathbb{C}^{n}$ be vectors of norm
$|a_{i}|^{2}\leqslant\delta$ for $i=1,\dots,m$, satisfying

$\alpha I\leqslant\sum_{i=1}^{m}a_{i}a_{i}^{*}\leqslant\beta I$ (25)

for some constants $\delta<\alpha\leqslant\beta$. Then there exists a
partition of $\\{1,\dots,m\\}$ into two sets $S_{1}$ and $S_{2}$ such that

$\frac{1-5\sqrt{\delta/\alpha}}{2}\alpha I\leqslant\sum_{i\in
S_{j}}a_{i}a_{i}^{*}\leqslant\frac{1+5\sqrt{\delta/\alpha}}{2}\beta I,\qquad
j=1,2.$ (26)

In Lemma 2 of [18] this result is applied inductively in order to find a
smaller set $J\subset\\{1,\dots,m\\}$ of cardinality $|J|\leqslant cn$ such
that

$C^{-1}I\leqslant\frac{m}{n}\sum_{i\in J}a_{i}a_{i}^{*}\leqslant CI,$ (27)

for some universal constants $c,C>1$. We adapt this approach in order to
obtain a complete partition of $\\{1,\dots,m\\}$ by sets having such
properties.

######  Lemma 5.

Let $a_{1},\dots,a_{m}\in\mathbb{C}^{n}$ be vectors of norm
$|a_{i}|^{2}=\frac{n}{m}$ for $i=1,\dots,m$, satisfying

$\frac{1}{2}I\leqslant\sum_{i=1}^{m}a_{i}a_{i}^{*}\leqslant\frac{3}{2}I.$ (28)

Then there exists an integer $L$ and a partition of $\\{1,\dots,m\\}$ into
$2^{L}$ sets $J_{1},\dots,J_{2^{L}}$ such that

$c_{0}I\leqslant\frac{m}{n}\sum_{i\in J_{k}}a_{i}a_{i}^{*}\leqslant
C_{0}I,\qquad 1\leqslant k\leqslant 2^{L},$ (29)

with universal constants $c_{0}$ and $C_{0}$. In addition, each set $J_{k}$
satisfies

$|J_{k}|\leqslant C_{0}n.$ (30)

###### Proof.

The cardinality estimate (30) follows from the upper inequality in (29) by
taking the trace

$nC_{0}={\rm tr}(C_{0}I)\geqslant\frac{m}{n}\sum_{i\in J_{k}}{\rm
tr}(a_{i}a_{i}^{*})=\frac{m}{n}|J_{k}|\frac{n}{m}=|J_{k}|.$ (31)

For the proof of (29), if $n/m\geqslant 1/200$, then the result holds with
$L=0$, $J_{1}=\\{1,\dots,m\\}$, $c_{0}=1/2$ and $C_{0}=300$. Now assuming
$\delta:=n/m<1/200$, define by induction $\alpha_{0}=\frac{1}{2}$,
$\beta_{0}=\frac{3}{2}$, and

$\alpha_{\ell+1}:=\alpha_{\ell}\frac{1-5\sqrt{\delta/\alpha_{\ell}}}{2},\qquad\beta_{\ell+1}:=\beta_{\ell}\frac{1+5\sqrt{\delta/\alpha_{\ell}}}{2},\qquad\ell\geqslant
0.$ (32)

As $\alpha_{\ell+1}\leqslant\frac{\alpha_{\ell}}{2}$, the minimal integer $L$
such that $\alpha_{L}\leqslant 100\delta$ is well defined, and satisfies

$\alpha_{L}=\alpha_{L-1}\frac{1-5\sqrt{\delta/\alpha_{L-1}}}{2}>100\delta\frac{1-5\sqrt{1/100}}{2}=25\delta.$
(33)

Moreover $\alpha_{\ell}\geqslant 2^{L-\ell-1}\alpha_{L-1}\geqslant
2^{L-\ell-1}\,100\delta$ for $\ell=0,\dots,L-1$, so

$\beta_{L}=3\alpha_{L}\prod_{\ell=0}^{L-1}\frac{1+5\sqrt{\delta/\alpha_{\ell}}}{1-5\sqrt{\delta/\alpha_{\ell}}}\leqslant
C\delta,$ (34)

with $C:=300\prod_{\ell\geqslant
2}\frac{1+\sqrt{2}^{-\ell}}{1-\sqrt{2}^{-\ell}}$.

Finally, we inductively define partitions
$\\{S^{\ell}_{1},\dots,S^{\ell}_{2^{\ell}}\\}$ for $0\leqslant\ell\leqslant
L$: start with $S^{0}_{1}=\\{1,\dots,m\\}$ and for any $\ell,j$, noticing that

$\alpha_{\ell}I\leqslant\sum_{i\in
S^{\ell}_{j}}a_{i}a_{i}^{*}\leqslant\beta_{\ell}I,$ (35)

apply Lemma 4 to split $S^{\ell}_{j}$ into subsets $S^{\ell+1}_{2j-1}$ and
$S^{\ell+1}_{2j}$ satisfying the same property. At the last step, we define

$J_{k}=S^{L}_{k}.$ (36)

The framing (29) thus holds with $c_{0}=\alpha_{L}/\delta\geqslant 25$ and
$C_{0}=\beta_{L}/\delta\leqslant 11000$. ∎

Proof of Theorem 1: Define

$a_{i}=\left(\sqrt{\frac{w(x^{i})}{m}}L_{j}(x^{i})\right)_{j=1,\dots,n}$ (37)

the normalised random vectors corresponding to the sample
$Y=\\{x^{1},\dots,x^{m}\\}$ introduced in the previous section. As

$\frac{1}{2}I\leqslant
G_{Y}=\sum_{i=1}^{m}a_{i}a_{i}^{*}\leqslant\frac{3}{2}I$ (38)

and

$|a_{i}|^{2}=\frac{1}{m}w(x^{i})\sum_{j=1}^{n}|L_{j}(x^{i})|^{2}=\frac{n}{m}$
(39)

thanks to the choice of weights (4), the assumptions of Lemma 5 are satisfied.
Applying this lemma, we obtain sets $J_{1},\dots,J_{2^{L}}$ partitioning
$\\{1,\dots,m\\}$. Let $\kappa$ be a random variable taking value
$k\in\\{1,\dots,2^{L}\\}$ with probability $p_{k}=|J_{k}|/m$, and create a
random subsampling $X$ of $Y$ through

$X=\\{x^{i}\in Y\,:\,i\in J_{\kappa}\\}.$ (40)

Then the budget condition $|X|=|J_{\kappa}|\leqslant C_{0}n$ is satisfied
according to (30). Here, we define the discrete norm as

$\|v\|_{X}^{2}:=\frac{1}{|X|}\sum_{i\in J_{\kappa}}w(x^{i})|v(x^{i})|^{2},$
(41)

and the associated Gram matrix

$G_{X}:=(\langle
L_{j},L_{k}\rangle_{X})_{j,k=1,\dots,n}=\frac{m}{|J_{\kappa}|}\sum_{i\in
J_{\kappa}}a_{i}a_{i}^{*}.$ (42)

The weighted least-squares estimate is now defined as

$\tilde{u}:=\underset{v\in
V_{n}}{\operatorname{argmin}}\frac{1}{|X|}\sum_{i\in
J_{\kappa}}w(x^{i})|u(x^{i})-v(x^{i})|^{2},$ (43)

and it thus depends on the random draws of both $Y$ and $\kappa$. Condition
(9) follows from the lower inequality in (29) with $\beta=\frac{C_{0}}{c_{0}}$
since

$G_{X}\geqslant\frac{m}{|J_{\kappa}|}\frac{n}{m}c_{0}I\geqslant\frac{c_{0}}{C_{0}}I.$
(44)

Finally, we have for any $v\in L^{2}(D,\mu)$

$\mathbb{E}_{X}(\|v\|_{X}^{2})=\mathbb{E}_{Y}\left(\sum_{k=1}^{2^{L}}\frac{p_{k}}{|J_{k}|}\sum_{i\in
J_{k}}w(x^{i})|v(x^{i})|^{2}\right)=\mathbb{E}_{Y}(\|v\|_{Y}^{2})\leqslant
2e_{n}(u)^{2},$ (45)

so condition (8) holds with $\alpha=2$. Applying Lemma 1, we conclude that (3)
holds with $C=C_{0}$ and $K=1+2\frac{C_{0}}{c_{0}}$. $\Box$

##  4 Comparison with related results

In order to compare Theorem 1 with several recent results [10, 12, 19, 23], we
consider its implication when the target function $u$ belongs to a certain
class of functions $\mathcal{K}$ that describes some prior information on $u$,
such as smoothness.

Recall that if $V$ is a Banach space of functions defined on $D$ and
$\mathcal{K}\subset V$ is a compact set, its Kolmogorov $n$-width is defined
by

$d_{n}(\mathcal{K})_{V}:=\inf_{\dim V_{n}=n}\sup_{u\in\mathcal{K}}\inf_{v\in
V_{n}}\|u-v\|_{V},$ (46)

where the first infimum is taken over all linear spaces $V_{n}\subset V$ of
dimension $n$. This quantity thus describes the best approximation error that
can be achieved uniformly over the class $\mathcal{K}$ by an $n$-dimensional
linear space.

On the other hand, building a best approximation of $u$ requires in principle
full knowledge on $u$, and we want to consider the situation where we only
have access to a limited number of point evaluations. This leads one to
consider the sampling numbers, also called optimal recovery numbers, both in
the deterministic and randomized settings.

For deterministic samplings, we define the (linear) sampling numbers

$\rho_{m}^{\det}(\mathcal{K})_{V}:=\inf_{X,R_{X}}\max_{u\in\mathcal{K}}\|u-R_{X}(u(x^{1}),\dots,u(x^{m}))\|_{V},$
(47)

where the infimum is taken over all samples $X=\\{x^{1},\dots,x^{m}\\}\in
D^{m}$ and linear reconstruction maps $R_{X}:\mathbb{C}^{m}\to V$. For random
samplings, we may define similar quantities by

$\rho_{m}^{\rm
rand}(\mathcal{K})_{V}^{2}:=\inf_{X,R_{X}}\max_{u\in\mathcal{K}}\mathbb{E}\left(\|u-R_{X}(u(x^{1}),\dots,u(x^{m}))\|_{V}^{2}\right),$
(48)

where the infimum is taken over all random variables
$X=\\{x^{1},\dots,x^{m}\\}\in D^{m}$ and linear reconstruction maps
$R_{X}:\mathbb{C}^{m}\to V$. Note that a deterministic sample can be viewed as
a particular choice of random sample following a Dirac distribution in
$D^{m}$, and therefore

$\rho_{m}^{\rm rand}(\mathcal{K})_{V}\leqslant\rho_{m}^{\rm
det}(\mathcal{K})_{V}.$ (49)

Sampling numbers may also be defined without imposing the linearity of
$R_{X}$, leading to smaller quantities. In what follows, we shall establish
upper bounds on the linear sampling numbers, which in turn are upper bounds
for the nonlinear ones. We refer to [20] for an introduction and study of
sampling numbers in the context of general linear measurements, and to [21]
that focuses on point evaluation, also termed as standard information.

By optimizing the choice of the space $V_{n}$ used in Theorem 1, we obtain as
a consequence that, for $V=L^{2}=L^{2}(D,\mu)$, the sampling numbers in the
randomized setting are dominated by the Kolmogorov $n$-widths.

######  Corollary 1.

For any compact set $\mathcal{K}\subset L^{2}$, one has

$\rho_{Cn}^{\rm rand}(\mathcal{K})_{L^{2}}\leqslant
Kd_{n}(\mathcal{K})_{L^{2}},$ (50)

where $C$ and $K$ are the same constants as in Theorem 1.

######  Remark 3.

The bound (50) cannot be attained with independent and identically distributed
sampling points $x^{1},\dots,x^{m}$. Indeed, consider the simple example,
already evoked in [22], where $D=[0,1]$, $\mu$ is the Lebesgue measure,

$V_{n}=\left\\{\sum_{i=1}^{n}a_{i}\chi_{\left[\frac{i-1}{n},\frac{i}{n}\right[},(a_{1},\dots,a_{n})\in\mathbb{C}^{n}\right\\}$
(51)

is a space of piecewise constant functions, and $\mathcal{K}=\\{u\in
V_{n},\|u\|_{L^{\infty}}\leqslant 1\\}$. Then $\mathcal{K}\subset V_{n}$ so
$d_{n}(\mathcal{K})_{L^{2}}=0$, and an exact reconstruction $R_{X}u=u$ is
possible if and only if $X$ contains at least one point in each interval
$\left[\frac{i-1}{n},\frac{i}{n}\right[$. Thus $\rho_{n}^{\rm
det}(\mathcal{K})_{L^{2}}=0$, but in the case of i.i.d measurements, $m$ has
to grow like $n\log n$ to ensure this constraint, due to the coupon
collector’s problem.

######  Remark 4.

In [11], a result similar to Theorem 1 is obtained under the extra assumption
of a uniform bound on $e_{n}(u)/e_{2n}(u)$, yielding the validity of (50)
assuming a uniform bound on
$d_{n}(\mathcal{K})_{L^{2}}/d_{2n}(\mathcal{K})_{L^{2}}$. The recovery method
used in [11] is not of least-square type, but rather an elaboration of the
pseudo-spectral approach that would simply approximate the inner products
$\langle u,L_{j}\rangle=\int_{D}uL_{j}d\mu$ by a quadrature, using a
hierarchical approach introduced in [24].

Ideally, one would like a “worst case” or “uniform” version of Theorem 1, in
the form

$\rho_{Cn}^{\det}(\mathcal{K})_{L^{2}}\leqslant Kd_{n}(\mathcal{K})_{L^{2}},$
(52)

but it is easily seen that such an estimate cannot be expected for general
compact sets of $L^{2}$, due to the fact that pointwise evaluations are not
continuous in $L^{2}$ norm.

It is however possible to recover such uniform estimates by mitigating the
non-achievable estimate (52) in various ways. One first approach, developed in
[14, 23], gives an inequality similar to (52), with
$d_{n}(\mathcal{K})_{L^{2}}$ replaced by $d_{n}(\mathcal{K})_{L^{\infty}}$. It
is based on the following lemma, see Theorem 2.1 in [23], which we recall for
comparison with our Lemma 1:

######  Lemma 6.

Assume that $\mu$ is a finite measure of mass $\mu(D)=M<\infty$, that the
constant functions belong to $V_{n}$, and that there exists a sample
$X=\\{x^{1},\dots,x^{m}\\}$ and weights $w^{i}$ such that the discrete norm

$\|v\|_{X}^{2}=\frac{1}{|X|}\sum_{i=1}^{m}w_{i}|v(x^{i})|^{2}$ (53)

satisfies a framing

$\beta^{-1}\|v\|^{2}\leqslant\|v\|_{X}^{2}\leqslant\alpha\|v\|^{2},\quad v\in
V_{n}.$ (54)

Then

$\|u-P^{X}_{V_{n}}u\|\leqslant\sqrt{M}\left(1+\sqrt{\alpha\beta}\right)e_{n}(u)_{L^{\infty}},$
(55)

where $e_{n}(u)_{L^{\infty}}=\min_{v\in V_{n}}\|u-v\|_{L^{\infty}}$.

###### Proof.

For any $v\in L^{2}$, we have $\|v\|^{2}\leqslant M\|v\|_{L^{\infty}}^{2}$,
and as $1\in V_{n}$,

$\|v\|_{X}^{2}\leqslant\|1\|_{X}^{2}\|v\|_{L^{\infty}}^{2}\leqslant\alpha\|1\|^{2}\|v\|_{L^{\infty}}^{2}=\alpha
M\|v\|_{L^{\infty}}^{2}.$ (56)

Hence

$\displaystyle\|u-P^{X}_{V_{n}}u\|$
$\displaystyle\leqslant\|u-v\|+\|v-P^{X}_{V_{n}}u\|$
$\displaystyle\leqslant\|u-v\|+\sqrt{\beta}\|v-P^{X}_{V_{n}}u\|_{X}$
$\displaystyle\leqslant\|u-v\|+\sqrt{\beta}\|v-u\|_{X}$
$\displaystyle\leqslant(\sqrt{M}+\sqrt{\alpha\beta M})\|u-v\|_{L^{\infty}},$

and we conclude by optimizing over $v\in V_{n}$. ∎

Here, in contrast to the derivation of (10) in Lemma 1, one only uses the
framing property (54), and does not need the condition
$\mathbb{E}(\|v\|_{X}^{2})\leqslant\alpha\|v\|^{2}$. For this reason, one may
achieve the above objective with a simpler sparsification approach proposed in
[4] and adapted in [14], which performs a greedy selection of the points
$x^{i}$ within the sample $Y$, together with the definition of weights $w_{i}$
associated with these points. If the initial sample $Y$ satisfies

$\frac{1}{2}I\leqslant G_{Y}\leqslant\frac{3}{2}I,$ (57)

then, for any $c>1$ the selection algorithm produces a sample $X$ of at most
$cn$ points such that (54) holds with
$\alpha=\frac{3}{2}\left(1+\frac{1}{\sqrt{c}}\right)^{2}$ and
$\beta^{-1}=\frac{1}{2}\left(1-\frac{1}{\sqrt{c}}\right)^{2}$.

Optimizing the choice of $V_{n}$ (but imposing that constant functions are
contained in this space), this leads to the following comparison result
between deterministic optimal recovery numbers in $L^{2}$ and $n$-widths in
$L^{\infty}$: for any compact set $\mathcal{K}\in\mathcal{C}(D)$, one has

$\rho_{cn}^{\det}(\mathcal{K})_{L^{2}}\leqslant
C\sqrt{M}d_{n-1}(\mathcal{K})_{L^{\infty}},$ (58)

where $C$ depends on $c>1$. For $c=2$, one can take $C=11$. We refer to [14,
23] where this type of result is established.

Another approach consists in making pointwise evaluations continuous by
restriction to the case where $\mathcal{K}=B_{H}$ is the unit ball of a
reproducing kernel Hilbert space $H\subset L^{2}$, and assuming that the
sequence $(d_{n}(B_{H})_{L^{2}})_{n\geqslant 1}$ is $\ell^{2}$-summable. The
following result from [19], also based on the sparsification techniques from
[17], improves on a bound found in [12]

$\rho_{Cn}^{\det}(B_{H})^{2}_{L^{2}}\leqslant K\frac{\log
n}{n}\sum_{k\geqslant n}d_{k}(B_{H})_{L^{2}}^{2},$ (59)

More general compact classes $\mathcal{K}$ of $L^{2}$, such that point
evaluations are well defined on functions of $\mathcal{K}$, are considered in
[13], where the following general result is established: if

$d_{n}(\mathcal{K})_{L^{2}}^{2}\leqslant Cn^{-\alpha}\ln(n+1)^{\beta},\quad
n\geqslant 0,$ (60)

for some $\alpha>1$ and $\beta\in\mathbb{R}$, then

$\rho_{n}^{\det}(\mathcal{K})_{L^{2}}^{2}\leqslant\tilde{C}n^{-\alpha}\ln(n+1)^{\beta+1},\quad
n\geqslant 0.$ (61)

In the above results, the additional logarithmic factor appears as a residual
of the result obtained before sparsification, contrarily to the bounds (50)
and (58), which do not explicitely depend on the size of the initial sample
$Y$. This results in a gap of a factor $\log n$ between (59), (61) and known
lower bounds for $\rho_{n}^{\det}(\mathcal{K})_{L^{2}}^{2}$, see [19].

##  5 Computational aspects

The various results (50), (58), (59), (61) ensure the existence of good
sampling and reconstructions algorithms in various settings. We end by a
discussion on the computational cost of these strategies.

For the weighted least-squares methods corresponding to samples $Z$ and $Y$,
the most expensive step consists in assembling the matrix $G_{Z}$ as a sum of
$m$ matrices of size $n$, so the algorithmic complexity is of order
$\mathcal{O}(mn^{2})=\mathcal{O}(n^{3}\log n)$. Besides, to obtain $G_{Y}$,
this step may need to be repeated a few times, as explained in Remark 1, but
this only affects the offline complexity by a small random factor.

Note that we assumed that an orthogonal basis $(L_{1},\dots,L_{n})$ of $V_{n}$
is explicitly known, which might not be the case for irregular domains $D$.
However, under reasonable assumptions on $D$ or $V_{n}$, one can compute an
approximately orthogonal basis $(\widetilde{L}_{1},\dots,\widetilde{L}_{n})$,
either by performing a first discretization of $D$ with a large number of
points, or by using a hierarchical method on a sequence of nested spaces
$V_{1}\subset\dots\subset V_{n}$, see [1, 3, 5, 8, 15, 16]. These additional
steps have complexities $\mathcal{O}(K_{n}n^{2})$ and $\mathcal{O}(n^{4})$
respectively, where $K_{n}$ is the maximal value of the inverse Christoffel
function $\sum_{j=1}^{n}|L_{j}|^{2}$ which might grow more than linearly with
$n$ for certain choices of spaces $V_{n}$. Results similar to Lemma 3 have
been obtained in the above references, with $(L_{j})_{j=1,\dots,n}$ replaced
by $(\widetilde{L}_{j})_{j=1,\dots,n}$.

One could stop at this point and compute the approximation
$\tilde{u}=P^{Y}_{V_{n}}u$, which satisfies error bounds both in expectation
when comparing to $e_{n}(u)$, see Lemma 3, or uniformly when comparing to
$e_{n}(u)_{L^{\infty}}$, see Theorem 1 (iii) in [6]. Once the measurements of
$u$ are performed, the computation of $\tilde{u}$ requires to solve a $n\times
n$ linear system as in Remark 2, so the online stage takes a time
$\mathcal{O}(\tau n\log n+n^{3})$, where $\tau$ is the cost of each
measurement of $u$.

However, in applications where the evaluation cost $\tau$ becomes very high
(for example when each evaluation $x\mapsto u(x)$ requires solving a PDE by
some numerical code, or running a physical experiment), further reduction of
the size of the sample may prove interesting, and justifies the interest for
sparsification methods. The greedy selection method from [4], which is used in
[23] and leads to (58), has a complexity in
$\mathcal{O}(mn^{3})=\mathcal{O}(n^{4}\log n)$, but it can only be applied to
the worst-case setting, with the uniform error bound $e_{n}(u)_{L^{\infty}}$.

On the other hand, the iterative splitting method that we have used in this
paper following the ideas from [17, 19] is not easily implemented, and one
obvious method consists in testing all partitions of $\\{1,\dots,m\\}$ into
sets $S_{1}$ and $S_{2}$ when applying Lemma 4. Note that this lemma is in
practice used $L$ times, with $L=\mathcal{O}(\log\log n)$ since
$2^{L}=\mathcal{O}(\frac{m}{n})=\mathcal{O}(\log n)$. The algorithm consisting
in subdividing the sample $L$ times, each time checking that the Gram matrices
corresponding to $S_{1}$ and $S_{2}$ are well conditioned, and keeping one
such subset at random, thus has an exponential complexity
$\mathcal{O}(2^{m}n^{3})=\mathcal{O}(n^{cn})$. Having a different strategy
that would produce the random sample in polynomial time is currently an open
problem to us. Note that the hierarchical Monte-Carlo approaches from [24, 11]
have similar optimal error bounds with an optimal sampling budget, and without
exponential complexity in the generation of samples, however under the
additional assumption that is described in Remark 4.

We summarise these computational observations in the following table, which
illustrates the conflicts between reducing the sampling budget, ensuring
optimal approximation results, and maintaining a reasonable cost for sample
generation.

sampling complexity sample cardinality $m$ offline complexity
$\mathbb{E}(\|u-\tilde{u}\|^{2})$ $\leqslant Ce_{n}(u)^{2}$
$\|u-\tilde{u}\|^{2}$ $\leqslant Ce_{n}(u)_{L^{\infty}}^{2}$ conditioned
$\rho^{\otimes m}\,|\,E$ $10n\log(4n)$ $\mathcal{O}(n^{3}\log n)$ ✓ ✓ $+$
deterministic sparsification [4] $(1+\varepsilon)n$ $\mathcal{O}(n^{4}\log n)$
✗ ✓ $+$ random sparsification [17] $Cn$
$\mathcal{O}(n^{cn})\to\mathcal{O}(n^{r})$ ? ✓ ✓

As a final remark, let us to emphasize that although the results presented in
our paper are mainly theorical and not practically satisfactory, due both to
the computational complexity of the sparsification, and to the high values of
the numerical constants $C$ and $K$ in Theorem 1, they provide some intuitive
justification to the boosted least-squares methods presented in [8], which
consist in removing points from the initial sample as long as the
corresponding Gram matrix $G_{X}$ remains well conditioned. For instance,
Lemma 4 allows to keep splitting the sample even after $L$ steps, if one still
has a framing $\frac{1}{2}I\leqslant G_{X}\leqslant\frac{3}{2}I$ and a
sufficiently large ratio $\frac{|X|}{n}$. Nevertheless, it would be of much
interest to find a randomized version of [4] giving a bound of the form (50),
since this would give algorithmic tractability, smaller values for $C$ and
$K$, and the possibility to balance these constants in Theorem 1.

## References

  * [1] B. Adcock and J.M. Cardenas, Near-optimal sampling strategies for multivariate function approximation on general domains, SIAM Journal on Mathematics of Data Science, 2 (2020), pp. 607-630.
  * [2] R. Ahlswede and A. Winter, Strong converse for identification via quantum channels, IEEE Transactions on Information Theory, vol. 48, no 3 (2002), pp. 569-579.
  * [3] B. Arras, M. Bachmayr and A. Cohen, Sequential sampling for optimal weighted least squares approximations in hierarchical spaces, SIAM Journal on Mathematics of Data Science, 1 (2019), pp. 189-207.
  * [4] J. Batson, D.A. Spielman, and N. Srivastava, Twice-Ramanujan Sparsifiers, SIAM Review, 56 (2014), pp. 315-334.
  * [5] A. Cohen and M. Dolbeault, Optimal sampling and Christoffel functions on general domains, arXiv: 2010.11040 (2020), to appear in Constructive Approximation.
  * [6] A. Cohen and G. Migliorati, _Optimal weighted least squares methods,_ SMAI Journal of Computational Mathematics, 3 (2017), pp. 181-203.
  * [7] A. Doostan and J. Hampton, Coherence motivated sampling and convergence analysis of least squares polynomial Chaos regression, Computer Methods in Applied Mechanics and Engineering, 290 (2015), pp. 73-97.
  * [8] C. Haberstich, A. Nouy, and G. Perrin, Boosted optimal weighted least-squares, arXiv:1912.07075 (2019).
  * [9] J.D. Jakeman, A. Narayan, and T. Zhou, A Christoffel function weighted least squares algorithm for collocation approximations, Mathematics of Computation, 86 (2017), pp. 1913-1947.
  * [10] L. Kaemmerer, T. Ullrich, and T. Volkmer Worst case recovery guarantees for least squares approximation using random samples, arXiv: 1911.10111 (2019), to appear in Constructive Approximation.
  * [11] D. Krieg, Optimal Monte Carlo methods for $L^{2}$-approximation, Constructive Approximation, 49 (2019), pp. 385-403.
  * [12] D. Krieg and M. Ullrich, Function values are enough for $L^{2}$-approximation, arXiv: 1905.02516 (2020), to appear in Foundations of Computational Mathematics.
  * [13] D. Krieg and M. Ullrich, Function values are enough for $L^{2}$-approximation: Part II, Journal of Complexity (2021), p. 101569, https://doi.org/10.1016/j.jco.2021.101569.
  * [14] I. Limonova and V. Temlyakov, On sampling discretization in $L_{2}$, arXiv:2009.10789 (2020).
  * [15] G. Migliorati, Adaptive approximation by optimal weighted least-squares methods, SIAM Journal on Numerical Analysis, 7 (2019), pp. 2217-2245.
  * [16] G. Migliorati, Multivariate approximation of functions on irregular domains by weighted least-squares methods, IMA journal of numerical analysis (2020), https://doi.org/10.1093/imanum/draa023.
  * [17] A. Marcus, D. Spielman and N. Srivastava, Interlacing families II: Mixed characteristic polynomials and the Kadison-Singer problem, Annals of Mathemathics (2015), pp. 327-350, http://dx.doi.org/10.4007/annals.2015.182.1.8.
  * [18] S. Nitzan, A. Olevskii and A. Ulanovskii, Exponential frames on unbounded sets, Proceedings of the American Mathematical Society, vol. 144, no 1 (2016), pp. 109-118, http://dx.doi.org/10.1090/proc/12868.
  * [19] N. Nagel, M. Schäfer, and T. Ullrich, A new upper bound for sampling numbers, arXiv: 2010.00327 (2020), to appear in Foundations of Computational Mathematics.
  * [20] E. Novak and H. Woźniakowski, Tractability of multivariate problems. Volume I: Linear information, EMS Tracts in Mathematics, vol. 6, European Mathematical Society (EMS), Zürich (2008).
  * [21] E. Novak and H. Woźniakowski, Tractability of multivariate problems. Volume III: Standard information for operators, EMS Tracts in Mathematics, vol. 18, European Mathematical Society (EMS), Zürich (2012).
  * [22] J. Tropp, User-Friendly tail bounds for sums of random matrices, Foundations of Computational mathematics, 12 (2012), pp. 389-434.
  * [23] V. N. Temlyakov, On optimal recovery in $L^{2}$, Journal of Complexity, 65 (2020) p. 101545, https://doi.org/10.1016/j.jco.2020.101545.
  * [24] G.W. Wasilkowski and H. Woźniakowski, The power of standard information for multivariate approximation in the randomized setting, Mathematics of Computation, 76 (2007), pp. 965-988.


# Iterative construction of conserved quantities in dissipative nearly
integrable systems

Iris Ulčakar  Jožef Stefan Institute, 1000 Ljubljana, Slovenia  University of
Ljubljana, Faculty for physics and mathematics, 1000 Ljubljana, Slovenia  Zala
Lenarčič  Jožef Stefan Institute, 1000 Ljubljana, Slovenia  Jožef Stefan
Institute, 1000 Ljubljana, Slovenia University of Ljubljana, Faculty for
physics and mathematics, 1000 Ljubljana, Slovenia Jožef Stefan Institute, 1000
Ljubljana, Slovenia

###### Abstract

Integrable systems offer rare examples of solvable many-body problems in the
quantum world. Due to the fine-tuned structure, their realization in nature
and experiment is never completely accurate, therefore effects of
integrability are observed only transiently. One way to surpass that is to
couple nearly integrable systems to baths and driving: these will stabilize
integrable effects up to arbitrary time, as encoded in the time dependent, and
eventually, the stationary state of form of a generalized Gibbs ensemble.
However, the description of such driven dissipative nearly integrable models
is challenging and no exact analytical methods have been proposed so far. Here
we develop an iterative scheme in which integrability breaking perturbations
(baths) determine the most necessary conserved quantities to be added into a
truncated generalized Gibbs ensemble description. Our scheme significantly
reduces the complexity of the problem, paving the way for thermodynamic
results.

Integrable models have played a paramount importance in our understanding of
non-equilibrium dynamics since for them, in certain cases and regimes,
dynamics can be followed exactly. One of the modern milestone descriptions of
integrable models has been the observation that in the steady state reached
after a sudden excitation, such models are locally describable with a
generalized Gibbs ensemble (GGE) [1, 2, 3, 4, 5]. This came as a direct
generalization of the equilibration of thermalizing generic models, where the
steady state is locally describable with a thermal Gibbs ensemble [6, 7]. Both
conclusions underlined the importance of conservation laws for the long time
description, being energy for generic systems and all the extra conservation
laws for integrable systems. The applicability of generalized Gibbs ensembles
on a finite timescale was first established experimentally with cold atoms
[8]. The latest advances showed that also spatially inhomogeneous dynamics in
integrable systems can be formulated using a local GGE description and
generalized hydrodynamics [9, 10, 11], recently receiving experimental
confirmations [12, 13, 14, 15, 16].

Due to their fine-tuned interactions, integrable models cannot be exactly
realized in nature: while on short to intermediate timescales, integrable
models are realizable with quantum simulators [17, 18, 19] and can
approximately describe some materials [20, 21, 22], on longer timescales,
additional Hamiltonian terms, such as longer range interaction in trapped ions
experiment, coupling to phonons in solid state experiment, trapping potential,
transverse couplings and loss of atoms in cold atomic setups eventually
prevent integrable effects from sustaining up to arbitrary times [23, 13, 12,
24]. As suggested by one of us, the only way to sustain integrable effects in
nearly integrable systems up to arbitrary times is to drive them out of
equilibrium [25, 26, 27, 28]. In order to prevent heating in such a protocol,
we must allow for some openness so that excessive energy can be damped out.
The numerical evidence shows that then the time evolution [26] as well as the
steady state [25, 29, 27, 28] is again approximately given by the GGE. In this
case, the equations of motion for the Lagrange parameters associated with the
conserved quantities entering the GGE and their steady-state values are given
deterministically by the integrability breaking perturbations, i.e., by the
drive and coupling to the baths. This opens the possibility for GGE
engineering [27] and stabilizing potentially technologically useful phenomena
in spin chain materials approximately described by Heisenberg models, for
example, for efficient spin and energy pumping [25]. So far time dependent GGE
descriptions have also been used to describe the effect of particle loss in
cold atoms [30] and non-interacting systems coupled to baths [27, 31, 28, 32,
33, 34].

While the expanded role of GGE description to the case of weakly open nearly
integrable systems is fundamentally important, in this case concrete
calculations of dynamics and steady states are much harder than in isolated
systems. Weak integrability breaking, which causes non-elastic scatterings and
a slow reshuffling of quasiparticle content, makes the usual quasiparticle
treatment of integrable models much harder, for interacting systems
analytically (probably) impossible [30]. A numerically feasible approach is to
approximate the GGE with macroscopically many Lagrange parameters,
corresponding to macroscopically many conserved quantities, with its truncated
version. Such an approximation has been used in the context of isolated
integrable systems [35, 4], as well as in driven dissipative, nearly
integrable systems [25, 26, 27]. For the latter case, calculating the steady
state Langange multipliers involves solving a set of coupled $N_{C}$
equations, where $N_{C}$ is the number of retained conserved quantities. Here
we propose a new iterative method, which adds the most necessary conserved
quantities to the truncated GGE iteratively, as suggested by the driving and
dissipation itself. We examine the convergence to the exact result and show
that it is typically achieved within a few steps. This significantly reduces
the complexity of the problem since we do not have to solve macroscopically
many coupled equations but just a small number of them.

Setup. We consider driven dissipative, nearly integrable setups with dominant
unitary dynamics given by an integrable Hamiltonian $H_{0}$. Weak driving and
dissipation could be due to a Floquet unitary drive and coupling to a
(thermal) bath; however, for simplicity, we will consider weak coupling to
non-thermal Markovian baths, whose action is described by Lindblad operators
$L_{i}$. As pointed out in previous works [25, 29], different types of non-
static integrability-breaking perturbations should stabilize a GGE, and the
formalism described below can treat them all. The Liouville equation gives the
dynamics of the density matrix operator

$\hat{\mathcal{L}}\rho=-i[H_{0},\rho]+\hat{\mathcal{D}}\rho,\
\hat{\mathcal{D}}\rho=\epsilon\sum_{i}L_{i}\rho
L_{i}^{\dagger}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},\rho\\}$ (1)

where $\epsilon\ll 1$ is the absolute strength of the coupling to baths. We
will consider homogeneous coupling to baths, where Lindblad operators $L_{i}$
of the same form act on every site $i$.

Perturbatively, the zeroth order approximation to the steady state is of a
diagonal form in terms of eigenstates of $H_{0}$, $H_{0}\ket{m}=E_{m}\ket{m}$,

$\displaystyle\lim_{\epsilon\to
0}\lim_{t\to\infty}\rho=\sum_{m}a_{m}\ket{m}\bra{m}\equiv\rho_{D},$ (2)

with weights $a_{m}$ that are obtained from the kernel of the dissipator
projected on the diagonal subspace, $\sum_{n}D_{mn}a_{n}=0$, $D_{mn}=\langle
m|(\hat{\mathcal{D}}\ket{n}\bra{n})|m\rangle$, [29]. In the case of
degeneracies of $H_{0}$, Eq. (2) should be extended to a block-diagonal form;
however, we expect that thermodynamically both descriptions are equivalent and
we will work with diagonal ensemble here.

As suggested in our previous works [25, 26, 29], for integrable $H_{0}$, the
zeroth order diagonal ensemble $\rho_{D}$ should be thermodinamically
equivalent to a generalized Gibbs ensemble,

$\lim_{\epsilon\to
0}\lim_{t\to\infty}\lim_{L\to\infty}\rho=\frac{e^{-\sum_{m}\lambda_{m}C_{m}}}{\text{Tr}[e^{-\sum_{m}\lambda_{m}C_{m}}]}\equiv\rho_{\bm{\lambda}},$
(3)

with $C_{m}$ the (quasi-)local conserved quantities of the underlying
integrable model, $[C_{m},H_{0}]=0$, and $\lambda_{m}$ the associated Lagrange
multipliers. In the steady state, latter are determined by the stationarity
condition for all conserved quantities

$\langle\dot{C}_{m^{\prime}}\rangle=\text{Tr}\left[C_{m^{\prime}}\hat{\mathcal{D}}\frac{e^{-\sum_{m}\lambda_{m}C_{m}}}{\text{Tr}[e^{-\sum_{m}\lambda_{m}C_{m}}]}\right]\stackrel{{\scriptstyle!}}{{=}}0\quad\forall
m^{\prime}.$ (4)

This equation is very instructive: (i) it tells us that the form of
integrability breaking dissipator (the form of Lindblad operators) will
determine the $\lambda_{m}$ values, and (ii) in order to find $\lambda_{m}$,
one must solve a set of coupled non-linear equations, which is a hard problem.
In order to reduce the complexity of step (ii), an approximate description in
terms of a truncated GGE (tGGE) with a finite number of included conserved
quantities has been used [25, 26, 27], (mostly) resorting to the most local
ones. To some extent, this was sufficient to describe the expectation values
of local observables.

As shown in Ref. [26], also the dynamics towards the steady state can be
approximated with a time-dependent GGE, with the equation of motion for
$\lambda_{m}(t)$ derived by the use of super-projector $\hat{P}$ onto slow
modes, which are for the GGE Ansatz naturally given by the operators
$\frac{\partial\rho_{\bm{\lambda}}}{\partial\lambda_{m}}$ tangential to the
GGE manifold,

$\displaystyle\hat{P}X$
$\displaystyle=-\sum_{m,n}\frac{\partial\rho_{\bm{\lambda}}}{\partial\lambda_{m}}(\chi^{-1})_{m,n}\text{Tr}[C_{n}X],$
(5)

where
$\chi_{m,n}=-\text{Tr}[C_{m}\frac{\partial\rho_{\bm{\lambda}}}{\partial\lambda_{n}}]=\langle
C_{m}C_{n}\rangle_{\rho_{\bm{\lambda}}}-\langle
C_{m}\rangle_{\rho_{\bm{\lambda}}}\langle C_{n}\rangle_{\rho_{\bm{\lambda}}}$
is the $\\{m,n\\}$ entry of matrix $\chi$ and $\langle
O\rangle_{\rho_{\bm{\lambda}}}=\text{Tr}[\rho_{\bm{\lambda}}O]$. Applying the
super-projector to the slow dynamics on the GGE manifold,

$\displaystyle\hat{P}\dot{\rho}_{\bm{\lambda}}$
$\displaystyle=-\sum_{m,n}\frac{\partial\rho_{\bm{\lambda}}}{\partial\lambda_{m}}(\chi^{-1})_{m,n}\text{Tr}[C_{n}\hat{\mathcal{D}}\rho_{\bm{\lambda}}]=\sum_{m}\frac{\partial\rho_{\bm{\lambda}}}{\partial\lambda_{m}}\dot{\lambda}_{m}$
(6)

gives the rate of change for the Lagrange multiplier associated with $C_{m}$,

$\dot{\lambda}_{m}=-\sum_{n}(\chi^{-1})_{m,n}\text{Tr}[C_{n}\hat{\mathcal{D}}\rho_{\bm{\lambda}}],$
(7)

in the super-projector language given by the flow along the corresponding
tangential direction. Here, the initial conditions $\lambda_{m}(0)$ are given
by the initial state, as in the prethermal state [36, 37].

Iteratively constructed truncated GGE. In this Letter, we use the above super-
projector technique to iteratively construct the conserved quantities most
necessary for a truncated description of the steady state for a given
dissipator $\hat{\mathcal{D}}$. Given that the steady state Lagrange
parameters are selected by the dissipator, Eq. (4), we will, in the first
place, use $\hat{\mathcal{D}}$ to select the conserved quantities that we
include in the truncated GGE Ansatz. If such an iterative truncated
description converges to the exact one quickly, the procedure reduces the
number of conditions (4) that need to be solved.

In the procedure, we iteratively construct conserved quantities
$\tilde{C}_{k}$ from the basis $Q_{m}$, which should ideally be the set of
known (quasi)-local conserved quantities $Q_{m}=C_{m}$ of the integrable model
$H_{0}$. If not known or if hard to work with, one can also use projectors
$Q_{m}=\ket{m}\bra{m}$ onto eigenstates of $H_{0}$,
$H_{0}\ket{m}=E_{m}\ket{m}$. The iterative procedure has the following steps:

Step 0: Start with a thermal state
$\rho^{(0)}_{{\bm{\tilde{\lambda}}}}=e^{-\tilde{\lambda}^{(0)}_{0}H_{0}}/Z$
and find $\tilde{\lambda}^{(0)}_{0}$ from condition (4),
$\text{Tr}[H_{0}\hat{\mathcal{D}}e^{-\tilde{\lambda}^{(0)}_{0}H_{0}}]\stackrel{{\scriptstyle!}}{{=}}0$.

Step 1: Add the first iterative conserved quantity of the form

$\displaystyle\tilde{C}_{1}$
$\displaystyle=\mathcal{N}_{1}^{-1}\sum_{m}w^{(1)}_{m}Q_{m},$ $\displaystyle
w^{(1)}_{m}$
$\displaystyle=-\sum_{n}(\chi_{(0)}^{-1})_{mn}\text{Tr}[Q_{n}\hat{\mathcal{D}}\rho_{{\bm{\tilde{\lambda}}}}^{(0)}],$
(8)

where, according to Eq. (7), weights $w^{(1)}_{m}$ are given by the flows
along additional directions
$\frac{\partial\rho_{\bm{\lambda}}}{\partial\lambda_{m}}|_{\rho^{(0)}_{{\bm{\tilde{\lambda}}}}}$,
if we allow for a GGE manifold that is not one dimensional (thermal) as in
Step 0, but is spanned by additional $Q_{m}$ conserved quantities. A new
direction
$\frac{\partial\rho_{\bm{\lambda}}}{\partial\lambda_{m}}|_{\rho^{(0)}_{{\bm{\tilde{\lambda}}}}}$
that causes a stronger correction to the existing solution is supposedly more
important and should be weighted by a stronger bias $w_{m}^{(1)}$. At the end,
we find $\\{\tilde{\lambda}_{0}^{(1)},\tilde{\lambda}_{1}^{(1)}\\}$ for
$\rho^{(1)}_{\bm{\tilde{{\lambda}}}}\propto
e^{-\tilde{\lambda}_{0}^{(1)}H_{0}-\tilde{\lambda}_{1}^{(1)}\tilde{C}_{1}}$
from the condition (4) for $H_{0}$ and $\tilde{C}_{1}$.

Step k: Add $k$th iterative conserved quantity

$\displaystyle\tilde{C}_{k}$
$\displaystyle=\mathcal{N}_{k}^{-1}\sum_{m}w^{(k)}_{m}Q_{m},$ $\displaystyle
w^{(k)}_{m}$
$\displaystyle=-\sum_{n}\left(\chi_{(k-1)}^{-1}\right)_{mn}\text{Tr}[Q_{n}\hat{\mathcal{D}}\rho_{\bm{\tilde{\lambda}}}^{(k-1)}]$
(9)

and find $\\{\tilde{\lambda}_{k^{\prime}}^{(k)}\\}$ for
$\rho_{\bm{\tilde{\lambda}}}^{(k)}\propto
e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}^{(k)}\tilde{C}_{k^{\prime}}}$
from the set of $k+1$ conditions (4) for
$\\{\tilde{C}_{k^{\prime}}\\}_{k^{\prime}=0}^{k}$, where we denote the
$\tilde{C}_{0}=H_{0}$. Normalization $\mathcal{N}_{k}$ is in principle
arbitrary and can be absorbed into the corresponding Lagrange parameters.
However, it can also be chosen such that $\tilde{C}_{k}$ scales properly with
the system size. The susceptibility matrix, $(\chi_{(k)})_{m,n}=\langle
Q_{m}Q_{n}\rangle_{\rho_{\bm{\tilde{\lambda}}}^{(k)}}-\langle
Q_{m}\rangle_{\rho_{\bm{\tilde{\lambda}}}^{(k)}}\langle
Q_{n}\rangle_{\rho_{\bm{\tilde{\lambda}}}^{(k)}}$, must be evaluated in each
iterative step. Matrix $\chi_{(k)}$ is not invertible for the non-local basis
$Q_{m}=\ket{m}\bra{m}$, however, one can regularize it as explained in the
Supplementary material (SM). In case of (additional) unitary (Floquet)
perturbations, the iterative procedure can be generalized to cover those as
well, see SM.

In order to evaluate the convergence of our procedure, we define the distance
between density matrices via [35]

$d\left(\rho_{1},\rho_{2}\right)=\sqrt{\frac{\operatorname{Tr}[\left(\rho_{1}-\rho_{2}\right)^{2}]}{\operatorname{Tr}[\rho_{1}^{2}]+\operatorname{Tr}[\rho_{2}^{2}]}}$
(10)

and use it to compare reduced density matrices of iterative tGGE description
(9) and exact diagonal description (2). As noted in Ref. [35], the distance
between two GGEs on reduced space of $\ell$ consecutive sites scales as
$d(\rho_{\bm{\lambda},\ell},\rho_{\bm{\lambda}^{\prime},\ell})\sim\ell\sum_{i}|\lambda_{i}-\lambda^{\prime}_{i}|$
for large enough $\ell$.

Figure 1: Scaled distances $d(\rho_{\bm{\lambda},\ell},\rho_{D,\ell})/\ell$,
Eq. (10), between the reduced density matrices of support $\ell$ for the
diagonal steady state ensemble, Eq. (2), and different (truncated) generalized
Gibbs ensembles. Solid lines correspond to the GGEs including all
$N_{C}^{\text{all}}=2L-2$ local conserved quantities on a given system size.
From top to bottom we consider tGGEs with: (i) $N_{C}$ local conserved
quantities $C_{m},m=0,...,N_{C}-1$, (ii) iterative scheme after the $k$th step
for the basis $Q_{m}=C_{m}$, $m=1,...,N_{C}^{\text{all}}-1$, (iii) iterative
scheme after the $k$th step for basis $Q_{m}=\ket{m}\bra{m}$. Parameters:
Ising model with $h_{x}=0.6$, $J=1$, system sizes $L=10,12,14$.

Example I: Transverse field Ising model. First, we test our approach on the
paradigmatic non-interacting integrable model, the transverse field Ising
model $H_{0}$,

$H_{0}=\sum_{i}J\sigma^{z}_{i}\sigma^{z}_{i+1}+h_{x}\sigma^{x}_{i},$ (11)

which preserves a series of local extensive operators $C_{m}$ [38], explicitly
given in the SM. We choose the following Lindblad operators with $a=0.2$,

$L_{i}=S_{i}^{+}S_{i+1}^{-}+iS_{i+1}^{-}S_{i+2}^{+}+a\sigma_{i}^{x}\sigma_{i+1}^{z}.$
(12)

In Fig. 1, we show the scaled distances
$d(\rho_{\bm{\lambda},\ell},\rho_{D,\ell})/\ell$, Eq. (10), between the
reduced density matrices on $\ell$ consecutive sites for different (truncated)
GGE descriptions and the diagonal solution $\rho_{D}$, Eq. (2), for system
sizes $L=10,12,14$. Distances $d(\rho_{\bm{\lambda},\ell},\rho_{D,\ell})$
quantify the error of (truncated) GGE description on length-scale $\ell$. For
$\rho_{D}$ close to a GGE, we expect $d$ to scale with $\ell$ for large enough
$\ell$ [35], therefore normalization by $\ell$ is chosen. (i) In the first
row, we show the distances for the tGGE based on $N_{C}$ most local conserved
quantities $C_{m},m=0,\dots N_{C}-1$, where $C_{0}$ is the Hamiltonian. To
find the solution, $N_{C}$ coupled Eqs. (4) must be solved. As $N_{C}$ is
increased on a given system size, a better description with smaller distances
$d$ is obtained, up to the saturated values $d_{s,\ell}$ at
$N_{C}^{\text{all}}=2L-2$ (solid lines) corresponding to a GGE including all
local conserved quantities on a system of size $L$. As $L$ is increased,
saturated $d_{s,\ell}$ are decreased, signaling one can get increasingly
better results with a GGE Ansatz; in the thermodynamic limit and a GGE with
all local $C_{m}$ included, we expect that $d_{s,\ell}$ vanish. (ii) In the
middle row, we show the distances
$d(\rho^{(k)}_{\bm{\tilde{\lambda}},\ell},\rho_{D,\ell})/\ell$ for the
iterative scheme using the basis of all local conserved quantities
$Q_{m}=C_{m},m=1,...,N_{C}^{\text{all}}-1$ on the given system size, excluding
the Hamiltonian, which has its own Lagrange multiplier. The distances are
plotted as a function of the number of iterative steps $k$. The saturated
distances $d_{s,\ell}$ (solid lines), corresponding to GGE with all local
conserved quantities, are approximately obtained after only cca. $k=1$
iterative step. In the SM we give further details on the convergence of the
iterative procedure. The crucial advantage is that we now solve for $k+1=2$
instead of $N_{C}^{\text{all}}\sim 2L$ conditions (4) in order to find the
steady state. (iii) In the last row, we present the distance after the $k$th
step of the iterative scheme, obtained using the non-local basis
$Q_{m}=\ket{m}\bra{m}$ of projectors onto the eigenstates of $H_{0}$,
$H_{0}\ket{m}=E_{m}\ket{m}$. The saturated distances $d_{s,\ell}$, achievable
with involving all local $C_{m}$ (solid lines), are obtained after cca $k=1$
iterative step. Further reduction of $d$ is probably mostly due to non-local
contributions in the non-local basis $Q_{m}=\ket{m}\bra{m}$, which are
expected to vanish in the thermodynamic limit. With $L=10,12,14$ results, we
sense only the trend towards this.

Figure 2: Scaled distances $d(\rho_{\bm{\lambda},\ell},\rho_{D,\ell})/\ell$,
Eq. (10), between the reduced density matrices of support $\ell$ of the
diagonal steady state ensemble, Eq. (2), and different truncated generalized
Gibbs ensembles. Solid lines correspond to the tGGEs including
$N_{C}^{\max}=8$ most local conserved quantities. From top to bottom we
consider tGGEs with: (i) $N_{C}$ local conserved quantities
$C^{\text{loc}}_{m},m=0,...,N_{C}-1$, (ii) iterative scheme after the $k$th
step for basis $Q_{m}=C^{\text{loc}}_{m},m=1,...,N_{C}^{\max}-1$, (iii)
iterative scheme after the $k$th step for basis $Q_{m}=\ket{m}\bra{m}$.
Parameters: Heisenberg model with $J=1$, system sizes $L=12,14,16$.

Example II: Heisenberg model. Next, we test our approach on the paradigmatic
interacting integrable model, the isotropic Heisenberg model,

$\displaystyle
H_{0}=\sum_{i}\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1}+\sigma^{z}_{i}\sigma^{z}_{i+1},$
(13)

for Lindblad operators as in Eq. (12), with $a=0$, such that they preserve
magnetization. Therefore each magnetization sector has its own unique steady
state and we consider below the largest, zero magnetization sector. Unlike the
Ising model, Eq. (11), the Heisenberg model has a few families of conserved
operators; a family of local extensive operators [39] and two families of
quasi-local extensive operators [2, 40, 41].

Fig. 2 shows the same analysis previously done for the Ising model. (i) The
first row shows the distances for tGGE based on the first $N_{C}$ local
conserved quantities [39], derived using the boost operator,
$B=-i\sum_{j}jh_{j}$, $C^{\text{loc}}_{m+1}=[B,C^{\text{loc}}_{m}]$ and
$C^{\text{loc}}_{0}=H_{0}$, where $h_{j}$ are Hamiltonian densities
$H_{0}=\sum_{j}h_{j}$. In this notation, $C^{\text{loc}}_{1}$ is the conserved
energy current. For this model, approximate saturation in $d$ happens for a
rather small $N_{C}$ number of local conserved quantities, therefore we
perform our best calculation at $N_{C}^{\text{max}}=8$ local conserved
quantities $C_{m}^{\text{loc}}$ to obtain $d_{s,\ell}$ denoted by solid lines.
(ii) The second row shows the distances for the iterative approach, using the
basis of first $N_{C}^{\text{max}}=8$ local conserved operators
$Q_{m}=C^{\text{loc}}_{m},m=1,...,N_{C}^{\text{max}}-1$, excluding the
Hamiltonian. The asymptotic distances $d_{s,\ell}$ (solid line) are
approximately obtained after cca. $k=1$ iterative step. (iii) The third row
shows the distances for the iterative approach using the non-local basis
$Q_{m}=\ket{m}\bra{m}$ of projectors onto the eigenstates of $H_{0}$. Already
after the first iterative step, distances $d$ are smaller then for the tGGE
including local conserved quantities because the non-local basis $Q_{m}$
captures also the impact of quasi-local conserved quantities of the Heisenberg
model [2, 40, 41]. Partially, the reduction of $d$ in the non-local basis is
again due to non-local contributions on small systems considered.

Figure 3:  Thermodynamic norms
$\lim_{L\to\infty}\frac{1}{L}||\tilde{C}_{k}^{(s)}||^{2}$ of the part
$\tilde{C}_{k}^{(s)}$ of $k$th iterative conserved quantities
$\tilde{C}_{k}=\sum_{s}\tilde{C}_{k}^{(s)}$ that acts non-trivially only on
support $s$ for the example with isotropic Heisenberg Hamiltonian $H_{0}$ and
non-local basis $Q_{m}=\ket{m}\bra{m}$. Analysis indicates that iterative
conserved operators $\tilde{C}_{k}$ have a quasi-local nature.

In order to asses the nature of iterative conserved quantities constructed in
the non-local basis $Q_{m}=\ket{m}\bra{m}$, we plot in Fig. 3 the norm
$\lim_{L\to\infty}\frac{1}{L}||\tilde{C}_{k}^{(s)}||^{2}$ of the part
$\tilde{C}_{k}^{(s)}$ of conserved quantities
$\tilde{C}_{k}=\sum_{s}\tilde{C}_{k}^{(s)}$ that acts non-trivially only on
support $s$. The thermodynamic quantities are extracted from $1/L$ scaling of
the finite system results, see SM. From Fig. 3 we can conclude that the
thermodynamic iterative conserved quantities have a quasi-local nature, with
their operatorial weight on support $s$ dropping exponentially with $s$, even
though they are constructed using the non-local basis $Q_{m}=\ket{m}\bra{m}$.
Finite size effects related to non-local contributions are discussed in SM. We
use the Hilbert-Schmitt norm
$||A||^{2}=\frac{\text{Tr}[AA^{\dagger}]}{N_{H}}-\frac{\text{Tr}[A]\text{Tr}[A^{\dagger}]}{N_{H}^{2}}$,
where $N_{H}$ is the Hilbert space dimension. Here, the choice of
normalization $\mathcal{N}_{k}$, Eq. (9), is relevant and we use
$\mathcal{N}_{k}=\sqrt{\frac{1}{LN_{H}}\sum_{m}|w_{m}^{(k)}|^{2}}$ such that
$\tilde{C}_{k}$ are scaling extensively, $||\tilde{C}_{k}||^{2}\sim L$. A
posteriori, Fig. 3 confirms that our choice of normalization $\mathcal{N}_{k}$
is consistent with numerical results. One might expect that more essential
conserved quantities are more localized, however, we do not have evidence for
that from our small system size analysis. Access to larger systems $L$ and
supports $s$ would be needed to put this on firmer grounds.

Conclusions and outlook. In this Letter, we propose a new scheme for an
iterative construction of the most necessary conserved quantities of relevance
to describe nearly integrable, driven dissipative systems. Our approach is
motivated by the fact that in weakly driven and dissipative nearly integrable
systems, the integrability-breaking perturbations (couplings to baths and
drives) determine the steady state. Latter is for nearly integrable systems
approximately described by a generalized Gibbs ensemble with Lagrange
parameters fixed by the integrability breaking perturbations. Here, we use the
dissipator itself to select the combination of conserved quantities that
contribute significantly to the truncated GGE. Such a physically motivated
construction of truncated GGEs significantly reduced the complexity of
calculating steady state parameters, i.e., the Lagrange parameters. Namely,
instead of solving, in the thermodynamic limit, ideally, macroscopically many
coupled conditions (4) for all (quasi)local conserved quantities, we need to
solve for just a few since the iterative procedure usually yields sufficiently
good results within a few steps.

Our current study evaluated the iterative steps using exact diagonalization,
which limited our analysis to relatively small systems. However, our
approach’s central promise is in considering thermodynamically large systems,
as explained above. In order to realize this potential, one has to be able to
evaluate the weights in Eq. (9) in the thermodynamic limit. This should be
possible for at least some perturbations (Lindblad operators) and some
integrable models. When matrix $\chi$ is not diagonal, some simplifications
will be needed; however, if all ingredients (matrix elements of $\chi$ and
decay rates) are calculated in the thermodynamic limit, truncation to a finite
but large number of conserved quantities in the basis $Q_{m}$ should yield
satisfiable results. We leave the exploration of this aspects to a future
study.

###### Acknowledgements.

We thank A. Rosch and E. Ilievski for several useful discussions. We
acknowledge the support by the projects J1-2463, N1-0318 and P1-0044 program
of the Slovenian Research Agency, the QuantERA grant QuSiED by MVZI, QuantERA
II JTC 2021, and ERC StG 2022 project DrumS, Grant Agreement 101077265. ED
calculations were performed at the cluster ‘spinon’ of JSI, Ljubljana.

## References

  * Rigol  _et al._ [2007] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a Completely Integrable Many-Body Quantum System: An Ab Initio Study of the Dynamics of the Highly Excited States of 1D Lattice Hard-Core Bosons, Phys. Rev. Lett. 98, 050405 (2007).
  * Ilievski  _et al._ [2015] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler, and T. Prosen, Complete Generalized Gibbs Ensembles in an Interacting Theory, Phys. Rev. Lett. 115, 157201 (2015).
  * Vidmar and Rigol [2016] L. Vidmar and M. Rigol, Generalized Gibbs ensemble in integrable lattice models, J. Stat. Mech. Theory Exp. 2016, 064007 (2016).
  * Essler and Fagotti [2016] F. H. L. Essler and M. Fagotti, Quench dynamics and relaxation in isolated integrable quantum spin chains, J. Stat. Mech. Theory Exp. 2016, 064002 (2016).
  * Caux [2016] J.-S. Caux, The Quench Action, J. Stat. Mech. Theory Exp. 2016, 064006 (2016).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Langen  _et al._ [2015a] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, Experimental observation of a generalized Gibbs ensemble, Science 348, 207 (2015a).
  * Castro-Alvaredo  _et al._ [2016] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Emergent Hydrodynamics in Integrable Quantum Systems Out of Equilibrium, Phys. Rev. X 6, 041065 (2016).
  * Bertini  _et al._ [2016] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in Out-of-Equilibrium $XXZ$ Chains: Exact Profiles of Charges and Currents, Phys. Rev. Lett. 117, 207201 (2016).
  * Doyon [2020] B. Doyon, Lecture notes on Generalised Hydrodynamics, SciPost Phys. Lect. Notes , 18 (2020).
  * Malvania  _et al._ [2021] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, Generalized hydrodynamics in strongly interacting 1D Bose gases, Science 373, 1129 (2021).
  * Schemmer  _et al._ [2019] M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, Generalized Hydrodynamics on an Atom Chip, Phys. Rev. Lett. 122, 090601 (2019).
  * Cataldini  _et al._ [2022] F. Cataldini, F. Møller, M. Tajik, J. a. Sabino, S.-C. Ji, I. Mazets, T. Schweigler, B. Rauer, and J. Schmiedmayer, Emergent Pauli Blocking in a Weakly Interacting Bose Gas, Phys. Rev. X 12, 041032 (2022).
  * Møller  _et al._ [2021] F. Møller, C. Li, I. Mazets, H.-P. Stimming, T. Zhou, Z. Zhu, X. Chen, and J. Schmiedmayer, Extension of the Generalized Hydrodynamics to the Dimensional Crossover Regime, Phys. Rev. Lett. 126, 090602 (2021).
  * Yang  _et al._ [2023] K. Yang, Y. Zhang, K.-Y. Li, K.-Y. Lin, S. Gopalakrishnan, M. Rigol, and B. L. Lev, Phantom energy in the nonlinear response of a quantum many-body scar state, arXiv:2308.11615 (2023).
  * Kinoshita  _et al._ [2006] T. Kinoshita, T. Wenger, and D. S. Weiss, A quantum Newton’s cradle, Nature 440, 900 (2006).
  * Meinert  _et al._ [2015] F. Meinert, M. Panfil, M. J. Mark, K. Lauber, J.-S. Caux, and H.-C. Nägerl, Probing the Excitations of a Lieb-Liniger Gas from Weak to Strong Coupling, Phys. Rev. Lett. 115, 085301 (2015).
  * Langen  _et al._ [2015b] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, Experimental observation of a generalized Gibbs ensemble, Science 348, 207 (2015b).
  * Hess [2007] C. Hess, Heat conduction in low-dimensional quantum magnets, Eur. Phys. J. Spec. Top. 151, 73 (2007).
  * Niesen  _et al._ [2014] S. K. Niesen, O. Breunig, S. Salm, M. Seher, M. Valldor, P. Warzanowski, and T. Lorenz, Substitution effects on the temperature versus magnetic field phase diagrams of the quasi-one-dimensional effective Ising spin-$\frac{1}{2}$ chain system ${\mathrm{BaCo}}_{2}{\mathrm{V}}_{2}{\mathrm{O}}_{8}$, Phys. Rev. B 90, 104419 (2014).
  * Scheie  _et al._ [2021] A. Scheie, N. E. Sherman, M. Dupont, S. E. Nagler, M. B. Stone, G. E. Granroth, J. E. Moore, and D. A. Tennant, Detection of Kardar–Parisi–Zhang hydrodynamics in a quantum Heisenberg spin-1/2 chain, Nat. Phys. 17, 726 (2021).
  * Tang  _et al._ [2018] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Thermalization near Integrability in a Dipolar Quantum Newton’s Cradle, Phys. Rev. X 8, 021030 (2018).
  * Bastianello  _et al._ [2020] A. Bastianello, A. De Luca, B. Doyon, and J. De Nardis, Thermalization of a Trapped One-Dimensional Bose Gas via Diffusion, Phys. Rev. Lett. 125, 240604 (2020).
  * Lange  _et al._ [2017] F. Lange, Z. Lenarčič, and A. Rosch, Pumping approximately integrable systems, Nat. Commun. 8, 1 (2017).
  * Lange  _et al._ [2018] F. Lange, Z. Lenarčič, and A. Rosch, Time-dependent generalized Gibbs ensembles in open quantum systems, Phys. Rev. B 97, 165138 (2018).
  * Reiter  _et al._ [2021] F. Reiter, F. Lange, S. Jain, M. Grau, J. P. Home, and Z. Lenarčič, Engineering generalized Gibbs ensembles with trapped ions, Phys. Rev. Res. 3, 033142 (2021).
  * Schmitt and Lenarčič [2022] M. Schmitt and Z. Lenarčič, From observations to complexity of quantum states via unsupervised learning, Phys. Rev. B 106, L041110 (2022).
  * Lenarčič  _et al._ [2018] Z. Lenarčič, F. Lange, and A. Rosch, Perturbative approach to weakly driven many-particle systems in the presence of approximate conservation laws, Phys. Rev. B 97, 024302 (2018).
  * Bouchoule  _et al._ [2020] I. Bouchoule, B. Doyon, and J. Dubail, The effect of atom losses on the distribution of rapidities in the one-dimensional Bose gas, SciPost Phys. 9, 044 (2020).
  * Rossini  _et al._ [2021] D. Rossini, A. Ghermaoui, M. B. Aguilera, R. Vatré, R. Bouganne, J. Beugnon, F. Gerbier, and L. Mazza, Strong correlations in lossy one-dimensional quantum gases: From the quantum Zeno effect to the generalized Gibbs ensemble, Phys. Rev. A 103, L060201 (2021).
  * Gerbino  _et al._ [2023] F. Gerbino, I. Lesanovsky, and G. Perfetto, Large-scale universality in Quantum Reaction-Diffusion from Keldysh field theory, arXiv:2307.14945 (2023).
  * Perfetto  _et al._ [2023] G. Perfetto, F. Carollo, J. P. Garrahan, and I. Lesanovsky, Reaction-Limited Quantum Reaction-Diffusion Dynamics, Phys. Rev. Lett. 130, 210402 (2023).
  * Mi  _et al._ [2023] X. Mi, A. Michailidis, S. Shabani, K. Miao, P. Klimov, J. Lloyd, E. Rosenberg, R. Acharya, I. Aleiner, T. Andersen, _et al._ , Stable Quantum-Correlated Many Body States via Engineered Dissipation, arXiv:2304.13878 (2023).
  * Fagotti and Essler [2013] M. Fagotti and F. H. L. Essler, Reduced density matrix after a quantum quench, Phys. Rev. B 87, 245107 (2013).
  * Moeckel and Kehrein [2008] M. Moeckel and S. Kehrein, Interaction Quench in the Hubbard Model, Phys. Rev. Lett. 100, 175702 (2008).
  * Bertini  _et al._ [2015] B. Bertini, F. H. L. Essler, S. Groha, and N. J. Robinson, Prethermalization and Thermalization in Models with Weak Integrability Breaking, Phys. Rev. Lett. 115, 180601 (2015).
  * Grady [1982] M. Grady, Infinite set of conserved charges in the ising model, Phys. Rev. D 25, 1103 (1982).
  * Grabowski and Mathieu [1995] M. Grabowski and P. Mathieu, Structure of the Conservation Laws in Quantum Integrable Spin Chains with Short Range Interactions, Annals of Physics 243, 299 (1995).
  * Prosen and Ilievski [2013] T. Prosen and E. Ilievski, Families of Quasilocal Conservation Laws and Quantum Spin Transport, Phys. Rev. Lett. 111, 057203 (2013).
  * Ilievski  _et al._ [2016] E. Ilievski, M. Medenjak, T. Prosen, and L. Zadnik, Quasilocal charges in integrable lattice systems, J. Stat. Mech. Theory Exp. 2016, 064008 (2016).

a

Iterative construction of conserved quantities in dissipative nearly
integrable systems Iris Ulčakar Zala Lenarčič

Supplemental Material:
Iterative construction of conserved quantities in dissipative nearly
integrable systems
Iris Ulčakar1,2 and Zala Lenarčič2

1Faculty for Mathematics and Physics, University of Ljubljana, Jadranska ulica
19, 1000 Ljubljana, Slovenia

2Department of Theoretical Physics, J. Stefan Institute, SI-1000 Ljubljana,
Slovenia

In the Supplemental Material, we (i) show how to regularize the inverse of
susceptibility matrix $\chi^{-1}$ in the case of non-local basis of projectors
$Q_{m}=\ket{m}\bra{m}$, (ii) write the explicit form for the local conserved
quantities of the transverse field Ising model, (iii) we discuss the
convergence of the iterative approach in terms of Lagrange parameters, (iv)
give additional information about the nature of iteratively constructed
conserved quantities, and (v) outline the iterative procedure for unitary
perturbations.

##  S0.1 Regularization of $\chi^{-1}$ in the non-local basis

Here we address the invertibility of the susceptibility matrix $\chi$,
appearing in the calculation of weights $w_{m}^{(k)}$ at basis elements
$Q_{m}$ in the iterative conserved quantities $\tilde{C}_{k}$, Eq. (9),

$\displaystyle\tilde{C}_{k}$
$\displaystyle=\mathcal{N}_{k}^{-1}\sum_{m}w^{(k)}_{m}Q_{m},$ $\displaystyle
w^{(k)}_{m}$
$\displaystyle=-\sum_{n}(\chi_{(k-1)}^{-1})_{mn}\text{Tr}[Q_{n}\hat{\mathcal{D}}\rho_{\bm{\tilde{\lambda}}}^{(k-1)}]$
(S1)

for different choices of operator basis $Q_{m}$. Let us first rewrite the
equation for weights $w_{m}^{(k)}$, Eq. (S0.1), in a vector form. To avoid the
cluttering of notation, we omit the step index $k$ in the vector
representation. We denote the vector of weights with
$\ket{w}=\ket{w_{1}^{(k)},\dots w_{M}^{(k)}}$, with $n$th element being the
weight at $Q_{n}$. The vector index runs from $1$ up to the size of the
$Q_{m}$ basis $M$. Additionally we denote with
$\ket{q}=\ket{q_{1}^{(k-1)},\dots q_{M}^{(k-1)}}$ the vector with entries
$q_{n}^{(k-1)}=\text{Tr}[Q_{n}\hat{\mathcal{D}}\rho_{\bm{\tilde{\lambda}}}^{(k-1)}]$.
Matrix $\chi$ has size $M\times M$ and it’s elements are the connected
correlation functions between the basis elements $Q_{m}$,
$\chi_{mn}=\text{Tr}[Q_{m}Q_{n}\rho_{{\bm{\tilde{\lambda}}}}^{(k-1)}]-\text{Tr}[Q_{m}\rho_{{\bm{\tilde{\lambda}}}}^{(k-1)}]\,\text{Tr}[Q_{n}\rho_{{\bm{\tilde{\lambda}}}}^{(k-1)}]$.
The equation (S0.1) for weights is then written as

$\ket{w}=\chi^{-1}\ket{q}.$ (S2)

and can be solved if $\chi$ is invertible, which is the case for the basis of
local conserved quantities $Q_{m}=C_{m}$. In the remaining part of this
section, we show that for the non-local basis of projectors onto eigenstates
of $H_{0}$, $Q_{m}=\ket{m}\bra{m}$, matrix $\chi$ is non-invertible and we
demonstrate how to calculate weights $\ket{w}$ by regularizing $\chi^{-1}$.

For the non-local $Q_{m}=\ket{m}\bra{m}$, matrix elements of $\chi$ can be
expressed with coefficients of $\rho_{\bm{\tilde{\lambda}}}^{(k-1)}$ in the
diagonal ensemble $a_{m}=\bra{m}\rho_{\bm{\tilde{\lambda}}}^{(k-1)}\ket{m}$,
such that $\chi_{mn}=a_{m}(\delta_{mn}-a_{n})$. Due to the normalization of
the density matrix $\text{Tr}[\rho_{\bm{\tilde{\lambda}}}^{(k-1)}]=1$,
$N_{H}-1$ coefficients $a_{m}$ are free parameters, while one depends on the
others. Here, $N_{H}$ is the Hilbert space dimension and also the number of
basis elements $M=N_{H}$. This implies that matrix $\chi$ has $N_{H}-1$
linearly independent rows (or columns), and therefore rank $N_{H}-1$. One
vector is in the kernel of $\chi$, meaning $\chi$ is non-invertible. The
kernel vector $\ket{v}$ is in the eigenbasis of $H_{0}$ of form
$\ket{v}=\frac{1}{\sqrt{N_{H}}}\ket{1,\dots,1}$, since the sum of any row (as
well as column) in matrix $\chi$ in this basis equals zero,

$\chi\ket{1,\dots,1},\qquad\bra{1,\dots,1}\chi=0.$ (S3)

Note also that
$\text{Tr}[\hat{\mathcal{D}}\rho_{\bm{\tilde{\lambda}}}^{(k-1)}]=\sum_{n}q_{n}=0$,
or

$\braket{1,\dots,1}{q}=0$ (S4)

due to the dissipator $\hat{\mathcal{D}}$ in Eq. (1) being of Lindblad form.

To avoid the issue of invertibility of $\chi$, one can instead of Eq. (S2)
solve

$\chi\ket{w}=\ket{q},\quad\ket{w}=\ket{w^{\prime}}+c\ket{v}$ (S5)

with a set of solutions parametrized by $c\in\mathbb{C}$, since
$\chi(\ket{w^{\prime}}+c\ket{v})=\chi\ket{w^{\prime}}=\ket{q}$. Below we show
that

$\ket{w^{\prime}}=\tilde{\chi}^{-1}\ket{q},\qquad\tilde{\chi}=\chi+\ket{v}\bra{v}$
(S6)

is a solution giving rise to a traceless iterative conserved quantity
$\text{Tr}[\tilde{C}_{k}]=0$ and that the choice of factor $c$ does not affect
the iterative steady state density matrix $\rho^{(k)}_{\bm{\tilde{\lambda}}}$.

In Eq. (S6) we introduced a regularized inverse $\tilde{\chi}^{-1}$, motivated
by the fact that we can choose an orthonormal basis, where one basis vector is
the kernel vector $\ket{v}$ (subspace $\tilde{Q}$) while others form a
complementary subspace (subspace $\tilde{P}$) on which $\chi$ acts
nontrivially. Using conditions in Eqs. (S3, S4) and ordering the basis such
that the kernel vector is the last one, $\chi$, $\ket{v}$ and $\ket{q}$ have
the following form in the new basis:

$\displaystyle\chi_{\mathrm{b}}$
<EMAIL_ADDRESS>N_{H}-1}\end{math}}}\end{array}}&\scalebox{1.2}{\begin{math}{0}\end{math}}_{\scalebox{0.6}{\begin{math}{N_{H}-1\times
1}\end{math}}}\\\ \\\ \ \
\scalebox{1.2}{\begin{math}{0}\end{math}}_{\scalebox{0.6}{\begin{math}{1\times
N_{H}-1}\end{math}}}&0\end{array}\right)}_{\mathrm{b}}$ (S11)
$\displaystyle\ket{v}_{\mathrm{b}}$
$\displaystyle={\left(\begin{array}[]{c}\scalebox{1.3}{\begin{math}{0}\end{math}}_{\scalebox{0.6}{\begin{math}{N_{H}-1\times
1}\end{math}}}\\\ \\\
1\end{array}\right)}_{\mathrm{b}},\quad\ket{q}_{\mathrm{b}}={\left(\begin{array}[]{c}\scalebox{1.3}{\begin{math}{q_{\tilde{P}}}\end{math}}_{\scalebox{0.6}{\begin{math}{N_{H}-1\times
1}\end{math}}}\\\ \\\ 0\end{array}\right)}_{\mathrm{b}},$ (S18)

where $\chi_{\tilde{P}}$ and $q_{\tilde{P}}$ are the matrix $\chi$ and vector
$\ket{q}$ projected onto subspace $\tilde{P}$. Subscript $b$ labels the new
basis. $\tilde{\chi}$ and its inverse are then

<EMAIL_ADDRESS>\\\ \ \
<EMAIL_ADDRESS>\\\ \ \
\scalebox{1.2}{\begin{math}{0}\end{math}}&1\end{array}\right)}_{\mathrm{b}}.$
(S19)

Now one can easily see that the Ansatz in Eq. (S6) is the solution to Eq.
(S5):
$\chi_{\mathrm{b}}\tilde{\chi}_{\mathrm{b}}^{-1}\ket{q}_{\mathrm{b}}=\ket{q}_{\mathrm{b}}$.
Moreover, it is also straightforward to show that the iterative conserved
quantity, constructed from $\ket{w^{\prime}}$, is traceless

$\displaystyle\text{Tr}[\tilde{C}_{k}]$
$\displaystyle=\mathcal{N}_{k}^{-1}\sum_{m}w^{(k)}_{m}=\mathcal{N}_{k}^{-1}\braket{v}{w^{\prime}}=$
$\displaystyle=\mathcal{N}_{k}^{-1}\bra{v}_{\mathrm{b}}\tilde{\chi}_{\mathrm{b}}^{-1}\ket{q}_{\mathrm{b}}=0.$
(S20)

Let us note that the transformation to a new basis is useful for analytical
arguments, however, to work with $\tilde{\chi}$ in practice, the original
basis can be used with
$\ket{v}\bra{v}=\frac{1}{N_{H}}\ket{1,...,1}\bra{1,...,1}$, as in Eq. (S6),
expressed in that basis.

Finally we show that shifting the solution $\ket{w^{\prime}}$ by $c\ket{v}$,
Eq. (S5), does not affect the steady state density matrix
$\rho^{(k)}_{\bm{\tilde{\lambda}}}$ at the corresponding iterative step.
Adding $c\ket{v}$ in the language of the iterative tGGE means adding a
constant weight $c$ to each basis element $Q_{m}$, i.e.,

$\displaystyle\rho^{(k)}_{\bm{\tilde{\lambda}}}(c)$
$\displaystyle\equiv\frac{e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}\mathcal{N}_{k^{\prime}}^{-1}\sum_{m}(w_{m}^{(k^{\prime})}+c)Q_{m}}}{\text{Tr}[{e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}\mathcal{N}_{k^{\prime}}^{-1}\sum_{m}(w_{m}^{(k^{\prime})}+c)Q_{m}}}]}$
(S21)
$\displaystyle=\frac{\sum_{m}e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}\mathcal{N}_{k^{\prime}}^{-1}(w_{m}^{(k^{\prime})}+c)}Q_{m}}{\sum_{m}[{e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}\mathcal{N}_{k^{\prime}}^{-1}(w_{m}^{(k^{\prime})}+c)}}]}$
$\displaystyle=\frac{e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}\mathcal{N}_{k^{\prime}}^{-1}c}\sum_{m}e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}\mathcal{N}_{k^{\prime}}^{-1}w_{m}^{(k^{\prime})}}Q_{m}}{e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}\mathcal{N}_{k^{\prime}}^{-1}c}\sum_{m}[e^{-\sum_{k^{\prime}=0}^{k}\tilde{\lambda}_{k^{\prime}}\mathcal{N}_{k^{\prime}}^{-1}w_{m}^{(k^{\prime})}}]}$
$\displaystyle=\rho^{(k)}_{\bm{\tilde{\lambda}}}.$

We see that the contribution from $c\ket{v}$ cancels out, so we can skip it in
the first place. Similarly we could take a different shift $c\rightarrow
c_{k^{\prime}}$ for every iterative conserved quantity
$\tilde{C}_{k^{\prime}}$ and get the same result as above.

##  S0.2 Conservation laws of transverse Ising model

For transverse field Ising model,

$H_{0}=\sum_{i}J\sigma^{z}_{i}\sigma^{z}_{i+1}+h_{x}\sigma^{x}_{i},$ (S22)

the conservation laws have the following form

$\displaystyle C_{0}$ $\displaystyle=H_{0}$ $\displaystyle C_{2}$
$\displaystyle=\sum_{j}JS^{zz}_{{}_{j,j+2}}-h_{x}\sigma^{y}_{j}\sigma^{y}_{j+1}-h_{x}\hat{\sigma}^{z}_{j}\sigma^{z}_{j+1}-J\sigma^{x}_{j}$
$\displaystyle C_{2\ell>2}$
$\displaystyle=\sum_{j}JS^{zz}_{{}_{j,j+\ell+1}}-h_{x}S^{yy}_{{}_{j,j+\ell}}-h_{x}S^{zz}_{{}_{j,j+\ell}}+JS^{yy}_{{}_{j,j+\ell-1}}$
$\displaystyle C_{2\ell-1}$
$\displaystyle=J\sum_{j}S^{yz}_{{}_{j,j+\ell}}-S^{zy}_{{}_{j,j+\ell}},$ (S23)

with
$S^{\alpha\beta}_{i,j}=\sigma^{\alpha}_{i}\sigma^{x}_{i+1}\dots\sigma^{x}_{j-1}\sigma^{\beta}_{j}$.
In our analysis we use $N_{C}$ most local ones from $H_{0}$ to $C_{N_{C}-1}$.

##  S0.3 Convergence of iterative procedure

Figure S1:  The weights
$\Lambda_{m}^{(k)}=|\sum_{k^{\prime}=1}^{k}\tilde{\lambda}_{k^{\prime}}^{(k)}w_{m}^{(k^{\prime})}/\mathcal{N}_{k^{\prime}}|$
at each basis element $Q_{m}$ after $k=1,2,3$ iterative steps, for the Ising
example using the basis with all local conserved quantities
$Q_{m}=C_{m},m=1,...,2L-3$. The Hamiltonian has it’s own weight
$\Lambda_{0}^{(k)}=|\tilde{\lambda}_{0}^{(k)}|$. For comparison, the absolute
values of Lagrange multipliers $|\lambda_{m}|$ for the GGE based on all local
conserved quantities are shown (crosses). Parameters: transverse field Ising
model with $h_{x}=0.6$ and $J=1$ on $L=14$ sites.  Figure S2: (a-d) Finite
size scaling of norms $\frac{1}{L}||\tilde{C}_{k}^{(s)}||^{2}$ for the part
$\tilde{C}_{k}^{(s)}$ of $k=1,2,3,4$ iterative conserved operator
$\tilde{C}_{k}=\sum_{s}\tilde{C}_{k}^{(s)}$ that acts non-trivially only on
support $s$. These results correspond to the example with isotropic Heisenberg
Hamiltonian $H_{0}$ and non-local basis of projectors $Q_{m}=\ket{m}\bra{m}$
onto eigenstates of $H_{0}$. Thermodynamic weights at support $s$ for $k$th
iterative conserved quantity are extracted from $1/L$ scaling of system sizes
$L=12,14,16$. Smaller system sizes cannot be used to obtain results on larger
supports $s$.

Fast convergence of the procedure for the iterative construction of conserved
quantities can be illustrated by comparing the steady state weights
$\Lambda_{m}^{(k)}$ at different basis elements $Q_{m}$ in the iteratively
constructed tGGE to the corresponding Lagrange parameters $\lambda_{m}$ of the
best GGE Ansatz. In Fig. S1, we do that for the Ising example from the main
test, using the basis with all local conserved quantities $Q_{m}=C_{m}$,
$m=1,...,2L-3$, and the best GGE Ansatz solution obtained including all local
conserved quantities $C_{m}$, $m=0,...,2L-3$ on the given system size $L=14$.
We compare the weights
$\Lambda_{m}^{(k)}=|\sum_{k^{\prime}=1}^{k}\tilde{\lambda}_{k^{\prime}}^{(k)}w_{m}^{(k^{\prime})}/\mathcal{N}_{k^{\prime}}|$
at each basis element $Q_{m}$ after $k=1,2,3$ iterative steps. The Hamiltonian
has it’s own Lagrange multiplier $\tilde{\lambda}_{0}^{(k)}$ so that
$\Lambda_{0}^{(k)}=|\tilde{\lambda}_{0}^{(k)}|$. We find that the asymptotic
weights (absolute values of Lagrange parameters $|\lambda_{m}|$ denoted by
crosses) are reached very quickly, to a good approximation already after a
single iterative step ($k=1$). In the subsequent iterative steps, the
approximate solution of the first step is further refined.

Weights $\Lambda_{m}^{(k)}$ also assist to interpret which local conserved
quantities $C_{m}$ more crucially contribute to the steady state for the given
dissipator $\hat{\mathcal{D}}$. For the transverse field Ising model, Eq.
(11), and dissipator Eq. (12), the even conserved quantities $C_{2\ell}$, Eq.
(S0.2), seem to be more important and are well estimated after a single $k=1$
iterative step. In the next steps, the smaller weights at less important odd
conserved quantities $C_{2\ell-1}$ are captured as well. One should note that
the absolute value of $\Lambda_{m}^{(k)}$ depends on the normalization of
$C_{m}$; here we show Lagrange parameters in case of extensive normalization
$||\tilde{C}_{k}||^{2}=L$, for conserved operators following the definitions
in Eqs. (S0.2).

##  S0.4 Quasilocal nature of iterative conserved quantities

As we reported in the main text, the nature of the $k$th iterative conserved
quantity, constructed for the Heisenberg model in the non-local basis
$Q_{m}=\ket{m}\bra{m}$ of projectors onto eigenstate of $H_{0}$, is accessed
in the following way: we consider the norm
$\lim_{L\to\infty}\frac{1}{L}||\tilde{C}_{k}^{(s)}||^{2}$ of the part
$\tilde{C}_{k}^{(s)}$ of iterative conserved quantities
$\tilde{C}_{k}=\sum_{s}\tilde{C}_{k}^{(s)}$ that acts non-trivially only on
support $s$. We use the Hilbert-Schmitt norm
$||A||^{2}=\frac{\text{Tr}[AA^{\dagger}]}{N_{H}}-\frac{\text{Tr}[A]\text{Tr}[A^{\dagger}]}{N_{H}^{2}}$.
The thermodynamic quantities are extracted from $1/L$ scaling of the finite
system results. In Fig. S2, we plot the values of the norms
$\frac{1}{L}||\tilde{C}_{k}^{(s)}||^{2}$ for different supports $s$ on a given
system size $L$. For a quasi-local operator one expects an exponentially
decreasing behaviour with $s$, however, we see that on a given small system of
size $L$, the norms on support $s$ can even increase with $s$. We anticipate
that large overlaps with operators with large support $s$ are due to the non-
local part of the iterative conserved quantities $\tilde{C}_{k}$, if we
perform the procedure in the non-local basis $Q_{m}=\ket{m}\bra{m}$. As we
increase the system size, the non-local part reduces and when we extrapolate
the result to the thermodynamic limit via $1/L$ scaling shown in Fig. S2, for
large enough $s$ the norms are decaying exponentially with $s$, as shown in
Fig. 3 of the main text. This is a behaviour, expected from quasi-local
conserved operators [41]. Our iteratively conserved quantities are thus a
quasi-local superposition of local conserved quantities [39] and a linear
superposition of quasi-local conserved quantities of Heisenberg model [41].
From rather limited analysis on $s\leq 9$ and $L=12,14,16$, the localization
properties of all conserved quantities seem comparable. Access to larger
systems $L$ and supports $s$ would be needed to put this on firmer grounds.
One might naively expect that more essential quantities are more localized,
however, we do not have evidence supporting that at the moment.

##  S0.5 Iterative procedure for unitary perturbation

In case of (additional) unitary perturbation $H_{1}$, the dissipator
$\hat{\mathcal{D}}$ in Eq. (9) should be replaced (or supplemented, if both
types of perturbations occur simultaneously) by
$\hat{\mathcal{L}}_{1}\hat{Q}\hat{\mathcal{L}}_{0}^{-1}\hat{Q}\hat{\mathcal{L}}_{1}$,
i.e.,

$\displaystyle\tilde{C}_{k}$
$\displaystyle=\mathcal{N}_{k}^{-1}\sum_{m}w^{(k)}_{m}Q_{m},$ (S24)
$\displaystyle w^{(k)}_{m}$
$\displaystyle=-\sum_{n}\left(\chi_{(k-1)}^{-1}\right)_{mn}\text{Tr}[Q_{n}(\hat{\mathcal{D}}+\hat{\mathcal{L}}_{1}\hat{Q}\hat{\mathcal{L}}_{0}^{-1}\hat{Q}\hat{\mathcal{L}}_{1})\rho_{\bm{\tilde{\lambda}}}^{(k-1)}]$

where $\hat{\mathcal{L}}_{j}\rho=-i[H_{j},\rho]$ for $j=0,1$,
$\hat{\mathcal{L}}_{0}^{-1}\rho=-\int_{0}^{\infty}dte^{t\hat{\mathcal{L}}_{0}}\rho$
and $\hat{Q}\rho=\mathbb{1}-\sum_{m}\ket{m}\bra{m}\bra{m}\rho\ket{m}$ is a
super-projector orthogonal to the diagonal subspace, which ensures that the
action of $\hat{\mathcal{L}}_{0}^{-1}$ is not singular. Unlike for the
dissipator, where the contribution to the flow in the GGE manifold occurs
already in the first order, for the unitary perturbation the first nontrivial
contribution happens only in the second order as
$\text{Tr}[Q_{n}\hat{\mathcal{L}}_{1}\rho_{\bm{\tilde{\lambda}}}^{(k-1)}]=0$
due to cyclicity of the trace [25, 29].