# Dissipative Gradient Descent Ascent Method: A Control Theory Inspired
Algorithm for Min-max Optimization

Tianqi Zheng1, Nicolas Loizou1, Pengcheng You2 and Enrique Mallada1 1T. Zheng,
and E. Mallada are with the Department of Electrical and Computer Engineering,
N. Loizou is with the Department of Applied Mathematics and Statistics at
Johns Hopkins University, Baltimore, MD 21218, USA
<EMAIL_ADDRESS>You is with the Department of Industrial
Engineering and Management at Peking University, Beijing, China
<EMAIL_ADDRESS>
###### Abstract

Gradient Descent Ascent (GDA) methods for min-max optimization problems
typically produce oscillatory behavior that can lead to instability, e.g., in
bilinear settings. To address this problem, we introduce a dissipation term
into the GDA updates to dampen these oscillations. The proposed Dissipative
GDA (DGDA) method can be seen as performing standard GDA on a state-augmented
and regularized saddle function that does not strictly introduce additional
convexity/concavity. We theoretically show the linear convergence of DGDA in
the bilinear and strongly convex-strongly concave settings and assess its
performance by comparing DGDA with other methods such as GDA, Extra-Gradient
(EG), and Optimistic GDA. Our findings demonstrate that DGDA surpasses these
methods, achieving superior convergence rates. We support our claims with two
numerical examples that showcase DGDA’s effectiveness in solving saddle point
problems.

##  I Introduction

In recent years, there has been a significant focus on solving saddle point
problems, namely min-max optimization problems [1, 2, 3, 4, 5, 6]. These
problems have garnered considerable attention, particularly in fields such as
Generative Adversarial Networks (GANs) [7, 6, 8], Reinforcement Learning (RL)
[9], and Constrained RL (C-RL) [10, 11]. However, a major challenge that
persists in these approaches is the instability of the training process. That
is, solving the min-max optimization problem via running the standard Gradient
Descent Ascent (GDA) algorithm often leads to unstable oscillatory behavior
rather than convergence to the optimal solution. This is particularly
illustrated in bilinear min-max problems, such as the training of Wasserstein
GANs [12] or solving C-RL problems in the occupancy measure space [13], for
which the standard GDA fails to converge [1, 2].

In order to understand the instability of the GDA method and further tackle
its limitation, we draw inspiration from the control-theoretic notions of
dissipativity [14], which enables the design of stabilizing controllers using
dynamic (state-augmented) components that seek to dissipate the energy
generated by the unstable process. This aligns with recent work that leverages
control theory tools in the analysis and design of optimization algorithms
[15, 16, 17, 18, 19]. From a dynamical system point of view, dissipativity
theory characterizes the manner in which energy dissipates within the system
and drives it towards equilibrium. It provides a direct way to construct a
Lyapunov function, which further relates the rate of decrease of this internal
energy to the rate of convergence of the algorithm.

We motivate our developments by looking first at a simple scalar bilinear
problem wherein, the energy of the system, expressed as the square 2-norm
distance to the saddle, strictly increases on every iteration, leading to
oscillations of increasing amplitude. To tackle this unstable oscillating
behavior, we propose the Dissipative GDA method, which, as the name suggests,
incorporates a simple friction term to GDA updates to dissipate the internal
energy and stabilize the system. Our algorithm can be seen as a discrete-time
version of [20], which has been applied to solve the C-RL problems [10]. In
this work, we build on this literature, making the following contributions:

_1\. Novel control theory inspired algorithm:_ We illustrate how to use
control theoretic concepts of dissipativity theory to design an algorithm that
can stabilize the unstable behavior of GDA. Particularly, we show that by
introducing a friction term, the proposed DGDA algorithm dissipates the stored
internal energy and converges toward equilibrium.

_2\. Theoretical analysis with better rates:_ We establish the linear
convergence of the DGDA method for bilinear and strongly convex-strongly
concave saddle point problems. In both settings, we show that the DGDA method
outperforms other state-of-the-art first-order explicit methods, surpassing
the standard known linear convergence rate (see Table I and II).

_3\. Numerical Validation:_ We corroborate our theoretical results with
numerical experiments by evaluating the performance of the DGDA method with
GDA, EG, and OGDA methods. When applied to solve bilinear and strongly convex-
strongly concave saddle point problems, the DGDA method systematically
outperforms other methods in terms of convergence rate.

_Outline:_ The rest of the paper is organized as follows. In Section II, we
provide some preliminary definitions and background. In Section III, we
leverage tools from dissipativity theory and propose the Dissipative GDA
(DGDA) algorithm to tackle the unstable oscillatory behavior of GDA methods.
In Section IV, we establish its linear convergence rate for bilinear and
strongly convex-strongly concave problems, which outperforms state-of-the-art
first-order explicit algorithms, including GDA, EG, and OGDA methods. In
Section V, we support our claims with two numerical examples. We close the
paper with concluding remarks and future research directions in Section VI.

Bil. |  Mokhtari, 2020 |  Azizian, 2020 |  This Work
---|---|---|---
EG |  $\frac{\kappa^{-1}}{20}$ |  $\frac{\kappa^{-1}}{64}$ |  -
OG |  $\frac{\kappa^{-1}}{800}$ |  $\frac{\kappa^{-1}}{128}$ |  -
DG |  - |  - |  $\frac{\kappa^{-1}}{4}$

TABLE I: Summary of the global convergence results for EG, OGDA, and DGDA
methods with bilinear objective functions. If a result shows that the iterates
converge as $\mathcal{O}((1-r)^{t})$, the quantity r is reported (the larger
the better). $\kappa$ represents the condition number.

##  II Problem Formulation

In this paper, we study the problem of finding saddle points in the min-max
optimization problem:

$\displaystyle\min_{x\in\mathbb{R}^{n}}\max_{y\in\mathbb{R}^{m}}f(x,y),$ (1)

where the function $f:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ is a
convex-concave function. Precisely, $f(\cdot,y)$ is convex for all
$y\in\mathbb{R}^{m}$ and $f(x,\cdot)$ is concave for all $x\in\mathbb{R}^{n}$.
We seek to develop a novel optimization algorithm that converges to some
saddle point $(x^{*},y^{*})$ of Problem 1.

######  Definition 1 (Saddle Point)

A point $(x^{*},y^{*})\in\mathbb{R}^{n}\times\mathbb{R}^{m}$ is a saddle point
of convex-concave function (1) if and only if it satisfies $f(x^{*},y)\leq
f(x^{*},y^{*})\leq f(x,y^{*})$ for all
$x\in\mathbb{R}^{n},y\in\mathbb{R}^{m}$.

Throughout this paper, we consider two specific instances of Problem 1
commonly studied in related literature: strongly convex-strongly concave and
bilinear functions. Herein, we briefly present some definitions and
properties.

######  Definition 2 (Strongly Convex)

A differentiable function $f:\mathbb{R}^{n}\to\mathbb{R}$ is said to be
$\mu$-strongly convex if $f(w)\geq f(w^{\prime})+\nabla
f(w)^{T}(w-w^{\prime})+\frac{\mu}{2}\|w-w^{\prime}\|^{2}$.

Notice that if $\mu=0$, then we recover the definition of convexity for a
continuously differentiable function and $f(w)$ is $\mu$-strongly concave if
$-f(w)$ is $\mu$-strongly convex. Another important property commonly used in
the convergence analysis of optimization algorithms is the Lipschitz-ness of
the gradient $\nabla f(w)$.

######  Definition 3 ($L$-Lipschitz)

A function $F:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is L-Lipschitz if $\forall
w,w^{\prime}\in\mathbb{R}^{n}$, we have $\|F(w)-F(w^{\prime})\|\leq
L\|w-w^{\prime}\|$.

Combining the above two properties leads to the first important class of
problem that has been extensively studied [21, 22, 2, 1].

###### Assumption 1

(Strongly strongly convex-concave functions with L-Lipschitz Gradient) The
function $f:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ is continuously
differentiable, $\mu$ strongly convex in $x$, and $\mu$ strongly concave in
$y$. Further, the gradient vector $(\nabla_{x}f(x,y);-\nabla_{y}f(x,y))$ is
$L$-Lipschitz.

It is also crucial to consider situations where the objective function is
bilinear. Such bilinear min-max problems often appear when solving constrained
reinforcement learning problems [10, 23], and training of WGANs [12].

######  Assumption 2 (Bilinear function)

The function $f:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ is a bilinear
function if it can be written in the form $f(x,y)=x^{T}Ay$. For simplicity, we
further assume that the matrix $A\in\mathbb{R}^{m\times n}$ is full rank, with
$m\leq n$.

As seen in Table I and II as well as in Section IV, the linear convergence
rates of existing algorithms are frequently characterized by the condition
number $\kappa$. Specifically, when the objective function is bilinear, the
condition number is defined as
$\kappa:=\sigma_{\max}^{2}(A)/\sigma_{\min}^{2}(A)$, where $\sigma_{\max}(M)$
and $\sigma_{\min}(M)$ denote the largest singular value and smallest singular
of a matrix $M$ respectively. When the objective function is strongly convex-
strongly concave with the $L$-Lipschitz gradient, the condition number of the
problem is defined as $\kappa:=L/\mu$.

##  III Dissipative Gradient Descent Ascent Algorithm

S.C |  Zhang, 2021 |  Mokhtari, 2020 |  Azizian, 2020 |  This Work
---|---|---|---|---
GD |  $\kappa^{-2}$ |  - |  - |  -
EG |  - |  $\frac{\kappa^{-1}}{4}$ |  $\frac{\kappa^{-1}}{4}+\epsilon$ |  -
OG |  - |  $\frac{\kappa^{-1}}{4}$ |  $\frac{\kappa^{-1}}{4}+\epsilon$ |  -
DG |  - |  - |  - |   $\kappa^{-1}\\!\\!\\!~{}-~{}\mathcal{O}(\kappa^{-2})$

TABLE II: Summary of the global convergence results for GDA, EG, OGDA, and
DGDA methods with strongly convex-strongly concave and $L$-Lipschitz objective
functions. The table reports the term $r$ of a $(1-r)$ linear rate. The
constant $\epsilon>0$ depends on the problem.

This section introduces the proposed first-order method for solving the min-
max optimization problem 1. The algorithm can be seen as a discretization of
the algorithm proposed by [20], where the authors introduce a regularization
framework for continuous saddle flow dynamics that guarantees asymptotic
convergence to a saddle point under mild assumptions. However, the continuous-
time analysis presented in [20] does generally extend to discrete time. In
this work, we will show the linear convergence of the discrete-time version of
this algorithm. Moreover, we illustrate a novel control-theoretic design
methodology for optimization algorithms that can broadly use dissipativity
tools to study the convergence of the resulting algorithm.

Our results build on gaining an intuitive understanding of the problems that
one encounters when applying the vanilla Gradient Descent Ascent method to
solve saddle point problems (1):
Gradient Descent Ascent (GDA)

$\displaystyle x_{k+1}=x_{k}-\eta\nabla_{x}f(x_{k},y_{k}),$ (2) $\displaystyle
y_{k+1}=y_{k}+\eta\nabla_{y}f(x_{k},y_{k}).$

When (1) is strongly convex-strongly concave, and has L-Lipschitz gradients,
the GDA method provides linear convergence, with step size $\eta=\mu/L^{2}$
and a know rate estimate of $1-1/\kappa^{2}$ [24]. However, when the problem
is bilinear, the standard GDA method fails to converge.

###  III-A Control Theory-Based Motivation

Before we dive into our algorithm, we would like to illustrate our key insight
from control theory about how to characterize the energy stored in the system
and how to introduce friction to dissipate energy and drive the system toward
equilibrium. To get started, consider the following motivating example with a
simple bilinear objective function:

$\displaystyle\min_{x}\max_{y}f(x,y):=xy$ (3)

We seek to analyze the behavior of (a possibly controlled) vanilla GDA
algorithm applied to the above simple bilinear objective function, that is:

$z_{k+1}=(I-\eta M)z_{k}+u_{k}$ (4)

where $z_{k}=(x_{k},y_{k})\in\mathbb{R}^{2}$ represent the system state,
$u_{k}\in\mathbb{R}^{2}$ a controlled input, and $M:=[\,0\,\,1\,;-1\,\,0\,]$.
In the absence of control ($u_{k}=0$, $\forall k$), the algorithm diverges as
illustrated by the upper subplot of Figure 1. To understand this unstable
process, we keep track of the system’s energy, expressed as the square 2-norm
distance to the unique saddle point $(0,0)$,
$V(z_{k})=\|z_{k}\|^{2}=x_{k}^{2}+y_{k}^{2}$. A simple algebraic manipulation
shows that

$\displaystyle V(z_{k+1})$
$\displaystyle=(1+\eta^{2})V(z_{k})+2u_{k}^{T}z_{k}+\|u_{k}\|^{2}.$ (5)

Notably, (5) shows how, without control, the energy grows exponentially with
each iteration.

Figure 1: Trajectories of states for GDA and DGDA for the simple bilinear
objective function $f(x,y):=xy$.

We to choose $u_{k}$ so that the right-hand side of (5) is upper-bounded by
$\alpha^{2}V(z_{k})$, with $|\alpha|<1$. However, it is usually difficult to
do without knowledge of the equilibrium.

The proposed approach uses state augmentation to incorporate a dissipation
component that receives some external input $\hat{u}_{k}$ and filters its
fluctuations, i.e.,

$\hat{z}_{k+1}=(1-\rho)\hat{z}_{k}+\rho\hat{u}_{k}.$ (6)

The dissipation interpretation of (6) can be readily seen, again, by looking
at the evolution of its energy

$\displaystyle
V(\hat{z}_{k+1})\\!=\\!(1\\!-\\!\rho)^{2}V(\hat{z}_{k})\\!+\\!2\rho(1\\!-\\!\rho)\hat{z}_{k}^{T}\hat{u}_{k}\\!+\\!\rho^{2}\|\hat{u}_{k}\|^{2}$
(7)

which has the natural tendency to dissipate energy when $|1-\rho|<1$ and
$\hat{u}_{k}=0$, $\forall k$.

The key innovation of the proposed approach is to combine (4) and (6) with a
proper choice of $u_{k}$ and $\hat{u}_{k}$ that satisfies:

  * •

The combined energy $V(\xi_{k})=\|z_{k}\|^{2}+\|\hat{z}_{k}\|^{2}$ decreases,
i.e., there is $|\alpha|<1$ such that:

$V(\xi_{k+1})\leq\alpha^{2}V(\xi_{k})\,.$ (8)

  * •

The iterates $u_{k}\rightarrow u^{*}=0$, so that the equilibrium of (4), that
is, the saddle point does not change.

To achieve these goals we select

$u_{k}=-\rho(z_{k}-\hat{z}_{k}),\quad\text{ and }\quad\hat{u}_{k}=z_{k}\,,$
(9)

thus leading to the combined algorithm

$\displaystyle z_{k+1}$ $\displaystyle=z_{k}-\eta
Mz_{k}-\rho(z_{k}-\hat{z}_{k})$ (10) $\displaystyle\hat{z}_{k+1}$
$\displaystyle=\hat{z}_{k}-\rho(\hat{z}_{k}-z_{k})$ (11)

The choice of $\hat{u}_{k}$ ensures that $\hat{z}_{k}$ acts as a filtering
process of $z_{k}$. Further, since for constant $z_{k}=\bar{z}$,
$\hat{z}_{k}\rightarrow\bar{z}$, the choice of $u_{k}$ ensures that the
equilibrium of (4) remains unchanged. We will show in Section IV that one can
choose $\eta$ and $\rho$ so that (8) holds. We refer to the bottom of Figure 1
for an illustration of its effectiveness.

###  III-B Dissipative GDA Algorithm

In the previous section, we showed how to design a dissipative component (6)
to stabilize the GDA updates on a simple bilinear objective. Since the design
approach is independent of the particular saddle function $f(x,y)=xy$, it can
be readily applied in more general settings. The algorithm presented below
seeks to incorporate a similar friction component into a standard GDA update.
This modification addresses the inherent limitation of the GDA update when
applied to general bilinear problems.

Dissipative gradient descent ascent (DGDA):

$\displaystyle\begin{bmatrix}x_{k+1}\\\ \hat{x}_{k+1}\\\ y_{k+1}\\\
\hat{y}_{k+1}\end{bmatrix}=\begin{bmatrix}x_{k}-\eta\nabla_{x}f(x_{k},y_{k})-\rho(x_{k}-\hat{x}_{k})\\\
\hat{x}_{k}-\rho(\hat{x}_{k}-x_{k})\\\
y_{k}+\eta\nabla_{y}f(x_{k},y_{k})-\rho(y_{k}-\hat{y}_{k})\\\
\hat{y}_{k}-\rho(\hat{y}_{k}-y_{k})\end{bmatrix}$ (12)

Particularly, for $f$ as in (1), in (12) we introduce two new sets of
variables $\hat{x}\in\mathbb{R}^{n}$ and $\hat{y}\in\mathbb{R}^{m}$ and a
damping parameter $\rho>0$. One important observation is that, once the system
reaches equilibrium, i.e.,
$x_{k+1}=x_{k},y_{k+1}=y_{k},\hat{x}_{k+1}=\hat{x}_{k},\hat{y}_{k+1}=\hat{y}_{k}$,
one necessarily has $\hat{x}_{k}=x_{k}$ and $\hat{y}_{k}=y_{k}$, which ensures
that the fixed point is necessarily a critical point of the saddle function.

In the remaining part of this section, we will first provide several key
observations and properties of the proposed DGDA updates. Then, in the next
section, we will formally prove that the proposed DGDA algorithm provides a
linear convergence guarantee for both bilinear and strongly convex-strongly
concave functions. Mainly, we will show that by simply adding a friction or
damping term, the DGDA updates succeed in dissipating the energy of the system
even in cases where GDA does not.

###  III-C Key Properties and Related Algorithms

The first important observation is that the above DGDA update could be
considered as applying a vanilla GDA update to the following regularized
surrogate for $f(x,y)$:

$\displaystyle
f(x,y,\hat{x},\hat{y})\\!:=\\!f(x,y)\\!+\\!\frac{\rho}{2}\|x\\!-\\!\hat{x}\|^{2}\\!-\\!\frac{\rho}{2}\|y\\!-\\!\hat{y}\|^{2}.$
(13)

Note that this is different from the Proximal Point Method [21, 25] or
introducing a $L_{2}$ regularization [26, 27].

Differences with $L_{2}$ regularization A commonly used method to ensure
convergence is to introduce a $L_{2}$ regularization term in $x$ and $y$ [26,
27]:

$\displaystyle\min_{x\in\mathbb{R}^{n}}\max_{y\in\mathbb{R}^{m}}f(x,y)+\frac{\rho}{2}\|x\|^{2}-\frac{\rho}{2}\|y\|^{2}.$
(14)

Although the augmented objective function becomes strongly convex-strongly
concave and the vanilla GDA updates will converge, this regularization changes
the saddle points. While our algorithm also introduces two regularizing terms,
the following Lemma verifies the fixed positions of saddle points between
$f(x,y)$ and $f(x,y,\hat{x},\hat{y})$ with virtual variables aligned with
original variables.

######  Lemma 1 (Saddle Point Invariance)

[20, Lemma 6] For problem 1, a point $(x^{*},y^{*})$ is a saddle point of
$f(x,y)$ if and only if $(x^{*},y^{*},\hat{x}^{*},\hat{y}^{*})$ is a saddle
point of $f(x,y,\hat{x},\hat{y})$, with $\hat{x}^{*}=x^{*}$ and
$\hat{y}^{*}=y^{*}$.

More interestingly, the regularization terms,
$\frac{\rho}{2}\|x-\hat{x}\|^{2}$ and $\frac{\rho}{2}\|y-\hat{y}\|^{2}$, do
not introduce extra strong convexity-stong concavity to the original problem.
Precisely, the augmented problem $f(x,y,\hat{x},\hat{y})$ is neither strongly
convex on ($x,\hat{x}$) nor strongly concave on $(y,\hat{y})$. Indeed, on the
hyperplane of $x=\hat{x}$ and $y=\hat{y}$, the augmented problem recovers the
original problem $f(x,y,\hat{x},\hat{y})=f(x,y)$. Additionally, the introduced
regularization terms of the DGDA method are separable and local, thus
preserving the distributed structure that original systems may have.
Consequently, it can be seamlessly integrated into a fully distributed
approach.

Differences with Proximal Point Method The Proximal Point Method for saddle
point problems [21] shares a similar structure with DGDA algorithm. In the
proximal method, the next iterates $(x_{k+1},y_{k+1})$ is the unique solution
to the saddle point problem

$\displaystyle(x_{k+1},y_{k+1})=\mathrm{prox}_{\eta}(x_{k},y_{k})$ (15)
$\displaystyle\\!:=\\!\arg\min_{x\in\mathbb{R}^{n}}\max_{y\in\mathbb{R}^{m}}f(x,y)\\!+\\!\frac{\eta}{2}\|x\\!-\\!x_{k}\|^{2}\\!-\\!\frac{\eta}{2}\|y\\!-\\!y_{k}\|^{2}$

Using the optimality conditions of (15), the update of the Proximal Point
method can be written as:

$\displaystyle x_{k+1}=x_{k}-\eta\nabla_{x}f(x_{k+1},y_{k+1}),$ (16)
$\displaystyle y_{k+1}=y_{k}+\eta\nabla_{y}f(x_{k+1},y_{k+1}).$

This expression shows that the Proximal Point method is an implicit method.
Although Implicit methods are known to be more stable and to benefit from
better convergence properties [25, 6], implementing the above updates requires
computing the operators $(I+\eta\nabla_{x}f)^{-1}$ and
$(I+\eta\nabla_{y}f)^{-1}$, and therefore may be computationally intractable.
In contrast, the DGDA algorithm is an explicit algorithm, which applies a
vanilla GDA update to the augmented objective function $\eqref{eq:augmented
saddle flow}$. Thus, as shown in Table I and II, and in the next section, we
choose to compare the convergence of DGDA only with other explicit algorithms
for saddle point problems, such as the GDA, Extra-Gradient (EG), and
Optimistic GDA (OGDA) methods which have comparable per-iteration
computational requirements; although the DGDA algorithm has twice as many
state variables, it only requires a single gradient computation per update.
Moreover, there is no need for retaining and reemploying the extrapolated
gradient, which also sets it apart from the OGDA method.

We finalize this section comparing DGDA with recent efforts leveraging Moreau-
Yosida smoothing techniques to solve nonconvex-concave [28, 29, 30],
nonconvex-nonconcave [31] min-max optimization problems.

Smoothed-GDA method The Smoothed-GDA consists of the following updates

$\displaystyle x_{k+1}=x_{k}-c\nabla_{x}K(x_{k},z_{k};y_{k}),$ (17)
$\displaystyle y_{k+1}=y_{k}+\alpha\nabla_{y}K(x_{k+1},z_{k};y_{k})$
$\displaystyle z_{k+1}=z_{k}+\beta(x_{k+1}-z_{k})$

where $K(x,z;y)=f(x,y)+\frac{p}{2}\|x-z\|^{2}.$ Smoothed-GDA was independently
introduced by [32]. It was originally motivated as an ADMM algorithm to solve
the linearly constrained nonconvex differentiable minimization problem [32].
The algorithm has been further extended to consider nonconvex-concave min-max
optimization problems [28]. Despite its similarities with DGDA, the iterates
of smoothed-GDA are sequential, requiring first updates in $x$ and
subsequently in $(y,z)$, whereas DGDA updates are fully parallel or
synchronized. Unfortunately, the setting where such algorithm has been studied
is different from the one considered in this paper, which difficult the
comparison with the present work. A thorough comparison of DGDA with
(LABEL:eq:_SGDA) is subject of current research.

##  IV Convergence Analysis

In this section, we provide a theoretical analysis of the proposed algorithm.
For the purpose of our analysis, we consider a quadratic Lyapunov function to
track the energy dissipation of the DGDA updates

$V_{k}:=\|x_{k}-x^{*}\|^{2}+\|y_{k}-y^{*}\|^{2}+\|\hat{x}_{k}-\hat{x}^{*}\|^{2}+\|\hat{y}_{k}-\hat{y}^{*}\|^{2},$

which denotes the square 2-norm distance to the saddle point at the $k$-th
iteration. The goal is, therefore, to find some $0to\leq\alpha<1$ such that
$V_{k+1}\leq\alpha V_{k}$, where $\alpha$ denotes the linear convergence rate.

###  IV-A Convergence Analysis for Bilinear Functions

When applied to the bilinear min-max optimization problem $f(x,y)=x^{T}Ay$,
the DGDA update (12) is equivalent to a linear dynamical system. Specifically,
denote $z=[x,y]^{T},\hat{z}=[\hat{x},\hat{y}]^{T}$ yields:

$\displaystyle\begin{bmatrix}z_{k+1}\\!-\\!z^{*}\\\
\hat{z}_{k+1}\\!-\\!\hat{z}^{*}\\\
\end{bmatrix}\\!=\\!\begin{bmatrix}(1\\!-\\!\rho)I\\!-\\!\eta
M\\!&\\!\\!\\!\rho I\\\ \rho I\\!&\\!\\!\\!(1\\!-\\!\rho)I\\\
\end{bmatrix}\\!\begin{bmatrix}z_{k}\\!-\\!z^{*}\\\
\hat{z}_{k}\\!-\\!\hat{z}^{*}\\\ \end{bmatrix},$ (18)

where $M=\begin{bmatrix}\mathbf{0}&A\\\ -A^{T}&\mathbf{0}\\\ \end{bmatrix}.$
As a result, the linear convergence of DGDA, as well as its convergence rate
can be derived from the analysis of the spectrum of the associated matrix that
defines the DGDA update in (18). This yields the following theorem.

###### Theorem 2

(Linear convergence of DGDA, Bilinear Case) Let Assumption 2 hold. Then the
updates 12 of DGDA with $0<\eta\leq\frac{2\rho}{\sigma_{\max}(A)}$ and
$\rho>0$ provide linearly converging iterates. That is, there is a constant
$\beta>0$, independent of $V_{0}$ such that

$\displaystyle\textstyle
V_{k}\\!\leq\\!\mathcal{O}\left(\left(1\\!-\\!2\rho\\!+\\!2\rho^{2}\\!+\\!(1\\!-\\!\rho)\sqrt{4\rho^{2}\\!-\\!\eta^{2}\sigma_{\min}^{2}(A)}\right)^{k}\right)V_{0},$

Particularly, setting $\rho=1/2$ and $\eta=1/\sigma_{\max}(A)$ we have

$\displaystyle\textstyle
V_{k}\leq\mathcal{O}\left(\left(1-\frac{1}{4\kappa}\right)^{k}\right)\,V_{0}.$
(19)

We remark that linear convergence requires $\rho>0$. This is not surprising
since GDA, which is known to diverge for bilinear functions, can be
interpreted as DGDA method when $\rho=0$. More importantly, by choosing the
optimal step size $\rho=1/2,\eta=1/\sigma_{\max}(A)$, DGDA method achieves a
better linear convergence rate than the EG and OGDA methods (see Table I).

###  IV-B Convergence Analysis for Strongly Convex Stronly Concave Functions

We now consider the case of strongly convex-strongly concave min-max problems.
Let $F(z_{k}):=(\nabla_{x}f(x_{k},y_{k}),-\nabla_{y}f(x_{k},y_{k}))$. The DGDA
updates can be written as follows:

$\displaystyle\begin{bmatrix}z_{k+1}\\\
\hat{z}_{k+1}\end{bmatrix}=\begin{bmatrix}z_{k}-\eta
F(z_{k})-\rho(z_{k}-\hat{z}_{k})\\\
\hat{z}_{k}-\rho(\hat{z}_{k}-z_{k})\end{bmatrix}$ (20)

Because of the existence of the nonlinear term $F(z_{k})$, we cannot analyze
the spectrum as in the previous bilinear case. This is indeed a common
challenge in analyzing most optimization algorithms beyond a neighborhood of
the fixed point. We circumvent this problem by leveraging recent results on
the analysis of variational mappings as $F(\cdot)$ via integral quadratic
constraint [18, 16, 17]. More details can be found in the Appendix where we
prove the following theorem.

###### Theorem 3

(Linear convergence of DGDA, Strongly Convex-Strongly Concave Case) Let
Assumption 1 hold, then the updates (12) with $\rho=1/2$ and $\eta=1/(L+\mu)$
of the DGDA algorithm provide linearly converging iterates:

$\displaystyle
V_{k}\leq\biggl{(}1-\kappa^{-1}+\mathcal{O}(\kappa^{-2})\biggr{)}^{k}V_{0}$
(21)

Similarly, as in the bilinear case, we remark on the importance of the
dissipation component. When $\rho=0$, a similar analysis as in the proof of
the theorem recovers the lower bound of the convergence rate of GDA
$(1-\kappa^{-2})$ as shown in [19, 3.1]. Thus, our DGDA method provides a
better convergence rate estimate than GDA, since clearly
$\kappa\in[1,\infty)$, and therefore $\kappa^{-2}\leq\kappa^{-1}$.

We remark that while the rate obtained in Theorem 3 is clearly better than
those of the EG and OGDA methods for large condition numbers $\kappa$ (see
Table II), the theorem fails to quantify the comparative performance of DGDA
for small values of $\kappa$. The following corollary shows that indeed, the
rate of DGDA is provably better for all $\kappa\geq 2$.

######  Corollary 4 (SCSC, comparison with known rates)

Let Assumption 1 hold, and suppose that $L\geq 2m$, i.e., $\kappa\geq 2$.
Then, the linear convergence rate estimate of DGDA (21) is smaller (better)
than that of of EG and OGDA, i.e., $1-\kappa^{-1}/4$ (Theorem $6\&7$ [2] and
Theorem $4\&7$ [1]).

##  V Numerical Experiments

Figure 2: Convergence of GDA, EG, OGDA, and DGDA in terms of the number of
gradient evaluations for the bilinear problem. GDA diverges and the error is
not shown. All other three algorithms converge linearly, where the DGDA method
provides the best performance.

In this section, we compare the performance of the proposed Dissipative
gradient descent (DGDA) method with the Extra-gradient (EG), Gradient descent
ascent (GDA), and Optimistic gradient descent ascent (OGDA) methods.

###  V-A Bilinear problem

We first consider the following bilinear min-max optimization problem:
$\min_{x\in\mathbb{R}^{n}}\max_{y\in\mathbb{R}^{m}}x^{T}Ay$, where
$A\in\mathbb{R}^{m\times n}$ is full-rank. The simulation results are
illustrated in Figure 2 and Figure 3. In this experiment, we set the dimension
of the problem to $m=n=10$ and the iterates are initialized at $x_{0},y_{0}$,
which are randomly drawn from the uniform distribution on the open interval
$(0,1)$.

We plot the errors (distance to saddle points) of DGDA, EG, and OGDA versus
the number of gradient evaluations for this problem in Figure 2. The solid
line and grey-shaded error bars represent the average trajectories and
standard deviations of 20 trials, where in each trial the randomly generated
matrix $A$ has a fixed condition number, i.e.,
$\kappa=\sigma^{2}_{\max}(A)/\sigma^{2}_{\min}(A)=25$. The key motivation is
that all three algorithms’ convergence rates critically depend on
$\kappa^{-1}$, and by fixing the condition number, we provide an explicit
comparison of their convergence speed.

We pick the step size for different methods according to theoretical findings.
That is, we select $\rho=1/2$ and $\eta=1/\sigma_{\max}(A)$ for DGDA (Theorem
2), $\eta=1/4L=1/4\sigma_{\max}(A)$ for EG and OGDA (Theorem 6&7 [2] and
Theorem 4&7 [1]). In Figure 2, we do not show the error of GDA since it
diverges for this bilinear saddle point problem. All other three algorithms
converge linearly, with the DGDA method providing the best performance.

Finally, to provide a qualitative demonstration of how DGDA fares with other
existing algorithms, we further plot the sample trajectories of GDA, EG, OGDA,
and EGDA on a simple 2D bilinear min-max problem, with $m=n=1$. In Figure 3,
we observe that while GDA diverges, the trajectories of all other three
algorithms converge linearly to the saddle point $(x^{*},y^{*})=(0,0)$.
Interestingly, our proposed algorithm (DGDA) despite taking larger steps,
exhibits faster linear convergence.

Figure 3: Trajectories of GDA, EG, OGDA, and DGDA for a 2d bilinear problem.
GDA diverges and all other three algorithms converge linearly, where the DGDA
method provides the best performance.

###  V-B Strongly convex-strongly concave problem

In the second numerical example, we focus on a strongly convex-strongly
strongly-concave quadratic problem of the following form:

$\displaystyle\min_{x\in\mathbb{R}^{n}}\max_{y\in\mathbb{R}^{m}}\frac{1}{2}x^{T}Ax-\frac{1}{2}y^{T}By+x^{T}Cy,$
(22)

where the matrices satisfy $\mu_{A}I\preceq A\preceq L_{A}I$, $\mu_{B}I\preceq
B\preceq L_{B}I$, $\mu_{c}^{2}I\preceq C^{T}C\preceq L_{c}^{2}I$. As a result,
the problem (22) satisfy Assumption 1. In this experiment, we set the
dimension of the problem to $n=50,m=10$, and the iterates are initialized at
$x_{0},y_{0}$, which are randomly drawn from the uniform distribution on the
open interval $(0,1)$. We plot the errors (distance to saddle points) of GDA,
DGDA, EG, and OGDA versus the number of gradient evaluations for this problem
in Figure 4. Again, the solid line and grey-shaded error bars represent the
average trajectories and standard deviations of 20 trials, where in each trial
the randomly generated matrix $\begin{bmatrix}A&C\\\ -C^{T}&B\end{bmatrix}$ is
chosen such that the condition number of (22) remains constant, i.e.,
$\kappa=L/\mu=31$. Similarly as in the bilinear problem in Section V-A, we
pick the step size for the DGDA method according to our theoretical finding in
Theorem 3. The step size of the GDA method is selected as $\eta=\mu/L^{2}$
(Theorem 5 [33]). The step sizes for EG and OGDA methods are selected as
$\eta=1/4L$ (Theorem 6&7 [2] and Theorem 4&7 [1]). According to the plots, all
algorithms converge linearly, and the DGDA method has the best performance.

Figure 4: Convergence of GDA, EG, OGDA, and DGDA in terms of the number of
gradient evaluations for problem 22. All algorithms converge linearly, and the
DGDA method has the best performance.

##  VI Conclusion and Future Work

In this work, we present the Dissipative GDA (DGDA) algorithm, a novel method
for solving min-max optimization problems. Drawing inspiration from
dissipativity theory and control theory, we address the challenge of diverging
oscillations in bilinear min-max optimization problems when using the Gradient
Descent Ascent (GDA) method. Particularly, we introduce a friction term into
the GDA updates aiming to dissipate the internal energy and drive the system
towards equilibrium. By incorporating a state-augmented regularization, our
proposed DGDA method can be seen as performing standard GDA on an extended
saddle function without introducing additional convexity. We further establish
the superiority of the convergence rate of the proposed DGDA method when
compared with other established methods including GDA, Extra-Gradient (EG),
and Optimistic GDA. The analysis is further supported by two numerical
examples, demonstrating its effectiveness in solving saddle point problems.
Our future work includes studying the DGDA method in a stochastic setting and
its application in solving Constrained Reinforcement learning problems in the
policy space.

## References

  * [1] A. Mokhtari, A. Ozdaglar, and S. Pattathil, “A unified analysis of extra-gradient and optimistic gradient methods for saddle point problems: Proximal point approach,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2020, pp. 1497–1507.
  * [2] W. Azizian, I. Mitliagkas, S. Lacoste-Julien, and G. Gidel, “A tight and unified analysis of gradient-based methods for a whole spectrum of differentiable games,” in _International conference on artificial intelligence and statistics_. PMLR, 2020, pp. 2863–2873.
  * [3] E. Gorbunov, H. Berard, G. Gidel, and N. Loizou, “Stochastic extragradient: General analysis and improved rates,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2022, pp. 7865–7901.
  * [4] N. Loizou, H. Berard, G. Gidel, I. Mitliagkas, and S. Lacoste-Julien, “Stochastic gradient descent-ascent and consensus optimization for smooth games: Convergence analysis under expected co-coercivity,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 19 095–19 108, 2021.
  * [5] A. Beznosikov, E. Gorbunov, H. Berard, and N. Loizou, “Stochastic gradient descent-ascent: Unified theory and new efficient methods,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2023, pp. 172–235.
  * [6] G. Gidel, H. Berard, G. Vignoud, P. Vincent, and S. Lacoste-Julien, “A variational inequality perspective on generative adversarial networks,” _arXiv preprint arXiv:1802.10551_ , 2018.
  * [7] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” _Advances in neural information processing systems_ , vol. 27, 2014.
  * [8] C. Daskalakis, A. Ilyas, V. Syrgkanis, and H. Zeng, “Training gans with optimism,” _arXiv preprint arXiv:1711.00141_ , 2017.
  * [9] D. Pfau and O. Vinyals, “Connecting generative adversarial networks and actor-critic methods,” _arXiv preprint arXiv:1610.01945_ , 2016.
  * [10] T. Zheng, P. You, and E. Mallada, “Constrained reinforcement learning via dissipative saddle flow dynamics,” in _2022 56th Asilomar Conference on Signals, Systems, and Computers_. IEEE, 2022, pp. 1362–1366.
  * [11] D. Ding, C.-Y. Wei, K. Zhang, and A. Ribeiro, “Last-iterate convergent policy gradient primal-dual methods for constrained mdps,” _arXiv preprint arXiv:2306.11700_ , 2023.
  * [12] J. Adler and S. Lunz, “Banach wasserstein gan,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [13] E. Altman, _Constrained Markov decision processes: stochastic modeling_. Routledge, 1999.
  * [14] S. Sastry, _Nonlinear systems: analysis, stability, and control_. Springer Science & Business Media, 2013, vol. 10.
  * [15] Z. E. Nelson and E. Mallada, “An integral quadratic constraint framework for real-time steady-state optimization of linear time-invariant systems,” in _2018 annual American control conference (ACC)_. IEEE, 2018, pp. 597–603.
  * [16] B. Hu and L. Lessard, “Dissipativity theory for nesterov’s accelerated method,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 1549–1557.
  * [17] M. Fazlyab, A. Ribeiro, M. Morari, and V. M. Preciado, “Analysis of optimization algorithms via integral quadratic constraints: Nonstrongly convex problems,” _SIAM Journal on Optimization_ , vol. 28, no. 3, pp. 2654–2689, 2018.
  * [18] L. Lessard, B. Recht, and A. Packard, “Analysis and design of optimization algorithms via integral quadratic constraints,” _SIAM Journal on Optimization_ , vol. 26, no. 1, pp. 57–95, 2016.
  * [19] G. Zhang, X. Bao, L. Lessard, and R. Grosse, “A unified analysis of first-order methods for smooth games via integral quadratic constraints,” _The Journal of Machine Learning Research_ , vol. 22, no. 1, pp. 4648–4686, 2021.
  * [20] P. You and E. Mallada, “Saddle flow dynamics: Observable certificates and separable regularization,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 4817–4823.
  * [21] R. T. Rockafellar, “Monotone operators and the proximal point algorithm,” _SIAM journal on control and optimization_ , vol. 14, no. 5, pp. 877–898, 1976.
  * [22] P. Tseng, “On linear convergence of iterative methods for the variational inequality problem,” _Journal of Computational and Applied Mathematics_ , vol. 60, no. 1-2, pp. 237–252, 1995.
  * [23] M. Wang, “Randomized linear programming solves the markov decision problem in nearly linear (sometimes sublinear) time,” _Mathematics of Operations Research_ , vol. 45, no. 2, pp. 517–546, 2020.
  * [24] B. Grimmer, H. Lu, P. Worah, and V. Mirrokni, “The landscape of the proximal point method for nonconvex–nonconcave minimax optimization,” _Mathematical Programming_ , vol. 201, no. 1-2, pp. 373–407, 2023.
  * [25] N. Parikh, S. Boyd _et al._ , “Proximal algorithms,” _Foundations and trends® in Optimization_ , vol. 1, no. 3, pp. 127–239, 2014.
  * [26] A. Krogh and J. Hertz, “A simple weight decay can improve generalization,” _Advances in neural information processing systems_ , vol. 4, 1991.
  * [27] T. Hastie, R. Tibshirani, J. H. Friedman, and J. H. Friedman, _The elements of statistical learning: data mining, inference, and prediction_. Springer, 2009, vol. 2.
  * [28] J. Zhang, P. Xiao, R. Sun, and Z. Luo, “A single-loop smoothed gradient descent-ascent algorithm for nonconvex-concave min-max problems,” _Advances in neural information processing systems_ , vol. 33, pp. 7377–7389, 2020.
  * [29] Z. Xu, H. Zhang, Y. Xu, and G. Lan, “A unified single-loop alternating gradient projection algorithm for nonconvex–concave and convex–nonconcave minimax problems,” _Mathematical Programming_ , pp. 1–72, 2023.
  * [30] J. Yang, A. Orvieto, A. Lucchi, and N. He, “Faster single-loop algorithms for minimax optimization without strong concavity,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2022, pp. 5485–5517.
  * [31] T. Zheng, L. Zhu, A. M.-C. So, J. Blanchet, and J. Li, “Universal gradient descent ascent method for nonconvex-nonconcave minimax optimization,” _Advances in Neural Information Processing Systems_ , vol. 36, 2024.
  * [32] J. Zhang and Z.-Q. Luo, “A proximal alternating direction method of multiplier for linearly constrained nonconvex minimization,” _SIAM Journal on Optimization_ , vol. 30, no. 3, pp. 2272–2302, 2020.
  * [33] A. Beznosikov, B. Polyak, E. Gorbunov, D. Kovalev, and A. Gasnikov, “Smooth monotone stochastic variational inequalities and saddle point problems: A survey,” _European Mathematical Society Magazine_ , no. 127, pp. 15–28, 2023.
  * [34] G. M. Korpelevich, “The extragradient method for finding saddle points and other problems,” _Matecon_ , vol. 12, pp. 747–756, 1976.
  * [35] R. T. Rockafellar, “Augmented lagrangians and applications of the proximal point algorithm in convex programming,” _Mathematics of operations research_ , vol. 1, no. 2, pp. 97–116, 1976.
  * [36] G. Zhang and Y. Yu, “Convergence of gradient methods on bilinear zero-sum games,” in _International Conference on Learning Representations_ , 2019.
  * [37] J. P. Hespanha, _Linear systems theory_. Princeton university press, 2018.

###  -A First Order Algorithms for Saddle Point Problems

In this section, we introduce several popular first-order methods for solving
the min-max problem 1 in the machine learning community. Precisely, we focus
on Gradient Descent-Ascent (GDA), Extra-gradient (EG), Optimistic Gradient
Descent-Ascent (OGDA), Proximal Point and Smoothed-GDA methods.

#### Gradient descent ascent (GDA)

$\displaystyle x_{k+1}=x_{k}-\eta\nabla_{x}f(x_{k},y_{k}),$ (23)
$\displaystyle y_{k+1}=y_{k}+\eta\nabla_{y}f(x_{k},y_{k}).$

When the problem is strongly convex-strongly concave and L-Lipschitz, the GDA
method provides linear convergence, with step size $\eta=\mu/L^{2}$ and a know
rate estimate of $1-1/\kappa^{2}$ [24, 19]. However, when the problem is
bilinear, the standard GDA method fails to converge. Therefore, variants of
the gradient method such as Extra-Gradient and Optimistic Gradient Descent-
Ascent methods have attracted much attention in recent literature because of
their superior empirical performance in solving min-max optimization problems
such as training GANs and solving C-RL problems.

#### Extra-gradient (EG)

$\displaystyle x_{k+1/2}=x_{k}-\eta\nabla_{x}f(x_{k},y_{k}),$ (24)
$\displaystyle y_{k+1/2}=y_{k}+\eta\nabla_{y}f(x_{k},y_{k}),$ $\displaystyle
x_{k+1}=x_{k}-\eta\nabla_{x}f(x_{k+1/2},y_{k+1/2}),$ $\displaystyle
y_{k+1}=y_{k}+\eta\nabla_{y}f(x_{k+1/2},y_{k+1/2}).$

Extra-gradient is a classical method introduced in [34], where its linear rate
of convergence for bilinear functions and smooth strongly convex-strongly
concave functions have been established in many recent works (see Table I and
II). The Extra-gradient method first computes an extrapolated point
($x_{k+1/2},y_{k+1/2}$) by performing a GDA update. Then the gradients
evaluated at the extrapolated point are used to compute the new iterates
$(x_{k+1},y_{k+1})$.

The linear convergence rate of EG for strongly convex-strongly concave is
established, with a standard known rate of $1-1/4\kappa$; see e.g.[2, 1]. One
issue with the Extra-gradient method is that, as the name suggests, each
update requires evaluation of extra gradients at the extrapolated point
($x_{k+1/2},y_{k+1/2}$), which doubles the computational complexity of EG
method compared to vanilla GDA method.

#### Optimistic gradient descent ascent (OGDA)

$\displaystyle
x_{k+1}\\!=\\!x_{k}-2\eta\nabla_{x}f(x_{k},y_{k})+\eta\nabla_{x}f(x_{k-1},y_{k-1}),$
(25) $\displaystyle
y_{k+1}\\!=\\!y_{k}+2\eta\nabla_{y}f(x_{k},y_{k})-\eta\nabla_{y}f(x_{k-1},y_{k-1}).$

The Optimistic gradient descent ascent (OGDA) method adds a ”negative-
momentum” term to each of the updates, which differentiates the OGDA method
from the vanilla GDA method. Meanwhile, the OGDA method stores and re-uses the
extrapolated gradient for the extrapolation, which only requires a single
gradient computation per update.

The convergence properties of OGDA were also recently investigated in (refer
to Table I and II), demonstrating linear convergence rates with smooth and
bilinear functions, as well as strongly convex-strongly concave functions.

#### Proximal Point (PP)

The proximal point method for convex minimization has been extensively studied
[35, 25] and extended to solve saddle point problems in [21]. Define the
iterates $\\{x_{k+1},y_{k+1}\\}$ as the unique solution to the saddle point
problem

$\displaystyle(x_{k+1},y_{k+1})=\mathrm{prox}_{\eta}(x_{k},y_{k})$ (26)
$\displaystyle:=\arg\min_{x\in\mathbb{R}^{n}}\max_{y\in\mathbb{R}^{m}}f(x,y)+\frac{\eta}{2}\|x-x_{k}\|^{2}-\frac{\eta}{2}\|y-y_{k}\|^{2}$

Using the optimality conditions, the update of the Proximal Point method can
be written as:

$\displaystyle x_{k+1}=x_{k}-\eta\nabla_{x}f(x_{k+1},y_{k+1}),$ (27)
$\displaystyle y_{k+1}=y_{k}+\eta\nabla_{y}f(x_{k+1},y_{k+1}).$

This expression shows that in contrast to explicit methods such as GDA, EG,
and OGDA methods, the Proximal Point method is an implicit method. Although
Implicit methods are known to be more stable and to benefit from better
convergence properties [25, 6], implementing the above updates requires
computing the operators $(I+\eta\nabla_{x}f)^{-1}$ and
$(I+\eta\nabla_{y}f)^{-1}$, and therefore may be computationally intractable.
Notably, in [1], the authors show the EG and OGDA methods can be interpreted
as an approximation of the PP method, and therefore exhibits similar
convergence performance.

#### Smoothed-GDA method

We finalize this section comparing DGDA with recent efforts leveraging Moreau-
Yosida smoothing techniques to solve nonconvex-concave [28, 29, 30],
nonconvex-nonconcave [31] min-max optimization problems.

$\displaystyle x_{k+1}=x_{k}-c\nabla_{x}K(x_{k},z_{k};y_{k}),$ (28)
$\displaystyle y_{k+1}=y_{k}+\alpha\nabla_{K}(x_{k+1},z_{k};y_{k})$
$\displaystyle z_{k+1}=z_{k}+\beta(x_{k+1}-z_{k})$

where $K(x,z;y)=f(x,y)+\frac{p}{2}\|x-z\|^{2}$.

The Smoothed-GDA was independently introduced by Jiawei et al. in [32] and
later [28]. It was originally motivated by ADMM to solve the linearly
constrained nonconvex differentiable minimization problem [32], where they
introduce an extra quadratic proximal term for the equality constraints and an
extra sequence $\\{z_{k}\\}$. They claim this smoothing or exponential
averaging scheme is necessary for the convergence of the proximal ADMM when
the objective is nonconvex. Later on, this scheme is further extended to solve
the nonconvex-concave min-max optimization problem [28].

###  -B Proof of Theorem 2

We consider, for ease of presentation, the case when $A\in\mathbb{R}^{m\times
m}$ is a square matrix. The extension for non-square matrices is
straightforward and has been covered in the literature [36, Appendix G].
Applying the updates 12 to $f(x,y)=x^{T}Ay$ and denoting
$z=[x,y]^{T},\hat{z}=[\hat{x},\hat{y}]^{T}$ yields:

$\displaystyle\begin{bmatrix}z_{k+1}-z^{*}\\\ \hat{z}_{k+1}-\hat{z}^{*}\\\
\end{bmatrix}$ $\displaystyle=\begin{bmatrix}z_{k}-\eta
Mz_{k}-\rho(z_{k}-\hat{z}_{k})\\\ \hat{z}_{k}-\rho(\hat{z}_{k}-z_{k})\\\
\end{bmatrix}$ (29) $\displaystyle=\begin{bmatrix}(1-\rho)I-\eta M&\rho I\\\
\rho I&(1-\rho)I\\\ \end{bmatrix}\begin{bmatrix}z_{k}-z^{*}\\\
\hat{z}_{k}-\hat{z}^{*}\\\ \end{bmatrix},$

where

$\displaystyle M=\begin{bmatrix}\mathbf{0}&A\\\ -A^{T}&\mathbf{0}\\\
\end{bmatrix}$

According to [2, Lemma 7]

$\displaystyle\mathrm{Sp}(M)=\\{\pm
i\sigma|\sigma^{2}\in\mathrm{Sp}(AA^{T})\\}.$

We will use $\sigma_{\max}$ and $\sigma_{\min}$ to denote the largest singular
value and smallest singular of matrix $A$, respectively. And according to
Assumption 2, we have $\sigma_{\min}>0$. Since $M$ is a normal matrix and
diagonalizable, we can compute the eigenvalues of the linear system (29) using
the following similarity transformation

$\displaystyle\begin{bmatrix}(1-\rho)I-\eta M&\rho I\\\ \rho I&(1-\rho)I\\\
\end{bmatrix}=$ (30) $\displaystyle\begin{bmatrix}U^{-1}&\mathbf{0}\\\
\mathbf{0}&U^{-1}\\\
\end{bmatrix}\begin{bmatrix}(1-\rho)I-\eta\boldsymbol{\Lambda}&\rho I\\\ \rho
I&(1-\rho)I\\\ \end{bmatrix}\begin{bmatrix}U&\mathbf{0}\\\ \mathbf{0}&U\\\
\end{bmatrix},$

where $\boldsymbol{\Lambda}=\mathrm{diag}(\lambda_{1},...,\lambda_{2m})$,
$\lambda_{2j-1}=+i\sigma_{j}$, $\lambda_{2j}=-i\sigma_{j}$, with $\pm
i\sigma_{j}\in\mathrm{Sp}(M),j=\\{1,...,m\\}$. In order to show linear
convergence, we want to show that $\max_{j\in[m]}{|\mu_{j}|^{2}}<1$ , where
$\mu_{j}$ are the eigenvalues of the above matrix (30) (i.e., the spectral
radius of matrix (30)). As straight forward computation leads to

$\displaystyle\mu_{j}=$
$\displaystyle\frac{1}{2}(2-2\rho-\eta\lambda_{j}\pm\sqrt{\eta^{2}\lambda_{j}^{2}+4\rho^{2}})$
$\displaystyle=$ $\displaystyle 1-\rho\pm
i(\frac{1}{2}\eta\sigma_{j})\pm\frac{1}{2}\sqrt{4\rho^{2}-\eta^{2}\sigma_{j}^{2}}$
(31)

Since, for complex number $c$, $|c|^{2}=c\bar{c}$, the magnitude of
eigenvalues $|\mu_{j}|^{2}$ are given by,

$\displaystyle|\mu_{j}|^{2}=\biggl{(}1-\rho+i\frac{1}{2}\eta\sigma_{j}\pm\frac{1}{2}\sqrt{4\rho^{2}-\eta^{2}\sigma_{j}^{2}}\biggr{)}\times$
$\displaystyle\qquad\quad\;\biggl{(}1-\rho-i\frac{1}{2}\eta\sigma_{j}\pm\frac{1}{2}\overline{\sqrt{4\rho^{2}-\eta^{2}\sigma_{j}^{2}}}\biggr{)}$
$\displaystyle=(1\\!-\\!\rho)^{2}\\!+\\!\frac{1}{4}\eta^{2}\sigma_{j}^{2}\\!+\\!\frac{1}{4}|4\rho^{2}\\!-\\!\eta^{2}\sigma_{j}^{2}\\!\pm\\!(1\\!-\\!\rho)\Re(\sqrt{4\rho^{2}\\!-\\!\eta^{2}\sigma_{j}^{2}})$
$\displaystyle\quad\pm\frac{i}{4}\eta\sigma_{j}\overline{\sqrt{4\rho^{2}-\eta^{2}\sigma_{j}^{2}}}\mp\frac{i}{4}\eta\sigma_{j}\sqrt{4\rho^{2}-\eta^{2}\sigma_{j}^{2}}$
$\displaystyle=\begin{cases}1\\!-\\!2\rho\\!+\\!2\rho^{2}\\!\pm\\!(1\\!-\\!\rho)\sqrt{4\rho^{2}\\!-\\!\eta^{2}\sigma_{j}^{2}},\;\text{if
$4\rho^{2}\\!-\\!\eta^{2}\sigma_{j}^{2}\\!\geq\\!0$}\\\
1\\!-\\!2\rho\\!+\\!\frac{1}{2}\eta^{2}\sigma_{j}^{2}\\!\pm\\!\frac{1}{2}\eta\sigma_{j}\sqrt{\eta^{2}\sigma_{j}^{2}\\!-\\!4\rho^{2}},\;\text{if
$\eta^{2}\sigma_{j}^{2}\\!-\\!4\rho^{2}\\!\geq\\!0$}\end{cases}$

Suppose that for all $j\in[m]$, we choose
$0<\eta\leq\frac{2\rho}{\sigma_{\max}}\leq\frac{2\rho}{\sigma_{j}}$ and
$\rho>0$, which implies $4\rho^{2}-\eta^{2}\sigma_{j}^{2}\geq 0$. From (-B),
$\forall j\in[m]$ we have,

$\displaystyle|\mu_{j}|^{2}$
$\displaystyle=1-2\rho+2\rho^{2}\pm(1-\rho)\sqrt{4\rho^{2}-\eta^{2}\sigma_{j}^{2}}$
$\displaystyle\leq
1-2\rho+2\rho^{2}+(1-\rho)\sqrt{4\rho^{2}-\eta^{2}\sigma_{j}^{2}}$
$\displaystyle\leq
1-2\rho+2\rho^{2}+(1-\rho)\sqrt{4\rho^{2}-\eta^{2}\sigma_{\min}^{2}}$
$\displaystyle<1-2\rho+2\rho^{2}+(1-\rho)\sqrt{4\rho^{2}}$
$\displaystyle=1\;\;.$

According to classical linear system theory, e.g. [37, Theorem 8.3], the above
spectral radius analysis of the linear system (29) results in the following
linear convergence rate estimate:

$\displaystyle
V_{k}\leq\mathcal{O}\left(\\!\left(1-2\rho+2\rho^{2}+(1-\rho)\sqrt{4\rho^{2}-\eta^{2}\sigma_{\min}^{2}}\right)^{k}\\!\right)V_{0},$

where
$V_{k}:=\|x_{k}-x^{*}\|^{2}+\|y_{k}-y^{*}\|^{2}+\|\hat{x}_{k}-\hat{x}^{*}\|^{2}+\|\hat{y}_{k}-\hat{y}^{*}\|^{2}$.

Furthermore, we want to select the optimal step size $\rho,\eta$. The
immediate step is to substitute the optimal
$\eta=\frac{2\rho}{\sigma_{\max}}$, which yields the following inequality:

$\displaystyle|\mu_{j}|^{2}$ $\displaystyle\leq
1\\!-\\!2\rho\\!+\\!2\rho^{2}\\!+\\!(1-\rho)\sqrt{4\rho^{2}-\frac{4\rho^{2}}{\sigma_{\max}^{2}}\sigma_{j}^{2}}\;\;,\forall
j\in[m].$ (32)

The spectral radius is therefore given by choosing $\sigma_{j}=\sigma_{\min}$
above, i.e.,

$\displaystyle\max_{j\in[m]}$
$\displaystyle{|\mu_{j}|^{2}}=2\rho^{2}-2\rho+1+(1-\rho)2\rho\sqrt{1-\frac{\sigma_{\min}^{2}}{\sigma_{\max}^{2}}}$
$\displaystyle=\rho^{2}\biggl{(}2\\!-\\!2\sqrt{1-\frac{\sigma_{\min}^{2}}{\sigma_{\max}^{2}}}\biggr{)}\\!-\\!\rho\biggl{(}2\\!-\\!2\sqrt{1-\frac{\sigma_{\min}^{2}}{\sigma_{\max}^{2}}}\biggr{)}\\!+\\!1$
$\displaystyle\leq
1-\frac{1}{2}(1-\sqrt{1-\frac{\sigma_{\min}^{2}}{\sigma_{\max}^{2}}})$
$\displaystyle=\frac{1}{2}+\frac{1}{2}\sqrt{1-\frac{\sigma_{\min}^{2}}{\sigma_{\max}^{2}}}\;\;,\forall
j\in[m]$

where the last inequality comes from selecting optimal $\rho=\frac{1}{2}$ of a
quadratic polynomial of $\rho$. Using the fact that $\sqrt{1-x}\leq 1-x/2$, we
have

$\displaystyle\max_{j\in[m]}{|\mu_{j}|^{2}}\leq
1-\frac{1}{4}\frac{\sigma_{\min}^{2}}{\sigma_{\max}^{2}}$ (33)

Again, this results in the following linear convergence rate estimate:

$\displaystyle\textstyle
V_{k}\leq\mathcal{O}\left(\left(1-\frac{1}{4\kappa}\right)^{k}\right)\,V_{0}.$
(34)

###### Remark 1

We could also choose $\eta=\frac{2\rho}{\sigma_{\min}}$ such that
$\eta^{2}\sigma_{j}^{2}-4\rho^{2}\geq 0$. And we could construct a similar
linear convergence rate by repeating the above process. However, in practice,
we found that the step sizes $\eta=\frac{2\rho}{\sigma_{\max}},\rho=1/2$
always perform better in numerical experiments. Therefore, we choose this pair
of step sizes by default.

###### Remark 2

Since GDA method could be interpreted as a special case of DGDA method when
selecting $\rho=0$, the above step proves that when $\eta>0$, the GDA method
diverges for a bilinear objective function. Specifically, when
$\rho=0,\eta>0$, we have $\eta^{2}\sigma_{j}^{2}-4\rho^{2}>0$ and

$\displaystyle|\mu_{j}|^{2}=1\pm\frac{1}{2}\eta^{2}\sigma^{2}_{j},$ (35)

###  -C Proof of Theorem 3

The proof relies on the application of dissipativity theory to construct
Lyapunov functions and establish linear convergence. For more detailed
information, refer to [16].

According to [16], a linear dynamical system of the form:

$\displaystyle{\xi}_{k+1}=A{\xi}_{k}+Bw_{k}$ (36)

Here, ${\xi}\in\mathbb{R}^{n_{\xi}}$ is the state,
$w_{k}\in\mathbb{R}^{n_{w}}$ is the input, $A$ is the state transition matrix
and $B$ is the input matrix. Suppose that there exist a (Lyapunov) function
$V$, satisfying $V({\xi})\geq 0,\forall{\xi}\in\mathbb{R}^{n_{\xi}}$, some
$0\leq\alpha<1$ and a supply rate function $S({\xi}_{k},w_{k})\leq 0,\forall
k$ such that

$\displaystyle V({\xi}_{k+1})-\alpha^{2}V({\xi}_{k})\leq S({\xi}_{k},w_{k}).$
(37)

This dissipation inequality (37) implies that
$V({\xi}_{k+1})\leq\alpha^{2}V({\xi}_{k})$, and the state will approach a
minimum value ate equilibrium no slower than the linear rate $\alpha^{2}$. The
flowing theorem states how to construct the dissipation inequality (37) by
solving a semidefinite programming problem.

###### Theorem 5

[16][Theorem 2] Consider the following quadratic supply rate with
$X\in\mathbb{R}^{(n_{\xi}+n_{w})\times(n_{\xi}+n_{w})}$ and $X^{T}=X$

$\displaystyle S({\xi},w):=\begin{bmatrix}{\xi}\\\
w\end{bmatrix}^{T}X\begin{bmatrix}{\xi}\\\ w\end{bmatrix}.$ (38)

If there exists matrix $P\in\mathbb{R}^{n_{\xi}\times n_{\xi}}$ with $P\succeq
0$ such that

$\displaystyle\begin{bmatrix}A^{T}PA-\alpha^{2}P&A^{T}PB\\\
B^{T}PA&B^{T}PB\end{bmatrix}-X\leq 0,$ (39)

then the dissipation inequality holds for all trajectories of (36) with
$V({\xi})={\xi}^{T}P{\xi}$.

A major benefit of the proposed constructive dissipation approach is that it
replaces the trouble some component of a dynamical system (e.g. the gradient
term $w=\nabla_{{\xi}}f({\xi})$) by a quadratic constraint on its inputs and
outputs that is always satisfied, namely the supply rate constraint
$S(\xi,w)\leq 0$. This leads to a two-step novel approach to the convergence
analysis of optimization algorithms.

  1. 1.

Choose a proper quadratic supply rate function $S$ such that
$S({\xi}_{k},w_{k})\leq 0,\forall k$, that depends on the specific nonlinear
term.

  2. 2.

Solve the Linear Matrix Inequality (39) to obtain a storage function $V$ and
finding the linear convergence rate $\alpha$.

We will apply this methodology to analyze the DGDA update (12). Let
$z=[x,y]^{T},\hat{z}=[\hat{x},\hat{y}]^{T}$ and
$F(z_{k})=(\nabla_{x}f(x_{k},y_{k});-\nabla_{y}f(x_{k},y_{k}))$, and rewrite
(12) as in the form of (36):

$\displaystyle\begin{bmatrix}z_{k+1}\\\ \hat{z}_{k+1}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}z_{k}-\eta F(z_{k})-\rho(z_{k}-\hat{z}_{k})\\\
\hat{z}_{k}-\rho(\hat{z}_{k}-z_{k})\end{bmatrix}$ (40)
$\displaystyle=\begin{bmatrix}1-\rho&\rho\\\
\rho&1-\rho\end{bmatrix}\begin{bmatrix}z_{k}\\\
\hat{z}_{k}\end{bmatrix}+\begin{bmatrix}-\eta\\\ 0\end{bmatrix}w_{k}$

where $w_{k}=F(z_{k})$.

According to the previous discussion, the first step would be to choose a
proper quadratic supply rate function $S$ such that $S({\xi}_{k},w_{k})\leq
0,\forall k$, where ${\xi}_{k}=(z_{k};\hat{z}_{k})$ that depends on the
specific nonlinear term $w_{k}=F(z_{k})$. According to the equations (7) in
the work by Hu et al. (2017) [16] and Lemma 6 from the research by Lessard et
al. (2016) [18], the following applies to the nonlinear operator $F(z_{k})$
that meets the conditions specified in Assumption 1:

$\displaystyle S(z_{k},w_{k})=\begin{bmatrix}z_{k}\\\
w_{k}\end{bmatrix}^{T}\begin{bmatrix}2\mu LI&(-\mu+L)I\\\
(-\mu+L)I&2I\end{bmatrix}\begin{bmatrix}z_{k}\\\ w_{k}\end{bmatrix}\leq 0$
(41)

The conditions in Assumption 1 are also commonly referred to as being L-smooth
and m-strongly monotone, as can be found in related literature on variational
inequality problems, such as the works by [19, 6]]. Therefore, we could easily
extend the above LMI into the following supply rate function for DGDA updates,
by augmenting the states ${\xi}_{k}=(z_{k};\hat{z}_{k})$:

$\displaystyle S(\xi_{k},w_{k})=$ (42) $\displaystyle\begin{bmatrix}z_{k}\\\
\hat{z}_{k}\\\ w_{k}\end{bmatrix}^{T}\begin{bmatrix}2\mu LI&0&(-\mu+L)I\\\
0&0&0\\\ (-\mu+L)I&0&2I\end{bmatrix}\begin{bmatrix}z_{k}\\\ \hat{z}_{k}\\\
w_{k}\end{bmatrix}\leq 0$

as a proper quadratic supply rate function $S({\xi}_{k},w_{k})\leq 0$,
whenever $w_{k}=F(z_{k})$.

Finally, according to Theorem 5 and the above discussion, proving linear
convergence reduces to finding a positive definite a matrix
$P\in\mathbf{R}^{2(n+m)\times 2(n+m)}$, $\alpha\in[0,1)$ such that (39) is
satisfied, where the problem parameters are given by

$\displaystyle A=\begin{bmatrix}1-\rho&\rho\\\ \rho&1-\rho\end{bmatrix}\otimes
I,B=\begin{bmatrix}-\eta\\\ 0\end{bmatrix}\otimes I,$ (43) $\displaystyle
X=\begin{bmatrix}2\mu L&0&(-\mu+L)\\\ 0&0&0\\\
(-\mu+L)&0&2\end{bmatrix}\otimes I,$

where $\otimes$ is the Kronecker product. Due to the Kronecker structure of
this problem, this is equivalent to solving an LMI problem of dimension 3 by
3, with design parameters $P=\bar{P}\otimes I$, with
$\bar{P}\in\mathbf{R}^{2\times 2}$, $\alpha^{2}\in[0,1)$, $\rho$ and $\eta$.
Because this Linear Matrix Inequality is simple (3 by 3), it can be solved
using analytical methods. This, in turn, results in a feasible solution for
the LMI, denoted as follows:

$\displaystyle\rho=\frac{1}{2},\qquad\eta=\frac{1}{L+\mu},$ (44)
$\displaystyle\alpha^{2}=\frac{3L^{2}+2L\mu+3\mu^{2}+\sqrt{(L+\mu)^{4}+16L^{2}\mu^{2}}}{4(L+\mu)^{2}},$
$\displaystyle P=\begin{bmatrix}(L+\mu)^{2}&0\\\
0&(L+\mu)^{2}\end{bmatrix}\otimes I.$

After substituting the definition for condition number $\kappa:=L/\mu$, the
convergence rate $\alpha^{2}$ simplifies to:

$\displaystyle\alpha^{2}=1-\kappa^{-1}+\mathcal{O}\bigl{(}(\frac{\mu}{L})^{2}\bigr{)}$
(45)

###  -D Proof of Corollary 4

According to Theorem 3, the linear convergence rate estimate of DGDA is

$\displaystyle\frac{3L^{2}+2L\mu+3\mu^{2}+\sqrt{(L+\mu)^{4}+16L^{2}\mu^{2}}}{4(L+\mu)^{2}}$
(46)
$\displaystyle=\frac{3\kappa^{2}+2\kappa+3+\sqrt{(\kappa+1)^{4}+16\kappa^{2}}}{4(\kappa+1)^{2}},$
(47)

where $\kappa=L/\mu$.

According to Theorem $6\&7$ [2] and Theorem $4\&7$ [1], the standard known
linear convergence rate estimate of EG and OGDA is

$\displaystyle 1-\frac{\mu}{4L}=1-\frac{1}{4\kappa}.$ (48)

By simple algebraic calculation, it can be shown that as a function of
$\kappa$, when $\kappa\geq 2$, the following polynomial is always nonnegative,
i.e.,

$\displaystyle
1-\frac{1}{4\kappa}-\frac{3\kappa^{2}+2\kappa+3+\sqrt{(\kappa+1)^{4}+16\kappa^{2}}}{4(\kappa+1)^{2}}\geq
0.$ (49)


functions are off-center, but also because they cannot fit inside one unit
cell, and this phenomenon cannot be captured by the usual Zak phase. This
situation is reminiscent to the results of a very recent work [114], in which
the authors propose a 3D topological material in which the topology is also
not captured by the standard tenfold classification. The reason for this is
that the maximally localized Wannier functions of the system cannot fit into a
single primitive unit cell, and this causes the apparition of topological
surface states. We believe the trimer chain and the delta CSSH ladder might
be, in this sense, 1D analogues of this phenomenon, termed by the authors
“multicellular topology”, although further research is needed.

##  Appendix E Edge state protection

The topological edge states in the second gap of the CSSH ladders, when
present, are protected by the chiral and PH symmetries (or only by the PH
symmetry, in the D class $\Delta$-CSSH ladder). These symmetries force the
spectrum to be symmetric around zero, and will typically relate two linear
combinations of the left and right states. These linear combinations will be
energy eigenstates, and then they will typically have small, finite energies
around zero, as long as hybridization of the left and right modes is allowed.
However, if hybridization is inhibited, whether because the system is large
enough or due to another mechanism (like AB caging), the left and right edge
modes will also be energy eigenstates. The protecting symmetries will then
force their energy to be exactly at zero. This is the sign of topological
protection in a 1D system.

In Fig. 10, the effect of hopping disorder in the terms $t_{i}$ of the CSSH
Hamiltonian (1) is studied. The disorder is modeled as:

$t_{i}=t_{i}^{c}+\delta JR,$ (65)

where $\delta J$ is the level of disorder, $t_{i}^{c}$ is the value of the
amplitude in a clean system, and $R$ is a uniformly drawn random number
between $-1/2$ and $1/2$, chosen individually for each term in the Hamiltonian
and for each value of $\delta J$. The level of disorder ranges from zero to
$\xi$ (with $\xi=1$) to better appreciate its effects, but it must be noted
that the level of disorder expected in a real experimental setup would be much
smaller, on the order of $0.1\xi$ [2, 3, 6].

In Fig. 10 (a) and (e), protected zero modes in the
$\mathbin{\rotatebox[origin={c}]{90.0}{\scalebox{1.5}{$\Join$}}}$\- and
$\Delta$-CSSH ladders are represented. Their energy is fixed at exactly zero,
even in the presence of intense disorder, because disorder in the hopping
amplitudes of the system does not break the protecting symmetries.

If both the chiral and the PH symmetries are absent, or if we consider edge
states that are not located in the gap around zero energy, then the situation
changes. The absence of each of these conditions must be studied separately.

If both protecting symmetries are absent, then the system belongs to a
symmetry class which cannot harbor topological phases in 1D, according to the
tenfold classification. However, the presence of crystalline symmetries can
expand and modify this classification, allowing for the presence of degenerate
topological end modes. The system is then a topological crystalline insulator.
As described in the text, this is the case in the Creutz and
$\mathbin{\rotatebox[origin={c}]{90.0}{\scalebox{1.5}{$\Join$}}}$-CSSH ladders
with $\phi\neq 0,\pi$ and $m\neq 0$: the central topological end modes of the
spectrum are degenerate due to inversion symmetry, at an energy value
different than zero. However, inversion symmetry will be broken by any
physically relevant disorder present in the system, and so we refrain from
categorizing these states as “protected” in this sense. States of this kind
are depicted in Fig. 10 (c). Their energy stays close to its value in a clean
system, but it starts drifting up and down as the disorder strength increases.

Figure 10: Spectra of hourglass (a-d) and delta (e-f) CSSH ladders with open
boundary conditions as a function of the level of hopping amplitude disorder
$\delta J$, for the parameters shown above each plot. All simulations use
$L=16$, $\xi=1$. The parameters were chosen so that the studied end mode in
each case has a distance of $\Delta E=0.35$ with the closest band, for the
effects of disorder to be comparable. The end states in the second (first or
third) gap are represented in red (green). (a) Protected zero modes. Their
energy is fixed, even for large levels of disorder. (b) Non-protected edge
states. For large levels of disorder, their energies drift up and down. (c,d)
Edge states in a topological crystalline insulator, not protected against
disorder. (e) Protected zero mode in the $\Delta$-CSSH ladder. (f) Non-
protected left edge states $\ket{\ell\pm}$, which do not follow the usual
bulk-boundary correspondence.

The second kind of non-protected topological states are those in the first and
third gaps of the system. The protecting symmetries can only cause their
energies to have opposite values, not to be fixed at a given value. For this
reason, they are not prevented from drifting into the bulk bands when a high
level of disorder is present.

In Fig. 10 (b), the response of non-protected states in the first and third
gaps of the
$\mathbin{\rotatebox[origin={c}]{90.0}{\scalebox{1.5}{$\Join$}}}$-CSSH ladder
spectrum can be seen for different values of $\delta J$. As disorder
increases, their energies start fluctuating around their value in the clean
system. This also happens in the topological crystalline insulator case, Fig.
10 (d), and for the left edge modes in the $\Delta$-CSSH ladder which do not
follow the usual bulk-boundary correspondence, Fig. 10 (f). In these two last
cases, the fluctuations are more pronounced than in Fig. 10 (b), probably due
to the absence of PH and chiral symmetries in (d), and the different nature of
the $\ket{\ell\pm}$ states in (f). In all the analyzed models, except in the
topological crystalline insulator phase, the spectrum is still symmetric
around zero energy, because the PH and chiral symmetries are not broken. This
means that when the state in the first gap drift upwards, the one in the third
gap will drift downwards, and vice versa.

##  Appendix F Rhombus chain implementation

We consider a generalized rhombus chain Hamiltonian with a 6-site unit cell
(see Fig. 8). We label the upper states $A$ and $C$, the lower states are $B$
and $D$, and the middle states are $X$ and $Z$:

$\displaystyle\mathcal{H}_{R}=-\sum_{n}\left\\{\vphantom{\frac{\epsilon}{2}}v_{1}^{+}\ket{n,A}\bra{n,X}+v_{2}^{+}\ket{n,X}\bra{n,B}\right.$
$\displaystyle+v_{3}^{+}\ket{n,Z}\bra{n,A}+v_{4}^{+}\ket{n,B}\bra{n,Z}$
$\displaystyle+w_{1}^{+}\ket{n,C}\bra{n,Z}+w_{2}^{+}\ket{n,Z}\bra{n,D}$
$\displaystyle+w_{3}^{+}\ket{n+1,X}\bra{n,C}+w_{4}^{+}\ket{n,D}\bra{n+1,X}$
$\displaystyle\left.+\frac{\epsilon}{2}\ket{n,X}\bra{n,X}+\frac{\epsilon}{2}\ket{n,Z}\bra{n,Z}+h.c.\right\\},$
(66)

where $v_{i}^{+}=v_{i}e^{i\phi/4}$ and $w_{i}^{+}=w_{i}e^{i\phi/4}$. The
phases implement a magnetic flux of $\phi$ through each rhombus of the chain,
and can be generated using Floquet driving [27].

If we consider the regime where $\epsilon\gg v_{i},w_{i}$, the Hilbert space
of the system will decouple into a low-energy subspace, with zero occupation
in the middle sites, and a high-energy subspace. As indicated in [27], the
low-energy subspace has the geometry of a Creutz ladder. Hence, changing the
tunneling amplitudes in certain ways will let us obtain the different CSSH
ladders. In particular, all topological flat band ladders and the SSH-like
topological zero modes of the
$\mathbin{\rotatebox[origin={c}]{90.0}{\scalebox{1.5}{$\Join$}}}$-CSSH ladder
can be implemented in a straight-forward way. This can allow experimentalists
to explore the various topological zero modes of Section 4 and Appendix C, and
the flat-band dynamics of the models. The hopping amplitudes can be tuned by
simply changing the distance between waveguides.

The low-energy effective Hamiltonian, obtained using a Schrieffer-Wolff
transformation [115], is:

$\displaystyle\mathcal{H}_{\mathrm{eff}}=\mathcal{H}_{C\\!S\\!S\\!H}+\sum_{n,\alpha}\epsilon_{\alpha}\ket{n,\alpha}\bra{n,\alpha},$
(67)

where $\alpha=A,B,C,D$ and $\mathcal{H}_{C\\!S\\!S\\!H}$ is the CSSH
Hamiltonian of equation (1), with the following parameters:

$\displaystyle J_{1}^{+}=\frac{v_{3}w_{1}}{\epsilon}e^{i\phi/2};\hskip
20.0ptJ_{2}^{+}=\frac{w_{2}v_{4}}{\epsilon}e^{i\phi/2}$ $\displaystyle
J_{3}^{+}=\frac{w_{3}v_{1}}{\epsilon}e^{i\phi/2};\hskip
20.0ptJ_{4}^{+}=\frac{v_{2}w_{4}}{\epsilon}e^{i\phi/2}$ $\displaystyle
t_{1}=\frac{v_{3}w_{2}}{\epsilon};\hskip
20.0ptt_{2}=\frac{w_{1}v_{4}}{\epsilon}$ $\displaystyle
t_{3}=\frac{w_{3}v_{2}}{\epsilon};\hskip
20.0ptt_{4}=\frac{v_{1}w_{4}}{\epsilon}$ $\displaystyle
m_{1}=\frac{1}{\epsilon}(v_{2}v_{1}e^{i\phi/2}+v_{3}v_{4}e^{-i\phi/2})$
$\displaystyle
m_{2}=\frac{1}{\epsilon}(w_{2}w_{1}e^{i\phi/2}+w_{3}w_{4}e^{-i\phi/2})$
$\displaystyle\epsilon_{A}=\frac{v_{1}^{2}+v_{3}^{2}}{\epsilon};\hskip
20.0pt\epsilon_{B}=\frac{v_{2}^{2}+v_{4}^{2}}{\epsilon}$
$\displaystyle\epsilon_{C}=\frac{w_{1}^{2}+w_{3}^{2}}{\epsilon};\hskip
20.0pt\epsilon_{D}=\frac{w_{2}^{2}+w_{4}^{2}}{\epsilon}.$

In the most general case, the vertical amplitudes are not real.

The hourglass CSSH ladder can be obtained by choosing the parameters:

$\displaystyle v_{1}=v_{2}=w_{3}=w_{4}=s^{\prime}$ $\displaystyle
v_{3}=v_{4}=w_{1}=w_{2}=s,$

which generates the amplitudes:

$\displaystyle J=\frac{s^{2}}{\epsilon};\hskip
20.0ptJ^{\prime}=\frac{s^{\prime 2}}{\epsilon}$ $\displaystyle
m_{1}=\frac{1}{\epsilon}(s^{\prime
2}e^{i\phi/2}+s^{2}e^{-i\phi/2})=m_{2}^{*},$

with an additional on-site energy of $(s^{2}+s^{\prime 2})/\epsilon$ on all
sites. This model is represented in Fig. 8 (b), with long (short) links to
represent weak (strong) tunneling amplitudes. If $\phi=\pi$, $m_{1}=m_{2}=0$
and the flat-band
$\mathbin{\rotatebox[origin={c}]{90.0}{\scalebox{1.5}{$\Join$}}}$-CSSH ladder
is obtained. If $\phi=0$, $m_{1}=m_{2}=(s^{2}+s^{\prime 2})/\epsilon$ and the
ladder has SSH-like topological zero modes if $s^{\prime}<s$. For other values
of $\phi$, the vertical links have the same magnitude and opposite phase. They
are always non-zero, so the rungless, non-$\pi$ flux topological case cannot
be obtained in this way.

The rhomboid CSSH ladder Hamiltonian is obtained from (67) if we choose:

$\displaystyle v_{2}=v_{3}=s^{\prime};\hskip 20.0ptv_{1}=v_{4}=s;\hskip
20.0ptw_{i}=w\,\,\,\forall i,$

with the parameters:

$\displaystyle J=\frac{sw}{\epsilon};\hskip
20.0ptJ^{\prime}=\frac{s^{\prime}w}{\epsilon}$ $\displaystyle
m_{1}=\frac{2ss^{\prime}}{\epsilon}\cos{\phi/2};\hskip
10.0ptm_{2}=\frac{2w^{2}}{\epsilon}\cos{\phi/2}$
$\displaystyle\epsilon_{A}=\epsilon_{B}=\frac{s^{2}+s^{\prime
2}}{\epsilon};\hskip 20.0pt\epsilon_{C}=\epsilon_{D}=\frac{2w^{2}}{\epsilon}$

This lattice is depicted in Fig. 8 (c). In the $\pi$-flux case,
$m_{1}=m_{2}=0$, and we can choose $w=\sqrt{s^{2}+s^{\prime 2}}$ to keep all
on-site energies $\epsilon_{\alpha}$ equal and obtain the topological flat-
band -CSSH ladder. For other values of $\phi$, $w=\sqrt{ss^{\prime}}$ might be
chosen to get $m_{1}=m_{2}$, and then the on-site energies $\epsilon_{\alpha}$
could be compensated using the appropriate chemical potentials on the rhombus
chain sites. Like in the hourglass case, the rungless non-$\pi$-flux case
cannot be achieved in this way.

Finally, to implement the delta CSSH ladder, we set:

$\displaystyle v_{1}=v_{3}=s^{\prime};\hskip 15.0ptv_{2}=v_{4}=s;\hskip
15.0ptw_{i}=\sqrt{ss^{\prime}}\,\,\,\forall i,$

and obtain [see Fig. 8 (d)]:

$\displaystyle J=\frac{s^{2}s^{\prime}}{\epsilon};\hskip
20.0ptJ^{\prime}=\frac{s^{\prime 2}s}{\epsilon}$ $\displaystyle
m_{1}=m_{2}=m=\frac{2ss^{\prime}}{\epsilon}\cos{\phi/2};\hskip
20.0pt\epsilon_{A}=\frac{2s^{\prime 2}}{\epsilon}$
$\displaystyle\epsilon_{B}=\frac{2s^{2}}{\epsilon};\hskip
20.0pt\epsilon_{C}=\epsilon_{D}=\frac{2ss^{\prime}}{\epsilon}.$

In this model, the energies $\epsilon_{\alpha}$ cannot be made equal by fixing
$w$, so they must be compensated using on-site potentials. If $\phi=\pi$, then
$m=0$ and the topological flat-band case is obtained. Like in the other cases,
the $\phi\neq\pi$ regime is not a strong topological phase, because $m$ is
finite there. The $\pi$-flux case with $m\neq 0$, which belongs to the D
class, cannot be obtained directly using this method. To implement it, an
additional tunneling amplitude of value $m$ could be established between the
top and bottom sites of each rhombi, to create the vertical links of the
ladder. This makes the lattice similar to the cold atom proposal in [12], and
both of them could be used to implement the rungless topological regime for
any value of $\phi$, and also the $\pi$-flux with $m\neq 0$ case. In the
$\Delta$-CSSH ladder, the latter corresponds to the D-class TI.

## References

  * Altland and Zirnbauer [1997] A. Altland and M. R. Zirnbauer, Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures, Physical Review B 55, 1142 (1997).
  * Lang and Büchler [2017] N. Lang and H. P. Büchler, Topological networks for quantum communication between distant qubits, npj Quantum Information 3, 47 (2017).
  * Mei  _et al._ [2018] F. Mei, G. Chen, L. Tian, S. L. Zhu, and S. Jia, Robust quantum state transfer via topological edge states in superconducting qubit chains, Physical Review A 98, 012331 (2018).
  * Qi  _et al._ [2020] L. Qi, G. L. Wang, S. Liu, S. Zhang, and H. F. Wang, Engineering the topological state transfer and topological beam splitter in an even-sized Su-Schrieffer-Heeger chain, Physical Review A 102, 022404 (2020).
  * Estarellas  _et al._ [2017] M. P. Estarellas, I. D’Amico, and T. P. Spiller, Topologically protected localised states in spin chains, Scientific Reports 7, 42904 (2017).
  * Boross  _et al._ [2019] P. Boross, J. K. Asbóth, G. Széchenyi, L. Oroszlány, and A. Pályi, Poor man’s topological quantum gate based on the Su-Schrieffer-Heeger model, Physical Review B 100, 045414 (2019).
  * Tovmasyan  _et al._ [2013] M. Tovmasyan, E. P. Van Nieuwenburg, and S. D. Huber, Geometry-induced pair condensation, Physical Review B 88, 220510(R) (2013).
  * Li and Chen [2015a] L. Li and S. Chen, Characterization of topological phase transitions via topological properties of transition points, Physical Review B 92, 085118 (2015a).
  * Padavić  _et al._ [2018] K. Padavić, S. S. Hegde, W. DeGottardi, and S. Vishveshwara, Topological phases, edge modes, and the Hofstadter butterfly in coupled Su-Schrieffer-Heeger systems, Physical Review B 98, 024205 (2018).
  * Velasco and Paredes [2019] C. G. Velasco and B. Paredes, Classification of topological ladder models (2019), arXiv:1907.11460 .
  * Piga  _et al._ [2017] A. Piga, M. Lewenstein, S.-J. Ran, J. Jünemann, A. Bermúdez, and M. Rizzi, Exploring Interacting Topological Insulators with Ultracold Atoms: The Synthetic Creutz-Hubbard Model, Physical Review X 7, 031057 (2017).
  * Sun and Lim [2017] N. Sun and L.-K. K. Lim, Quantum charge pumps with topological phases in a Creutz ladder, Physical Review B 96, 035139 (2017).
  * Kuno [2020] Y. Kuno, Extended flat band, entanglement, and topological properties in a Creutz ladder, Physical Review B 101, 184112 (2020).
  * van Dalum  _et al._ [2021] G. A. R. van Dalum, C. Ortix, and L. Fritz, Magnetic impurities along the edge of a quantum spin Hall insulator: Realizing a one-dimensional AIII insulator, Physical Review B 103, 075115 (2021).
  * Liu and Hersam [2019] X. Liu and M. C. Hersam, 2D materials for quantum information science, Nature Reviews Materials 4, 669-684 (2019).
  * Gholizadeh  _et al._ [2018] S. Gholizadeh, M. Yahyavi, and B. Hetényi, Extended Creutz ladder with spin-orbit coupling: A one-dimensional analog of the Kane-Mele model, Epl 122, 27001 (2018).
  * Creutz [1999] M. Creutz, End States, Ladder Compounds, and Domain-Wall Fermions, Physical Review Letters 83, 2636-2639 (1999).
  * Vidal  _et al._ [2000] J. Vidal, B. Douçot, R. Mosseri, and P. Butaud, Interaction induced delocalization for two particles in a periodic potential, Phys. Rev. Lett. 85, 3906 (2000).
  * Douçot and Vidal [2002] B. Douçot and J. Vidal, Pairing of Cooper pairs in a fully frustrated Josephson-junction chain, Phys. Rev. Lett. 88, 227005 (2002).
  * Creffield and Platero [2010] C. E. Creffield and G. Platero, Coherent control of interacting particles using dynamical and Aharonov-Bohm phases, Physical Review Letters 105, 086804 (2010).
  * Pelegrí  _et al._ [2019] G. Pelegrí, A. M. Marques, R. G. Dias, A. J. Daley, J. Mompart, and V. Ahufinger, Topological edge states and Aharanov-Bohm caging with ultracold atoms carrying orbital angular momentum, Phys. Rev. A 99, 023613 (2019).
  * Vidal  _et al._ [1998] J. Vidal, R. Mosseri, and B. Douçot, Aharonov-Bohm cages in two-dimensional structures, Phys. Rev. Lett. 81, 5888 (1998).
  * Vidal  _et al._ [2001] J. Vidal, P. Butaud, B. Douçot, and R. Mosseri, Disorder and interactions in Aharonov-Bohm cages, Phys. Rev. B 64, 155306 (2001).
  * Abilio  _et al._ [1999] C. C. Abilio, P. Butaud, T. Fournier, B. Pannetier, J. Vidal, S. Tedesco, and B. Dalzotto, Magnetic field induced localization in a two-dimensional superconducting wire network, Physical Review Letters 83, 5102-5105 (1999).
  * Naud  _et al._ [2001] C. Naud, G. Faini, and D. Mailly, Aharonov-Bohm cages in 2D normal metal networks, Physical Review Letters 86, 5104-5107 (2001).
  * Shinohara  _et al._ [2002] K. Shinohara, T. Aoki, and A. Morinaga, Scalar Aharonov-Bohm effect for ultracold atoms, Physical Review A 66, 042106 (2002).
  * Mukherjee  _et al._ [2018] S. Mukherjee, M. Di Liberto, P. Öhberg, R. R. Thomson, and N. Goldman, Experimental Observation of Aharonov-Bohm Cages in Photonic Lattices, Physical Review Letters 121, 075502 (2018).
  * Kremer  _et al._ [2020] M. Kremer, I. Petrides, E. Meyer, M. Heinrich, O. Zilberberg, and A. Szameit, A square-root topological insulator with non-quantized indices realized with photonic Aharonov-Bohm cages, Nature Communications 11, 907 (2020).
  * Jörg  _et al._ [2020] C. Jörg, G. Queraltó, M. Kremer, G. Pelegrí, J. Schulz, A. Szameit, G. von Freymann, J. Mompart, and V. Ahufinger, Artificial gauge field switching using orbital angular momentum modes in optical waveguides, Light: Science & Applications 9, 150 (2020).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43-50 (2018).
  * Takayoshi  _et al._ [2013] S. Takayoshi, H. Katsura, N. Watanabe, and H. Aoki, Phase diagram and pair Tomonaga-Luttinger liquid in a Bose-Hubbard model with flat bands, Physical Review A 88, 063613 (2013).
  * Mondaini  _et al._ [2018] R. Mondaini, G. G. Batrouni, and B. Grémaud, Pairing and superconductivity in the flat band: Creutz lattice, Physical Review B 98, 155142 (2018).
  * Zurita  _et al._ [2019] J. Zurita, C. E. Creffield, and G. Platero, Topology and interactions in the photonic Creutz and Creutz-Hubbard ladders, Advanced Quantum Technologies 3, 1900105 (2019).
  * Hügel and Paredes [2014] D. Hügel and B. Paredes, Chiral ladders and the edges of quantum Hall insulators, Physical Review A 89, 023619 (2014).
  * Zhou and Du [2019] L. Zhou and Q. Du, Floquet topological phases with fourfold-degenerate edge modes in a driven spin-1/2 Creutz ladder, Physical Review A 101, 033607 (2019).
  * Molignini [2020] P. Molignini, Edge mode manipulation through commensurate multifrequency driving, Physical Review B 102, 235143 (2020).
  * Liu  _et al._ [2019] Z. X. Liu, Z. H. Li, and A. M. Wang, Fractional charged edge states in ladder topological insulators, Journal of Physics Condensed Matter 31, 125402 (2019).
  * Sticlet  _et al._ [2014] D. Sticlet, L. Seabra, F. Pollmann, and J. Cayssol, From fractionally charged solitons to Majorana bound states in a one-dimensional interacting model, Physical Review B 89, 115430 (2014).
  * Budich and Ardonne [2013] J. C. Budich and E. Ardonne, Equivalent topological invariants for one-dimensional Majorana wires in symmetry class $D$, Phys. Rev. B 88, 075419 (2013).
  * Lin  _et al._ [2017] Y. Lin, W. Hao, M. Wang, J. Qian, and H. Guo, Topological superconductors from one-dimensional periodically modulated Majorana chains, Scientific Reports 7, 9210 (2017).
  * Martinez Alvarez and Coutinho-Filho [2019] V. M. Martinez Alvarez and M. D. Coutinho-Filho, Edge states in trimer lattices, Phys. Rev. A 99, 013833 (2019).
  * Asbóth  _et al._ [2015] J. K. Asbóth, L. Oroszlány, and A. Pályi,  _A Short Course on Topological Insulators_ (Springer, 2015).
  * Bello  _et al._ [2016] M. Bello, C. E. Creffield, and G. Platero, Long-range doublon transfer in a dimer chain induced by topology and ac fields, Scientific Reports 6, 22562 (2016).
  * Ivanov and Richter [1997] N. Ivanov and J. Richter, Competition between plaquette and dimer phases in Heisenberg chains, Physics Letters A 232, 308-312 (1997).
  * Richter  _et al._ [1998] J. Richter, N. B. Ivanov, and J. Schulenburg, The antiferromagnetic spin-1/2 chain with competing dimers and plaquettes: Numerical versus exact results, Journal of Physics Condensed Matter 10, 3635-3649 (1998).
  * Čanová  _et al._ [2004] L. Čanová, J. Strečka, and M. Jaščur, Exact results of the Ising-Heisenberg model on the diamond chain with spin-1/2, Czechoslovak Journal of Physics 54, 579-582 (2004).
  * Paulinelli  _et al._ [2013] H. G. Paulinelli, S. M. De Souza, and O. Rojas, Thermal entanglement in an orthogonal dimer-plaquette chain with alternating Ising-Heisenberg coupling, Journal of Physics Condensed Matter 25, 306003 (2013).
  * Verkholyak and Strečka [2013] T. Verkholyak and J. Strečka, Exact solution for a quantum spin-12 Ising-Heisenberg orthogonal-dimer chain with Heisenberg intradimer and Ising interdimer interactions, Physical Review B 88, 134419 (2013).
  * Leykam  _et al._ [2018] D. Leykam, A. Andreanov, and S. Flach, Artificial flat band systems: From lattice models to experiments, Advances in Physics: X 3, 1473052 (2018).
  * Nandy [2019] A. Nandy, Analytical study of quasi-one dimensional flat band networks and slow light analogue, Acta Physica Polonica A 136, 164-173 (2019).
  * Rojas  _et al._ [2019] O. Rojas, J. Strečka, M. L. Lyra, and S. M. De Souza, Universality and quasicritical exponents of one-dimensional models displaying a quasitransition at finite temperatures, Physical Review E 99, 042117 (2019).
  * Gálisová  _et al._ [2021] L. Gálisová, J. Strečka, T. Verkholyak, and S. Havadej, Magnetization plateaus and bipartite entanglement of an exactly solved spin-1/2 Ising-Heisenberg orthogonal-dimer chain, Physica E: Low-Dimensional Systems and Nanostructures 125, 114089 (2021).
  * Yan and Wan [2014] Z. Yan and S. Wan, Topological phases, topological flat bands, and topological excitations in a one-dimensional dimerized lattice with spin-orbit coupling, EPL 107, 47007 (2014).
  * Han  _et al._ [2014] S. Han, H. Jiang, X. Li, and S. Xie, Spontaneous spin polarization in organic thiophene oligomers, Organic Electronics 15, 240-244 (2014).
  * Jiang  _et al._ [2015] H. Jiang, C. Zhang, X. Hu, G. Hu, and S. Xie, Spin polarization of polaron in quasi-one dimensional organic system, Modern Physics Letters B 29, 1450266 (2015).
  * Bahari and Hosseini [2020] M. Bahari and M. V. Hosseini, Topological properties of a generalized spin–orbit-coupled Su–Schrieffer–Heeger model, Physica E: Low-Dimensional Systems and Nanostructures 119, 113973 (2020).
  * Volkova and Marinin [2017] L. M. Volkova and D. V. Marinin, Frustrated Antiferromagnetic Spin Chains of Edge-Sharing Tetrahedra in Volcanic Minerals K3Cu3(Fe0.82Al0.18)O2(SO4)4 and K4Cu4O2(SO4)4MeCl, Journal of Superconductivity and Novel Magnetism 30, 959-971 (2017).
  * Brenig and Becker [2001] W. Brenig and K. W. Becker, Magnetism of a tetrahedral cluster spin chain, Physical Review B 64, 214413 (2001).
  * Ootomo and Ono [2000] F. Ootomo and Y. Ono, Acoustic polaron in a ladder-type electron-lattice system, Journal of the Physical Society of Japan 69, 3320-3327 (2000).
  * Zhang and Zhou [2017] S. L. Zhang and Q. Zhou, Two-leg Su-Schrieffer-Heeger chain with glide reflection symmetry, Physical Review A 95, 061601 (2017).
  * Nersesyan [2020] A. A. Nersesyan, Phase diagram of an interacting staggered Su-Schrieffer-Heeger two-chain ladder close to a quantum critical point, Physical Review B 102, 045108 (2020).
  * Jangjan and Hosseini [2020] M. Jangjan and M. V. Hosseini, Floquet engineering of topological metal states and hybridization of edge states with bulk states in dimerized two-leg ladders, Scientific Reports 10, 14256 (2020).
  * Ryu  _et al._ [2020] J.-W. Ryu, S. Woo, N. Myoung, and H. C. Park, Topological edge states in bowtie ladders with different cutting edges (2020), arXiv:2012.00270 .
  * Drescher and Mielke [2017] M. Drescher and A. Mielke, Hard-core bosons in flat band systems above the critical density, European Physical Journal B 90, 217 (2017).
  * Schnyder  _et al._ [2008] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).
  * Kitaev [2009] A. Kitaev, Periodic table for topological insulators and superconductors, AIP Conference Proceedings 1134, 22-30 (2009).
  * Stanescu [2016] T. D. Stanescu,  _Introduction to Topological Quantum Matter and Quantum Computation_ (CRC Press, 2016).
  * Alaeian  _et al._ [2019] H. Alaeian, C. W. S. Chang, M. V. Moghaddam, C. M. Wilson, E. Solano, and E. Rico, Creating lattice gauge potentials in circuit QED: The bosonic Creutz ladder, Physical Review A 99, 053834 (2019).
  * Leonforte  _et al._ [2021] L. Leonforte, A. Carollo, and F. Ciccarello, Vacancy-like Dressed States in Topological Waveguide QED, Phys. Rev. Lett. 126, 063601 (2021).
  * Jackiw [2012] R. Jackiw, Fractional and Majorana fermions: the physics of zero-energy modes, Physica Scripta T146, 014005 (2012).
  * Herviou [2017] M. L. Herviou,  _Topological Phases and Majorana Fermions_ , Ph.D. thesis, Université Paris-Saclay (2017).
  * Pérez-González  _et al._ [2019a] B. Pérez-González, M. Bello, Á. Gómez-León, and G. Platero, Interplay between long-range hopping and disorder in topological systems, Physical Review B 99, 035146 (2019a).
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in Polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * Zak [1989] J. Zak, Berry’s phase for energy bands in solids, Phys. Rev. Lett. 62, 2747 (1989).
  * Ryu and Hatsugai [2002] S. Ryu and Y. Hatsugai, Topological Origin of Zero-Energy Edge States in Particle-Hole Symmetric Systems, Physical Review Letters 89, 077002 (2002).
  * Lau and Ortix [2018] A. Lau and C. Ortix, Novel topological insulators from crystalline symmetries, European Physical Journal: Special Topics 227, 1309 (2018).
  * Cariglia [2014] M. Cariglia, Hidden symmetries of dynamics in classical and quantum physics, Reviews of Modern Physics 86, 1283-1336 (2014).
  * Hou [2013] J. M. Hou, Hidden-Symmetry-Protected Topological Semimetals on a Square Lattice, Physical Review Letters 111, 130403 (2013).
  * Fukui  _et al._ [2013] T. Fukui, K. I. Imura, and Y. Hatsugai, Symmetry protected weak topological phases in a superlattice, Journal of the Physical Society of Japan 82, 073708 (2013).
  * Li and Chen [2015b] L. Li and S. Chen, Hidden-symmetry-protected topological phases on a one-dimensional lattice, EPL 109, 40006 (2015b).
  * Hou and Chen [2016] J. M. Hou and W. Chen, Hidden-symmetry-protected quantum pseudo-spin Hall effect in optical lattices, Physical Review A 93, 063626 (2016).
  * Hou and Chen [2017] J. M. Hou and W. Chen, Hidden symmetry-protected Z2 topological insulator in a cubic lattice, Physical Review B 96, 235108 (2017).
  * Xiao and An [2018] J. Xiao and J. An, Boundary states in the chiral symmetric systems with a spatial symmetry, Journal of Physics Condensed Matter 30, 055002 (2018).
  * Marques and Dias [2019] A. M. Marques and R. G. Dias, One-dimensional topological insulators with noncentered inversion symmetry axis, Physical Review B 100, 041104 (2019).
  * Wang  _et al._ [2018] L. Wang, W. Cai, M. Bie, X. Zhang, and J. Xu, Zak phase and topological plasmonic Tamm states in one-dimensional plasmonic crystals, Optics Express 26, 28963 (2018).
  * Wang  _et al._ [2019] H.-X. Wang, G.-Y. Guo, and J.-H. Jiang, Band topology in classical waves: Wilson-loop approach to topological numbers and fragile topology, New Journal of Physics 21, 093029 (2019).
  * Zhu  _et al._ [2018] W. Zhu, Y.-q. Ding, J. Ren, Y. Sun, Y. Li, H. Jiang, and H. Chen, Zak phase and band inversion in dimerized one-dimensional locally resonant metamaterials, Phys. Rev. B 97, 195307 (2018).
  * Deng  _et al._ [2013] S. Deng, G. Ortiz, and L. Viola, Multiband $s$-wave topological superconductors: Role of dimensionality and magnetic field response, Phys. Rev. B 87, 205414 (2013).
  * Li  _et al._ [2016] L. Li, C. Yang, and S. Chen, Topological invariants for phase transition points of one-dimensional Z2 topological systems, European Physical Journal B 89, 195 (2016).
  * Ke  _et al._ [2016] Y. Ke, X. Qin, F. Mei, H. Zhong, Y. S. Kivshar, and C. Lee, Topological phase transitions and Thouless pumping of light in photonic waveguide arrays, Laser & Photonics Reviews 10, 995 (2016).
  * Ganeshan  _et al._ [2013] S. Ganeshan, K. Sun, and S. Das Sarma, Topological zero-energy modes in gapless commensurate Aubry-André- Harper models, Physical Review Letters 110, 180403 (2013).
  * Kang  _et al._ [2018] J. H. Kang, J. H. Han, and Y. Shin, Realization of a Cross-Linked Chiral Ladder with Neutral Fermions in a 1D Optical Lattice by Orbital-Momentum Coupling, Physical Review Letters 121, 150403 (2018).
  * Hyoun Kang  _et al._ [2020] J. Hyoun Kang, J. Ho Han, and Y. Shin, Creutz ladder in a resonantly shaken 1D optical lattice, New Journal of Physics 22, 013023 (2020).
  * Hung  _et al._ [2021] J. S. C. Hung, J. H. Busnaina, C. W. S. Chang, A. M. Vadiraj, I. Nsanzineza, E. Solano, H. Alaeian, E. Rico, and C. M. Wilson, Quantum Simulation of the Bosonic Creutz Ladder with a Parametric Cavity, Phys. Rev. Lett. 127, 100503 (2021).
  * Han  _et al._ [2020] J. Han, A. A. Sukhorukov, and D. Leykam, Disorder-protected quantum state transmission through helical coupled-resonator waveguides, Photonics Research 8, 10, B15-B24 (2020).
  * Barbarino  _et al._ [2019] S. Barbarino, D. Rossini, M. Rizzi, R. Fazio, G. E. Santoro, and M. Dalmonte, Topological Devil’s staircase in atomic two-leg ladders, New Journal of Physics 21, 043048 (2019).
  * Resta [1998] R. Resta, Quantum-Mechanical Position Operator in Extended Systems, Physical Review Letters 80, 1800 (1998).
  * Vanderbilt [2018] D. Vanderbilt,  _Berry phase in electronic structure theory_ (Cambridge University Press, 2018).
  * Bena and Montambaux [2009] C. Bena and G. Montambaux, Remarks on the tight-binding model of graphene, New Journal of Physics 11, 095003 (2009).
  * Cayssol and Fuchs [2021] J. Cayssol and J. N. Fuchs, Topological and geometrical aspects of band theory, Journal of Physics: Materials 4, 034007 (2021).
  * Fuchs and Piéchon [2021] J. N. Fuchs and F. Piéchon, Electric polarization of one-dimensional inversion-symmetric two-band insulators (2021), arXiv:2106.03595 .
  * Delplace  _et al._ [2011] P. Delplace, D. Ullmo, and G. Montambaux, Zak phase and the existence of edge states in graphene, Physical Review B 84, 195452 (2011).
  * Longhi [2018] S. Longhi, Probing one-dimensional topological phases in waveguide lattices with broken chiral symmetry, Optics Letters 43, 4639 (2018).
  * Bello  _et al._ [2019] M. Bello, G. Platero, J. I. Cirac, and A. González-Tudela, Unconventional quantum optics in topological waveguide QED, Science Advances 5, eaaw0297 (2019).
  * Aihara  _et al._ [2020] Y. Aihara, M. Hirayama, and S. Murakami, Anomalous dielectric response in insulators with the $\pi$ Zak phase, Phys. Rev. Research 2, 033224 (2020).
  * Pérez-González  _et al._ [2019b] B. Pérez-González, M. Bello, G. Platero, and Á. Gómez-León, Simulation of 1D Topological Phases in Driven Quantum Dot Arrays, Physical Review Letters 123, 126401 (2019b).
  * Jiao  _et al._ [2021] Z.-Q. Jiao, S. Longhi, X.-W. Wang, J. Gao, W.-H. Zhou, Y. Wang, Y.-X. Fu, L. Wang, R.-J. Ren, L.-F. Qiao, and X.-M. Jin, Experimentally Detecting Quantized Zak Phases without Chiral Symmetry in Photonic Lattices, Physical Review Letters 127, 147401 (2021).
  * Cooper  _et al._ [2019] N. R. Cooper, J. Dalibard, and I. B. Spielman, Topological bands for ultracold atoms, Reviews of Modern Physics 91, 015005 (2019).
  * Atala  _et al._ [2013] M. Atala, M. Aidelsburger, J. Barreiro, D. Abanin, T. Kitagawa, E. Demler, and I. Bloch, Direct measurement of the Zak phase in topological Bloch bands, Nature Physics , 795 (2013).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Maximally localized generalized Wannier functions for composite energy bands, Physical Review B 56, 12847 (1997).
  * Creffield [2020] C. E. Creffield, Controlling spin without magnetic fields: The Bloch-Rashba rotator, Phys. Rev. B 102, 085436 (2020).
  * Jafari  _et al._ [2019] R. Jafari, H. Johannesson, A. Langari, and M. A. Martin-Delgado, Quench dynamics and zero-energy modes: the case of the Creutz model, Physical Review B 99, 054302 (2019).
  * Tovmasyan  _et al._ [2018] M. Tovmasyan, S. Peotta, L. Liang, P. Törmä, and S. D. Huber, Preformed pairs in flat Bloch bands, Phys. Rev. B 98, 134513 (2018).
  * Nelson  _et al._ [2021] A. Nelson, T. Neupert, T. c. v. Bzdušek, and A. Alexandradinata, Multicellularity of Delicate Topological Insulators, Phys. Rev. Lett. 126, 216404 (2021).
  * Hofmann and Potthoff [2012] F. Hofmann and M. Potthoff, Doublon dynamics in the extended Fermi-Hubbard model, Physical Review B 85, 205127 (2012).
  * Raoux  _et al._ [2014] A. Raoux, M. Morigi, J.-N. Fuchs, F. Piéchon, and G. Montambaux, From Dia- to Paramagnetic Orbital Susceptibility of Massless Fermions, Phys. Rev. Lett. 112, 026402 (2014).
  * Illes and Nicol [2016] E. Illes and E. J. Nicol, Magnetic properties of the $\alpha-T_{3}$ model: Magneto-optical conductivity and the Hofstadter butterfly, Physical Review B 94, 125435 (2016).
  * Illes and Nicol [2017] E. Illes and E. J. Nicol, Klein tunneling in the $\alpha-T_{3}$ model, Physical Review B 95, 235432 (2017).