1]Moscow Institute of Physics and Technology, Dolgoprudny, Russia 2]Skoltech,
Moscow, Russia 3]Institute for Information Transmission Problems, Moscow,
Russia

# Decentralized optimization with affine constraints over time-varying
networks

Demyan Yarmoshik<EMAIL_ADDRESS>Alexander Rogozin
<EMAIL_ADDRESS>Alexander Gasnikov<EMAIL_ADDRESS>[  [  [

###### Abstract

The decentralized optimization paradigm assumes that each term of a finite-sum
objective is privately stored by the corresponding agent. Agents are only
allowed to communicate with their neighbors in the communication graph. We
consider the case when the agents additionally have local affine constraints
and the communication graph can change over time. We provide the first
linearly convergent decentralized algorithm for time-varying networks by
generalizing the optimal decentralized algorithm ADOM to the case of affine
constraints. We show that its rate of convergence is optimal for first-order
methods by providing the lower bounds for the number of communications and
oracle calls.

###### keywords:

Convex optimization Decentralized optimization Affine constraints Time-varying
networks.

##  1 Introduction

Decentralized optimization is a popular approach for high-dimensional machine
learning problems and control of distributed systems, such as drone swarms. It
is common for these problems to have interconnections between variables,
usually posed as affine constraints, e.g., direct-current(DC) power flow
constraints in control problems related to electrical energy systems. It is
also often the case that the communication graph between computation nodes is
subject to change during the optimization process.

Over the past decade, constrained distributed optimization has attracted the
attention of researchers. Among the first applications of first-order methods
to constrained decentralized optimization was the projected subgradient
algorithm in [1], where the time-varying case was also analyzed. A systematic
review of main problem classes falling into the definition of distributed
constrained optimization along with algorithms working on various levels of
decentralization was given in [2]. More recent works use first-order methods
to deal with a broad range of problem variants, including the nonconvex
objectives [3, 4], the composite objectives [5, 6], the inequality constraints
[7, 5, 6, 3, 8, 9] and other assumptions on problem’s structure [5, 10]. The
ADMM-based approaches are also popular [11, 12, 13].

However, to the best of our knowledge, no decentralized linear convergent
first-order algorithms for affine-constrained problems have been proposed. In
this work, we close this gap by providing a linearly convergent dual algorithm
for decentralized affine-constrained optimization of the sum of smooth
strongly convex functions over time-varying networks. We build on the recently
developed optimal algorithms for decentralized optimization over time-varying
networks [14, 15], and extend these results to the affine-constrained case.
This paper could also be seen as a generalization of [16] to the time-varying
networks.

We also show that our new algorithm inherits the optimality of ADOM by
constructing lower bounds on the number of communications and oracle calls in
time-varying case. During this analysis we also prove lower bounds for the
static communication graph setup, thus showing the optimality of algorithms in
[16].

###  1.1 Basic definitions and assumptions

  * •

Differentiable function $f$ is $\mu$-strongly convex if

$\displaystyle f(y)\geq f(x)+\langle\nabla
f(x),y-x\rangle+\frac{\mu}{2}\left\|y-x\right\|^{2}~{}\forall
x,y\in\mathbb{R}.$ (1)

  * •

Differentiable function $f$ is $L$-smooth (or has $L$-Lipschitz continuous
gradient) if

$\displaystyle f(y)\leq f(x)+\langle\nabla
f(x),y-x\rangle+\frac{L}{2}\left\|y-x\right\|^{2}~{}\forall x,y\in\mathbb{R}.$
(2)

  * •

$\lambda_{\min}(A)$, $\lambda_{\min}^{+}(A),~{}\lambda_{\max}(A)$ are the
minimal, the minimal positive and the maximal eigenvalues of a matrix $A$
respectively.

  * •

$\sigma_{\min}(A)$, $\sigma_{\min}^{+}(A),~{}\sigma_{\max}(A)$ are minimal,
minimal positive and maximal singular values
$\sigma_{i}(A)=\sqrt{\lambda_{i}(A^{\top}A)}$ of a matrix $A$ respectively.

Notation

  * •

“Dimension-lifted” vectors and matrices are written in bold: ${\bf x},{\bf
A}$.

  * •

$x_{[i]}$ denotes $i$-th component of a vector $x$.

  * •

The identity matrix of order $d$ is denoted by $I_{d}$. Sometimes the
subscript is omitted if the order of $I$ is clear from the context.

  * •

$\mathbf{1}_{d}$ denotes the column vector of ones in $\mathbb{R}^{d}$.

  * •

$\operatorname{Im}A$ denotes the image of a linear operator $A$.

##  2 Problem formulation

###  2.1 Objective and constraints

We consider the following affine constrained optimization problem, where
$f_{i}$ are assumed to be $L_{F}$-smooth (have Lipschitz-continuous gradient
with constant $L_{F}$) and $\mu_{F}$-strongly convex:

$\displaystyle\min_{x_{1},\ldots,x_{n}\in\mathbb{R}^{d}}~{}$
$\displaystyle\sum_{i=1}^{n}f_{i}(x_{i})$ (3) s.t. $\displaystyle
A_{i}x_{i}=b_{i},~{}i=1,\ldots,n$ (4) $\displaystyle x_{1}=\ldots=x_{n}.$ (5)

Practical examples of this type of finite-sum affine constrained optimization
problems include constrained estimation problems, such as constrained least
squares problems [17].

The matrix-vector form of the problem is

$\displaystyle\min_{{\bf x}\in\mathbb{R}^{nd}}~{}$ $\displaystyle F({\bf x})$
(6) s.t. $\displaystyle{\bf A}{\bf x}={\bf b}$ $\displaystyle{\bf
x}\in{\mathcal{L}},$

where ${\bf x}$ is the column vector of $x_{1},\ldots,x_{n}$,
$x_{i}\in\mathbb{R}^{d}$, $F({\bf x})=\sum_{i=1}^{n}f_{i}(x_{i})$, and
${\mathcal{L}}$ denotes the consensus hyperlane defined by constraint (5).
Note, that $F({\bf x})$ is also $\mu_{F}$-strongly convex and $L_{F}$-smooth.

In case when $A_{i}$ and $b_{i}$ are all different, we set ${\bf
A}=\operatorname{diag}\left(A_{1},\ldots,A_{n}\right)$, ${\bf
b}=\mathbf{1}\otimes b$, where $\otimes$ denotes the Kronecker product. If all
$A_{i}=A$ and $b_{i}=b$ are the same, then ${\bf A},{\bf b}$ can be chosen in
various ways, e.g.:

  * •

${\bf A}={\bf I}\otimes A=\operatorname{diag}\left(A,\ldots,A\right)$, ${\bf
b}=\mathbf{1}\otimes b$, or

  * •

${\bf A}=(A~{}0~{}\ldots~{}0),~{}{\bf b}=b$.

For definiteness, we will assume that the first variant is chosen. This logic
also applies if there are clusters of agents with same affine constraints in
each group.

###  2.2 Decentralized communication

We assume that the problem is distributed over a computational network
consisting of $n$ agents (or nodes). Each agent privately holds $f_{i}$,
$A_{i}$ and $b_{i}$. Agents are connected through a communication network,
represented by a time-varying undirected graph, i.e. a sequence of undirected
graphs. Agents are only allowed to exchange information with their current
neighbors in the communication graph The limitations imposed on the
communication process are formally described in Defintion 1.

In further developments, we will heavily rely on the notion of gossip matrix.
$W$ is a gossip matrix of an undirected graph $G=(V,E)$, $|V|=n$ if it
satisfies following properties

  1. 1.

${W}$ is symmetric and positive semi-definite;

  2. 2.

(Network compatibility) $[{W}]_{ij}\neq 0$ if and only if $(i,j)\in E$ or
$i=j$;

  3. 3.

(Kernel property) For any $v=[v_{1},\ldots,v_{n}]^{\top}\in\mathbb{R}^{n}$,
${W}v=0$ if and only if $v_{1}=\ldots=v_{n}$.

A typical example of a gossip matrix is the Laplacian matrix of a graph:
$W\in\mathbb{R}^{m\times m}$,

$\displaystyle[W]_{ij}=\begin{cases}-1,&\text{if }(i,j)\in E,\\\
\text{deg}(i),&\text{if }i=j,\\\ 0,&\text{otherwise.}\end{cases}$ (7)

Later we will use the dimension-lifted gossip matrix $W\otimes I_{d}$. From
the third property of gossip matrices, we have that $\left(W\otimes
I_{d}\right){\bf x}=0$ if and only if ${\bf x}\in{\mathcal{L}}$.

Since we assume that the communication network is time-varying, we denote the
communication graph at the $k$-th step by $G(k)=(V,E(k))$ and the associated
gossip matrix by $W(k)$. Note that the existence of a gossip matrix for each
step implies, by the kernel property, that $G(k)$ is connected for all $k$.
According to the second property of gossip matrices, multiplication
$\left(W(k)\otimes I_{d}\right){\bf x}$ can be performed in a decentralized
way at step $k$.

The convergence rate of decentralized optimization algorithms depends on the
spectrum of the gossip matrices. Therefore we assume that

$\tilde{\lambda}_{\min}^{+}\leq\lambda_{\min}^{+}(W(k))\leq\lambda_{\max}(W(k))\leq\tilde{\lambda}_{\max}.$
(8)

##  3 Lower bounds

We consider the class of first-order decentralized algorithms defined as
follows.

######  Definition 1.

Denote ${\mathcal{M}}_{i}(k)$, where ${\mathcal{M}}_{i}(0)=\\{x_{i}^{0}\\}$,
as the local memory of the $i$-th agent at step $k$. The set of allowed
actions of a first order decentralized algorithm at step $k$ is restricted to
the three options

  1. 1.

Local computation: ${\mathcal{M}}_{i}(k)=\operatorname{span}\left(\\{x,\nabla
f_{i}(x),\nabla f^{*}_{i}(x):x\in{\mathcal{M}}_{i}(k)\\}\right)$;

  2. 2.

Decentralized communication:
${\mathcal{M}}_{i}(k)=\operatorname{span}\left(\\{{\mathcal{M}}_{j}(k):\text{edge}~{}(i,j)\in
E(k)\\}\right)$.

  3. 3.

Matrix multiplication:
${\mathcal{M}}_{i}(k)=\operatorname{span}\left(\\{b_{i},A_{i}^{\top}A_{i}x:x\in{\mathcal{M}}_{i}(k)\\}\right)$.

After each step $k$, an algorithm must provide a current approximate solution
$x^{k}_{i}\in{\mathcal{M}}_{i}(k)$ and set
${\mathcal{M}}_{i}(k+1)={\mathcal{M}}_{i}(k)$.

Using the standard approach for constructing lower bounds in smooth strongly
convex optimization [18], we consider $d=\infty$. Then, following [19], we
consider the affine-constrained problem with constraint $(W^{\prime}\otimes
I_{d/m})x=0$, where $W^{\prime}\in\mathbb{R}^{m\times m}$ is a gossip matrix
of some (static) communication graph, thus interpreting an affine-constrained
problem as a decentralized optimization problem.

Let $W^{\prime}$ be a gossip matrix associated with the graph
$G^{\prime}(V^{\prime},E^{\prime})$, $V^{\prime}=\\{1,\ldots,m\\}$. Set
$A_{i}=A=\sqrt{W^{\prime}\otimes I_{m}}$ and $b_{i}=0$. This leads to a two-
level decentralized optimization problem. On the upper level we have the
conventional decentralized optimization problem over the communication network
$G$:

$\displaystyle\min_{x\in\mathbb{R}^{d}}\sum_{i=1}^{n}f_{i}(x),$

but each $f_{i}=\sum_{j=1}^{m}f_{ij}$ is distributed among the subnodes of the
inner computational network $G^{\prime}$ located inside the node $i$. This
forms the inner level of our problem, as shown in Fig. 1. Thus, instead of
thinking about affine constraints, we can think about the inner decentralized
computational network. From this perspective, subnodes exchange the
$d/m$-dimension vectors, and nodes exchange the $d$-dimension vectors, which
are stacked from the vectors of their subnodes.

We use this construction to obtain lower bounds for both static and time-
varying setups.

######  Theorem 1.

For every $L_{F}\geq\mu_{F}>0$, $\chi_{W}>0$ and $\chi_{A}>0$ there exists a
decentralized optimization problem (3) ((6)) over a static communication graph
$G$ with a gossip matrix
$W:\frac{\lambda_{\max}(W)}{\lambda_{\min}^{+}(W)}=\chi_{W}$, and a matrix
$\mathbf{A}:\frac{\lambda_{\max}(\mathbf{A}^{\top}\mathbf{A})}{\lambda_{\min}^{+}(\mathbf{A}^{\top}\mathbf{A})}=\chi_{A}$,
such that any decentralized first-order algorithm, as per Definition 1,
requires at least

$\displaystyle
N=\Theta\left(\sqrt{\frac{L_{F}}{\mu_{F}}}\log\frac{1}{\varepsilon}\right)\text{sequential
local computations},$
$\displaystyle\Theta\left(N\sqrt{\chi_{W}}\right)\text{communications},$
$\displaystyle\Theta\left(N\sqrt{{\chi_{A}}}\right)\text{multiplications by
$\mathbf{A}^{\top}\mathbf{A}$}.$

to achieve $\left\|x_{i}^{N}-x^{*}\right\|^{2}\leq\varepsilon~{}\forall
i\in[1,n]$.

$G^{\prime}$$G$ Figure 1: Two-level computational network

The proof is based on the technique from [20] and is provided in the Appendix.

We now present a variant of Theorem 1 for the time-varying communication
networks.

######  Theorem 2.

For every $L_{F}\geq\mu_{F}>0$, $\chi_{A}>0$, $\tilde{\lambda}_{\min}^{+}$ and
$\tilde{\lambda}_{\max}:$
$\frac{\tilde{\lambda}_{\max}}{\tilde{\lambda}_{\min}^{+}}=\chi_{W}\geq 3$
there exists a decentralized optimization problem (3) ((6)), a sequence of
gossip matrices $W(k)$ satisfying (8), and a matrix
$\mathbf{A}:\frac{\lambda_{\max}(\mathbf{A}^{\top}\mathbf{A})}{\lambda_{\min}^{+}(\mathbf{A}^{\top}\mathbf{A})}=\chi_{A}$,
such that any decentralized first-order algorithm, as per Definition 1,
requires at least

$\displaystyle
N=\Theta\left(\sqrt{\frac{L_{F}}{\mu_{F}}}\log\frac{1}{\varepsilon}\right)\text{sequential
local computations},$
$\displaystyle\Theta\left(N\chi_{W}\right)\text{communications},$
$\displaystyle\Theta\left(N\sqrt{{\chi_{A}}}\right)\text{multiplications by
$\mathbf{A}^{\top}\mathbf{A}$}.$

to achieve $\left\|x_{i}^{N}-x^{*}\right\|^{2}\leq\varepsilon~{}\forall
i\in[1,n]$.

The proof is based on [15] and is provided in the Appendix.

##  4 Application of ADOM

By standard duality arguments we rewrite problem (6) as

$\displaystyle\min_{{\bf x}}$ $\displaystyle\max_{{\bf
s}\in{\mathcal{L}}^{\perp},{\bf p}}\left[F({\bf x})-\left\langle{\bf s},{\bf
x}\right\rangle-\left\langle{\bf p},{\bf A}{\bf x}-{\bf
b}\right\rangle\right]$ $\displaystyle=\max_{{\bf
s}\in{\mathcal{L}}^{\perp},{\bf p}}\left[-\max_{\bf x}\left[\left\langle{\bf
s}+{\bf A}^{\top}{\bf p},{\bf x}\right\rangle-F({\bf
x})\right]+\left\langle{\bf p},{\bf b}\right\rangle\right]$
$\displaystyle=\max_{{\bf s}\in{\mathcal{L}}^{\perp},{\bf p}}\left[-F^{*}({\bf
s}+{\bf A}^{\top}{\bf p})+\left\langle{\bf p},{\bf b}\right\rangle\right]$
$\displaystyle=-\min_{{\bf s}\in{\mathcal{L}}^{\perp},{\bf p}}\left[F^{*}({\bf
s}+{\bf A}^{\top}{\bf p})-\left\langle{\bf p},{\bf b}\right\rangle\right],$

where $F^{*}$ is the convex (Fenchel) conjugate of $F$. Therefore, problem (6)
is equivalent to

$\displaystyle\min_{{\bf s}\in{\mathcal{L}}^{\perp},{\bf p}}\left[F^{*}({\bf
s}+{\bf A}^{\top}{\bf p})-\left\langle{\bf p},{\bf b}\right\rangle\right].$

Introduce

$\displaystyle{\bf z}=\begin{pmatrix}{\bf p}\\\ {\bf s}\end{pmatrix},~{}{\bf
q}=\begin{pmatrix}{\bf b}\\\ 0\end{pmatrix},~{}{\bf B}=\begin{pmatrix}{\bf
A}\\\ I_{nd}\end{pmatrix},~{}{\bf P}=\begin{pmatrix}I_{\dim{({\bf b})}}&0\\\
0&I_{nd}-\frac{1}{n}\left(\mathbf{1}_{n}\mathbf{1}_{n}^{\top}\right)\otimes
I_{d}\end{pmatrix},~{}$ $\displaystyle{\bf W}(k)=\begin{pmatrix}I_{\dim({\bf
b})}&0\\\ 0&W(k)\otimes I_{d}\end{pmatrix}.$

Also denote $H({\bf z})=F^{*}({\bf B}^{\top}{\bf z})-\left\langle{\bf z},{\bf
q}\right\rangle$. Now we can equivalently rewirte optimization problem (6) as

$\displaystyle\min_{{\bf z}\in\operatorname{Im}{\bf P}}~{}H({\bf z}).$ (9)

After that, we apply ADOM [14] to the problem (9).

Algorithm 1 ADOM with local affine constraints

0: ${\bf z}^{0}\in\operatorname{Im}{\bf P}{\bf B},~{}{\bf
m}^{0}\in\operatorname{Im}{\bf P}{\bf
B},~{}\alpha,\eta,\theta,\sigma>0,~{}\tau\in(0,1)$

1: Set ${\bf z}_{f}^{0}={\bf z}^{0}$

2: for $k=0,1,\ldots$ do

3: ${\bf z}_{g}^{k}=\tau{\bf z}^{k}+(1-\tau){\bf z}_{f}^{k}$

4: $\Delta^{k}=\sigma{\bf W}(k)({\bf m}^{k}-\eta\nabla H({\bf z}_{g}^{k}))$

5: ${\bf m}^{k+1}={\bf m}^{k}-\eta\nabla H({\bf z}_{g}^{k})-\Delta^{k}$

6: ${\bf z}^{k+1}={\bf z}^{k}+\eta\alpha({\bf z}_{g}^{k}-{\bf
z}^{k})+\Delta^{k}$

7: ${\bf z}_{f}^{k+1}={\bf z}_{g}^{k}-\theta{\bf W}(k)\nabla H({\bf
z}_{g}^{k})$

8: end for

Note that $b_{i}\in\operatorname{Im}A_{i}$ and therefore ${\bf
q}\in\operatorname{Im}{\bf B}$. Therefore, for any ${\bf z}\in\mathbb{R}^{nd}$
we have $\nabla H({\bf z})={\bf B}\nabla F^{*}({\bf B}^{\top}{\bf z})-{\bf
q}\in\operatorname{Im}{\bf B}$. We conclude that the iterates ${\bf
z}^{k},\Delta^{k},{\bf z}_{f}^{k},{\bf z}_{g}^{k}$ of Algorithm 1 lie in
$\operatorname{Im}{\bf P}{\bf B}=\operatorname{Im}{\bf
A}\otimes{\mathcal{L}}^{\bot}$, and ${\bf m}^{k}\in\operatorname{Im}{\bf
B}=\operatorname{Im}{\bf A}\otimes\mathbb{R}^{nd}$.

On the subspace $\operatorname{Im}{\bf P}{\bf B}$ we estimate strong convexity
and smoothness constants as

$\displaystyle\mu_{H}=\frac{1+(\sigma_{\min}^{+}(A))^{2}}{L_{F}},~{}L_{H}=\frac{1+\sigma_{\max}^{2}(A)}{\mu_{F}}.$

From (8) we have

$\lambda_{\min}^{+}\vcentcolon=\min(1,\tilde{\lambda}_{\min}^{+})\leq\lambda_{\min}^{+}(\mathbf{W}(k))\leq\lambda_{\max}(\mathbf{W}(k))\leq\max(1,\tilde{\lambda}_{\max})=\vcentcolon\lambda_{\max}.$
(10)

Now we can formulate the key convergence result.

######  Theorem 3.

Set parameters $\alpha,\eta,\theta,\sigma,\tau$ of Algorithm 1 to
$\alpha=\frac{\mu_{H}}{2}$,
$\eta=\frac{2\lambda_{\min}^{+}}{7\lambda_{\max}\sqrt{\mu_{H}L_{H}}}$,
$\theta=\frac{1}{L_{H}\lambda_{\max}}$, $\sigma=\frac{1}{\lambda_{\max}}$, and
$\tau=\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu_{H}}{L_{H}}}$.
Then there exists $C>0$, such that for Algorithm 1 applied for problem (6) it
holds

$\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf z}_{g}^{N})-{\bf
x}^{*}\right\|^{2}\leq
C\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu_{H}}{L_{H}}}\right)^{N},$
(11)

where ${\bf x}^{*}$ is the solution of the problem 6.

The proof is a minor modification of original ADOM convergence proof in [14],
and is provided in the appendix for reader’s convenience.

As a corollary of Theorem 3 we have the following communication, dual-oracle
call and matrix multiplication complexity:

$\displaystyle N\leq
O\left(\frac{\max(1,\tilde{\lambda}_{\max})}{\min(1,\tilde{\lambda}_{\min}^{+})}\sqrt{\frac{L_{F}}{\mu_{F}}}\sqrt{\frac{1+\sigma_{\max}^{2}(A)}{1+(\sigma_{\min}^{+}(A))^{2}}}\log\left(\frac{1}{\varepsilon}\right)\right),$

where $\varepsilon$ is the desired accuracy of the approximate solution:
$\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf z}_{g}^{N})-{\bf
x}^{*}\right\|^{2}\leq\varepsilon$.

##  5 Separating complexities

As shown in [15] one can separate oracle (computation) and communication
complexities in the time-varying setup by using the multi-consensus procedure.
This will not change (up to a $\ln 2$ factor) the number of communications but
decreases the number of oracle calls.

In the time-varying setup, acceleration over
$\chi=\frac{\lambda_{\max}}{\lambda_{\min}^{+}}$ is not applicable, as stated
by Theorem 2, so we can separate communication complexity but can not improve
it. However, since $\mathbf{A}$ is constant, we can use Chebyshev acceleration
over $\mathbf{A}$ to separate matrix multiplication complexity and decrease
the number of multiplications by $\mathbf{A}$, as was done for communication
complexity in the static communication graph setup in [20].

For this section we will assume that $W(k)$ are divided by their maximum
eigenvalue and therefore $\lambda_{\max}=\tilde{\lambda}_{\max}=1$ and
$\lambda_{\min}^{+}=\tilde{\lambda}_{\min}^{+}$. In practice this could be
achieved by using $I-M(k)$ instead of $W(k)$, there $M(k)$ are the Metropolis
weight matrices [15].

Let $P(x)$ be a polynomial such that $P(0)=0$ and $P(\lambda_{i})\neq 0$ for
all eigenvalues $\lambda_{i}$ of $A$. Then, using the fact that $\ker M=\ker
M^{T}M$, we can do the following sequence of equivalent reformulations

$\displaystyle\mathbf{A}{\bf x}={\bf b}$
$\displaystyle\Leftrightarrow\mathbf{A}{\bf x}=\mathbf{A}{\bf
x}_{0}\Leftrightarrow\mathbf{A}^{\top}\mathbf{A}{\bf
x}=\mathbf{A}^{\top}\mathbf{A}{\bf x}_{0}$ $\displaystyle\Leftrightarrow
P(\mathbf{A}^{\top}\mathbf{A}){\bf
x}=P^{\prime}(\mathbf{A}^{\top}\mathbf{A})\mathbf{A}^{\top}\mathbf{A}{\bf
x}_{0}$ $\displaystyle\Leftrightarrow P(\mathbf{A}^{\top}\mathbf{A}){\bf
x}=P^{\prime}(\mathbf{A}^{\top}\mathbf{A})\mathbf{A}^{\top}{\bf b}.$

where ${\bf x}_{0}$ is any vector satisfying $\mathbf{A}{\bf x}_{0}={\bf b}$
(here we used consistency of constraints), and $P^{\prime}(x)=P(x)/x$ is
correctly defined since $P(0)=0$.

Thus the idea is to replace $\mathbf{A}$ with $P(\mathbf{A}^{\top}\mathbf{A})$
to improve the spectral properties of the matrix. The polynomials of choice
are shifted and scaled Chebyshev polynomials [20], because Chebyshev
polynomials increase magnitude more quickly than any other polynomials of the
same degree satisfying $|P(x)|\leq 1~{}\forall x\in[-1,1]$. This allows to
significantly compress the spectrum, using polynomials of a relatively low
degree. In particular, let $P(x)$ be defined as

$\displaystyle
P(x)=1-\frac{T_{K}\left(-\nu+\frac{2}{\lambda_{\max}-\lambda_{\min}^{+}}x\right)}{T_{K}(-\nu)},~{}K=\left\lfloor{\sqrt{\chi_{A}}}\right\rfloor,$

where
${\chi_{A}}=\sigma_{\max}(\mathbf{A}^{T}\mathbf{A})/\sigma_{\min}^{+}(\mathbf{A}^{T}\mathbf{A})$,
$\nu=\frac{\chi+1}{\chi-1}$ and $T_{K}$ are Chebyshev polynomials defined by
$T_{0}(x)=1$, $T_{1}(x)=x$ and $T_{K+1}=2xT_{K}(x)-T_{K-1}(x)$ for $K\geq 1$.
If $P$ has degree $K=\left\lfloor{\sqrt{\chi_{A}}}\right\rfloor$, then
$\chi\left(P({\mathbf{A}^{\top}}\mathbf{A})\right)\leq 4$ [20], which allows
to quadratically improve the dependence of convergence rate on $\chi_{A}$ by
replacing $\mathbf{A}$ with $P(\mathbf{A}^{\top}\mathbf{A})$ and ${\bf b}$
with $P^{\prime}({\mathbf{A}^{\top}}\mathbf{A}){\mathbf{A}^{\top}}{\bf b}$.

To separate communication complexity, a multi-consensus procedure should be
used, i.e. $W(k)$ should be replaced with $D(W(k))=I-(I-W(k))^{K}$, where
$K=\left\lceil{\chi\ln 2}\right\rceil$, what makes $\chi=O(1)$ at the cost of
$K$ communication rounds [15].

Finally, by replacing $\mathbf{A}\to P({\mathbf{A}^{\top}}\mathbf{A})$, ${\bf
b}\to P^{\prime}({\mathbf{A}^{\top}}\mathbf{A}){\mathbf{A}^{\top}}{\bf b}$ and
$W(k)\to D(W(k))$ we obtain following complexity estimates to reach
$\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf z}_{g}^{k})-{\bf
x}^{*}\right\|^{2}\leq\varepsilon$ for ADOM algorithm with affine constraints,
Chebyshev acceleration, and multi-consensus:

$\displaystyle
N=O\left(\sqrt{\frac{L_{F}}{\mu_{F}}}\log\frac{1}{\varepsilon}\right)\text{
oracle calls at each node},$ $\displaystyle O\left(N\chi\right)\text{
communications},$ $\displaystyle O\left(N\sqrt{{\chi_{A}}}\right)\text{
multiplications by $\mathbf{A}^{\top}\mathbf{A}$}.$

These upper bounds match the lower bounds of Theorem 2, thus the obtained
algorithm is optimal among first-order decentralized algorithms for strongly
convex problems with affine constraints on time-varying networks.

##  6 Numeric validation

We verify Theorem 3 with numeric experiments on problems with the quadratic
objective:

$\displaystyle\sum_{i}^{n}f_{i}(x)\to\min_{x}$ s.t. $\displaystyle Ax=b,$
$\displaystyle f_{i}(x)$
$\displaystyle=\frac{1}{2}x^{T}C_{i}x+d_{i}^{T}x,\;\mu_{F}I\preceq
C_{i}\preceq L_{F}I.$

From the results of Section 5, the influence of the gossip matrix’s spectrum
and the affine constraint matrix spectrum on the convergence rates of
Algorithm 1 is straightforward to comprehend. Therefore, we focus only on the
impact of the objective’s condition number on the convergence rates.

Our numeric experiments do not pretend to trial on real-world problems, but
rather illustrations of the theoretical properties of the algorithm, because
quadratic objects are a good representative of the smooth and strongly convex
problem class.

It is not difficult to implement exact dual oracle for this objective, but we
do not want to exploit the simplicity of quadratic problem and, following
[14], we obtain an approximation of ${\bf g}^{k}\approx\nabla
F^{*}({\mathbf{B}^{\top}}{\bf z}_{g}^{k})$ by using few gradient steps at each
iteration: ${\bf g}^{k}={\bf g}^{k-1}-\frac{1}{L_{F}}(\nabla F({\bf
g}^{k-1})-{\mathbf{B}^{\top}}{\bf z}_{g}^{k})$. So in fact the implemented
algorithm uses a primal oracle because dual oracle call in Algorithm 1 is
replaced with primal oracle call.

Experiment parameters are $d=20$, $A\in\mathbb{R}^{10\times d}$,
${\chi_{A}}=20$, $\mu_{F}=1$. Communication graphs $G(k)$ are random ring
graphs at each iteration. We run Algorithm 1 for $N=2500$ iterations for
different values of $L_{F}$, and do the linear regression to obtain the
coefficient $\kappa$ in the dependence $\|{\bf g}^{k}-{\bf
x}^{*}\|_{2}=C_{1}\exp(-\kappa N)$ using only last $N/2$ iterations. This is
illustrated in Figure 2. In all cases a steady linear convergence to the
solution is present.

Then we do the linear regression in the log-log scale to obtain the
coefficient $\nu$ in the dependence
$\kappa=C_{2}\left(\frac{L_{F}}{\mu_{F}}\right)^{\nu}$, as shown in Figure 3.
The resulting value is $\nu\approx 0.54$ with the standard error being
$\approx 0.03$, what is rather close to the value $\nu=\frac{1}{2}$ in the
Theorem 3.

Figure 2: Algorithm 1 convergence for different values of $L_{F}$ Figure 3:
Regression for $\nu$

##  7 Conclusion

By viewing the affine-constrained problem as a decentralized optimization
problem alike [19], and combining constructions of lower bound for static [20]
and time-varying [15] setups, lower bounds for decentralized optimization with
affine constraints over static and time-varying networks via first-order
methods were obtained.

As we found, the ADOM algorithm can be straightforwardly extended to the
affine-constrained case. For this problem class, we were also able to apply
Chebyshev acceleration over $A$, and the resulting complexity estimates match
the lower bounds.

However, a lot of questions are left for the future work. We did not succeed
to provide an extension of the ADOM+ algorithm [15] to the affine-constrained
problems, thus no linearly convergent primal algorithm is known for this
problem class . It is also of interest to obtain optimal algorithms for time-
varying networks in case of shared affine inequality constraints ($\sum
A_{i}x_{i}\leq b_{i}$, where $A_{i}$ and $b_{i}$ are held privately by $i$-th
agent). This problem variant has more practical applications [5, 21, 2], but
also brings additional difficulties to the theoretical analysis, e.g. in this
case $b_{i}$ might not belong to $\operatorname{Im}A_{i}$. This means that
cannot apply our approach, because it needs a gradient method to stay in the
subspace where the objective is strongly convex.

##  8 Acknowledgements

This work was supported by a grant for research centers in the field of
artificial intelligence, provided by the Analytical Center for the Government
of the Russian Federation in accordance with the subsidy agreement (agreement
identifier 000000D730321P5Q0002) and the agreement with the Moscow Institute
of Physics and Technology dated November 1, 2021 No. 70-2021-00138.

##  9 Appendix

###  9.1 Proof of Theorem 1

###### Proof.

Let the affine constraint in problem 3 be $A_{i}x_{i}=0$, with
$A_{i}=A=\sqrt{W^{\prime}\otimes I_{m}}$. Then the affine constrained
decentralized problem can be seen as two-level decentralized problem, as
explained above.

Select sets of subnodes $S_{1}$, $S_{2}$ and $S_{3}$ such that $S_{1}$,
$S_{2}$ are at the distance $\geq\Delta_{A}$ through the inner graph, and
$S_{2}$, $S_{3}$ are at the distance $\geq\Delta_{W}$ through the outer graph.
Consider the following splitting of the Nesterov’s “bad” function

$\displaystyle
f_{ij}(x)=\frac{\alpha}{2mn}\left\|x\right\|^{2}+\frac{\beta-\alpha}{8}\cdot\begin{cases}\frac{1}{|S_{1}|}\left(x^{\top}M_{1}x-2x_{[1]}\right),~{}(i,j)\in
S_{1},\\\ \frac{1}{|S_{2}|}x^{\top}M_{2}x,~{}(i,j)\in S_{2},\\\
\frac{1}{|S_{3}|}x^{\top}M_{3}x,~{}(i,k)\in S_{3},\\\
0,~{}\text{otherwise,}\end{cases}$ (12)

where

$\displaystyle M_{1}$
$\displaystyle=\operatorname{diag}{(1,0,M_{0},0,M_{0},\ldots)},$
$\displaystyle M_{2}$
$\displaystyle=\operatorname{diag}{(M_{0},0,M_{0},0,\ldots)},$ $\displaystyle
M_{3}$ $\displaystyle=\operatorname{diag}{(0,M_{0},0,M_{0},\ldots)},$
$\displaystyle M_{0}$ $\displaystyle=\begin{pmatrix}1&-1\\\
-1&1\end{pmatrix}.$

Then increasing the number of nonzero components in $x^{k}$ on any subnode by
three requires one local computation on a node in $S_{1}$, $\Delta_{A}$ inner
communications a.k.a. multiplications by $A^{\top}A$, one local computation on
a node in $S_{2}$, $\Delta_{W}$ communications in the outer graph and one
local computation on a node in $S_{3}$. Denote by
$\kappa_{g}=\frac{\beta}{\alpha}$ the “global” condition number of
$\sum_{ij}^{mn}f_{ij}$. Since the solution is
$x^{*}_{k}=\left(\frac{\sqrt{\kappa}_{g}-1}{\sqrt{\kappa}_{g}+1}\right)^{k}$,
we have

$\left\|x^{N}-x^{*}\right\|^{2}\geq\sum_{k=N+2}^{\infty}(x^{*}_{k})^{2}\geq\left(\frac{\sqrt{\kappa}_{g}-1}{\sqrt{\kappa}_{g}+1}\right)^{N+2},$
(13)

where $N$ is the number of iterations, each including 3 sequential
computational steps, $\Delta_{A}$ multiplications by $A^{T}A$ and $\Delta_{W}$
communications.

To finish the proof we need to construct communication graphs $G$,
$G^{\prime}$, where distances between $S_{1}$, $S_{2}$ and $S_{2}$, $S_{3}$
are close to $\Delta_{A},\Delta_{W}$, and equip the graphs with gossip
matrices with given condition numbers $\chi_{A},\chi_{W}$.

We should also choose $\alpha$ and $\beta$ such that $f_{i}$ are
$L_{F}$-smooth and $\mu_{F}$-strongly convex, and choose $S_{1}$, $S_{2}$,
$S_{3}$ so that $\kappa_{g}$ is similar to $\frac{L_{F}}{\mu_{F}}$.

Denote $\gamma(M)=\sigma_{\min}^{+}(M)/\sigma_{\max}(M)$,
$\gamma_{W}=1/\chi_{W}$, $\gamma_{A}=1/\chi_{A}$. Let
$\gamma_{n}=\frac{1-\cos\left(\frac{\pi}{n}\right)}{1+\cos\left(\frac{\pi}{n}\right)}$
be a decreasing sequence of positive numbers. Since $\gamma_{2}=1$ and
$\lim_{n}\gamma_{n}=$ 0, there exists $n\geq 2$ such that
$\gamma_{n}\geq\gamma>\gamma_{n+1}$ and $m\geq 2$ such that
$\gamma_{m}\geq\gamma>\gamma_{m+1}$.

First, construct graph $G$. The cases $n=2$ and $n\geq 3$ are treated
separately. If $n\geq 3$, let $G$ be the linear graph of size $n$ ordered from
node $1$ to $n$, and weighted with $w_{i,i+1}=\begin{cases}1-a,i=1\\\
1,\text{otherwise}\end{cases}.$ Then set $S_{2G}=\left\\{1,\ldots,\lceil
n/32\rceil\right\\}$ and $\Delta_{W}=(1-1/16)n-1$, so that $S_{3G}=\\{\lceil
n/32\rceil+\lceil\Delta_{W}\rceil,\ldots,n\\}$.

Take $W_{a}$ as the Laplacian of the weighted graph $G$. A simple calculation
gives that, if $a=0$, $\gamma\left(W_{a}\right)=\gamma_{n}$ and, if $a=1$, the
network is disconnected and $\gamma\left(W_{a}\right)=0$. Thus, by continuity
of the eigenvalues of a matrix, there exists a value $a\in[0,1]$ such that
$\gamma\left(W_{a}\right)=\gamma_{W}$. Finally, by definition of $n$, one has
$\gamma_{W}>\gamma_{n+1}\geq\frac{2}{(n+1)^{2}}$, and
$\Delta_{W}\geq\frac{15}{16}\left(\sqrt{\frac{2}{\gamma_{W}}}-1\right)-1\geq\frac{1}{5\sqrt{\gamma_{W}}}$
when $\gamma_{W}\leq\gamma_{3}=\frac{1}{3}$.

For the case $n=2$, we consider the totally connected network of 3 nodes,
reweight only the edge $\left(1,3\right)$ by $a\in[0,1]$, and let $W_{a}$ be
its Laplacian matrix. If $a=1$, then the network is totally connected and
$\gamma\left(W_{a}\right)=1$. If, on the contrary, $a=0$, then the network is
a linear graph and $\gamma\left(W_{a}\right)=\gamma_{3}$. Thus, there exists a
value $a\in[0,1]$ such that $\gamma\left(W_{a}\right)=\gamma$. Set
$S_{2G}=\left\\{1\right\\},S_{3G}=\left\\{2\right\\}$, then
$\Delta_{W}=1\geq\frac{1}{\sqrt{3\gamma_{W}}}$.

Second, do the same for graph $G^{\prime}$, obtaining $m$,
$S_{1G^{\prime}},S_{2G^{\prime}}$ and
$\Delta_{A}\geq\frac{1}{5\sqrt{\gamma_{A}}}$.

Define $S_{1}=S_{2G}\times S_{1G^{\prime}}$, $S_{2}=S_{2G}\times
S_{2G^{\prime}}$ and $S_{3}=S_{3G}\times S_{2G^{\prime}}$, see Fig. 4. In all
cases we have
$|S_{k}|\geq|S_{2}|\geq\lceil\frac{n}{32}\rceil\lceil\frac{m}{32}\rceil$ for
$k\in\\{1,3\\}$.

Because $\mu_{F}=\frac{\alpha}{n}$, we set $\alpha=\mu_{F}n$. Since $0\preceq
M_{k}\preceq 2I$ for $k\in\\{1,2,3\\}$,
$L_{F}=\frac{\alpha}{n}+\frac{(\beta-\alpha)m}{2|S_{2}|}$, thus set
$\beta=2|S_{2}|(L_{F}-\mu_{F})/m+\mu_{F}n$ to make all $f_{i}$ be
$L_{F}$-smooth and $\mu_{F}$-strongly convex. Then
$\kappa_{g}=\frac{\beta}{\alpha}=1+\frac{2|S_{2}|(L_{F}-\mu_{F})}{\mu_{F}mn}\geq\frac{L_{F}}{512\mu_{F}}$.
Combining this with (13) and the inequalities between $\Delta_{A},\gamma_{A}$
and $\Delta_{W},\gamma_{W}$ we conclude the proof. ∎

………………………⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮$S_{2}$$S_{1}$$S_{3}$$n$$m$ Figure 4: Splitting of
Nesterov’s “bad” function in the main case $m,n\geq 3$, namely $m=20,n=80$.
Vertical dimension corresponds to inner graph $G^{\prime}$, and horizontal
dimension — to outer graph $G$

###  9.2 Proof of Theorem 2

###### Proof.

As in the proof of Theorem 1 we set the affine constraint in problem 3 to be
$A_{i}x_{i}=0$, with $A_{i}=A=\sqrt{W^{\prime}\otimes I_{m}}$, where
$W^{\prime}$ is a gossip matrix of some inner communication graph
$G^{\prime}$. Let the sequence of outer communication graphs $G(k)$ is the
same as in the proof of Theorem 1 in [15]:
$n=3\left\lfloor{\chi_{W}/3}\right\rfloor$ nodes are split into three disjoint
sets $V_{1},V_{2},V_{3}$ of equal size, and $G(k)=(V,E(k))$ are star graphs
with the center nodes cycling through $V_{2}$. Choose the inner communicaiton
graph $G^{\prime}=(V^{\prime},E^{\prime})$ as in the proof of Theorem 1. Use
Nesterov’s function splitting given by (12), choose $S_{1G^{\prime}}$ and
$S_{2G^{\prime}}$ as in the proof of Theorem 1. Set $S_{1}=V_{1}\times
S_{1G^{\prime}}$, $S_{2}=V_{1}\times S_{2G^{\prime}}$ and $S_{3}=V_{3}\times
E^{\prime}$. Setting $W(k)$ to be the Laplacian of the star graph $G(k)$ we
have $\frac{\lambda_{\max}(W(k))}{\lambda_{\min}(W(k))}=n\leq\chi_{W}$. Also
(Lemma 2 [15] and proof of Theorem 1), increasing the number of nonzero
components of $x_{k}$ on any subnode requires local computation on a node in
$S_{1}$, $\Theta\left(\sqrt{\chi_{A}}\right)$ communications in the inner
graph $G^{\prime}$ (i.e. multiplications by $A^{\top}A$), one local
computation on a node in $S_{2}$, $\Theta\left(\chi_{W}\right)$ communications
in the outer graph $G$ and one local computation on a node in $S_{3}$. Same
reasoning as in the proof of the previous theorem gives
$\kappa_{g}=\Theta\left(\frac{L_{F}}{\mu_{F}}\right)$, then using (13) we
conclude the proof.

∎

###  9.3 Auxiliary lemmas for Theorem 3

######  Lemma 1.

For $\theta\leq\frac{1}{L_{H}\lambda_{\max}}$ we have the inequality

$H({\bf z}_{f}^{k+1})\leq H({\bf
z}_{g}^{k})-\frac{\theta\lambda_{\min}^{+}}{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}.$ (14)

###### Proof.

We start with $L_{H}$-smoothness of $H$ on
$\operatorname{Im}\mathbf{P}\mathbf{B}$:

$\displaystyle H({\bf z}_{f}^{k+1})\leq H({\bf z}_{g}^{k})+\langle\nabla
H({\bf z}_{g}^{k}),{\bf z}_{f}^{k+1}-{\bf
z}_{g}^{k}\rangle+\frac{L_{H}}{2}\left\|{\bf z}_{f}^{k+1}-{\bf
z}_{g}^{k}\right\|^{2}.$

Using line 7 of Algorithm 1 together with (10) we get

$\displaystyle H({\bf z}_{f}^{k+1})$ $\displaystyle\leq H({\bf
z}_{g}^{k})-\theta\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{W}(k)}+\frac{L_{H}\theta^{2}}{2}\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{W}^{2}(k)}$ $\displaystyle\leq H({\bf
z}_{g}^{k})-\frac{\theta\lambda_{\min}^{+}}{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}-\frac{\theta}{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{W}(k)}+\frac{L_{H}\theta^{2}\lambda_{\max}}{2}\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{W}(k)}$ $\displaystyle=H({\bf
z}_{g}^{k})-\frac{\theta\lambda_{\min}^{+}}{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}+\frac{\theta}{2}\left(\theta
L_{H}\lambda_{\max}-1\right)\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{W}(k)}.$

Using condition $\theta\leq\frac{1}{L_{H}\lambda_{\max}}$ we get

$\displaystyle H({\bf z}_{f}^{k+1})$ $\displaystyle\leq H({\bf
z}_{g}^{k})-\frac{\theta\lambda_{\min}^{+}}{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}.$

∎

######  Lemma 2.

For $\sigma\leq\frac{1}{\lambda_{\max}}$ we have the inequality

$\begin{split}&\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}\leq\left(1-\frac{\sigma\lambda_{\min}^{+}}{4}\right)\frac{4}{\sigma\lambda_{\min}^{+}}\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}\\\
&-\frac{4}{\sigma\lambda_{\min}^{+}}\left\|{\bf
m}^{k+1}\right\|^{2}_{\mathbf{P}}+\frac{8\eta^{2}}{(\sigma\lambda_{\min}^{+})^{2}}\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}.\end{split}$ (15)

###### Proof.

Using $\mathbf{P}=\mathbf{P}^{2}$ and
$\mathbf{P}\mathbf{W}(k)=\mathbf{W}(k)\mathbf{P}=\mathbf{W}(k)$ together with
lines 4 and 5 of Algorithm 1 we obtain

$\displaystyle\left\|{\bf m}^{k+1}\right\|^{2}_{\mathbf{P}}$
$\displaystyle=\left\|{\bf m}^{k}-\eta\nabla H({\bf
z}_{g}^{k})-\Delta^{k}\right\|^{2}_{\mathbf{P}}$
$\displaystyle=\left\|(\mathbf{P}-\sigma\mathbf{W}(k))({\bf m}^{k}-\eta\nabla
H({\bf z}_{g}^{k}))\right\|^{2}$ $\displaystyle=\left\|{\bf m}^{k}-\eta\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}-2\sigma\left\|{\bf
m}^{k}-\eta\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{W}(k)}+\sigma^{2}\left\|{\bf
m}^{k}-\eta\nabla H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{W}^{2}(k)}.$

Using (10) we obtain

$\displaystyle\left\|{\bf m}^{k+1}\right\|^{2}_{\mathbf{P}}$
$\displaystyle\leq\left\|{\bf m}^{k}-\eta\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}-\sigma\lambda_{\min}^{+}\left\|{\bf
m}^{k}-\eta\nabla H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle-\sigma\left\|{\bf m}^{k}-\eta\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{W}(k)}+\sigma^{2}\lambda_{\max}\left\|{\bf
m}^{k}-\eta\nabla H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{W}(k)}$
$\displaystyle=\left\|{\bf m}^{k}-\eta\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}-\sigma\lambda_{\min}^{+}\left\|{\bf
m}^{k}-\eta\nabla H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle+\sigma(\sigma\lambda_{\max}-1)\left\|{\bf m}^{k}-\eta\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{W}(k)}.$

Using condition $\sigma\leq\frac{1}{\lambda_{\max}}$ we get

$\displaystyle\left\|{\bf m}^{k+1}\right\|^{2}_{\mathbf{P}}$
$\displaystyle\leq(1-\sigma\lambda_{\min}^{+})\left\|{\bf m}^{k}-\eta\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}.$

Using Young’s inequality we get

$\displaystyle\left\|{\bf m}^{k+1}\right\|^{2}_{\mathbf{P}}$
$\displaystyle\leq(1-\sigma\lambda_{\min}^{+})\left(\left(1+\frac{\sigma\lambda_{\min}^{+}}{2(1-\sigma\lambda_{\min}^{+})}\right)\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+\left(1+\frac{2(1-\sigma\lambda_{\min}^{+})}{\sigma\lambda_{\min}^{+}}\right)\left\|\eta\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}\right)$
$\displaystyle=\left(1-\frac{\sigma\lambda_{\min}^{+}}{2}\right)\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+\eta^{2}\frac{(1-\sigma\lambda_{\min}^{+})(2-\sigma\lambda_{\min}^{+})}{\sigma\lambda_{\min}^{+}}\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle\leq\left(1-\frac{\sigma\lambda_{\min}^{+}}{2}\right)\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+\frac{2\eta^{2}}{\sigma\lambda_{\min}^{+}}\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}.$

Rearranging concludes the proof. ∎

######  Lemma 3.

Let

$\alpha=\frac{\mu_{H}}{2},$ (16)
$\eta=\frac{2\lambda_{\min}^{+}}{7\lambda_{\max}\sqrt{\mu_{H}L_{H}}},$ (17)
$\theta=\frac{1}{L_{H}\lambda_{\max}},$ (18)
$\sigma=\frac{1}{\lambda_{\max}},$ (19)
$\tau=\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu_{H}}{L_{H}}}.$
(20)

Define the Lyapunov function

$\Psi^{k}\vcentcolon=\left\|\hat{{\bf z}}^{k}-{\bf
z}^{*}\right\|^{2}+\frac{2\eta(1-\eta\alpha)}{\tau}(F^{*}({\bf
z}_{f}^{k})-F^{*}({\bf z}^{*}))+6\left\|{\bf m}^{k}\right\|^{2}_{\mathbf{P}},$
(21)

where $\hat{{\bf z}}^{k}$ is defined as

$\hat{{\bf z}}^{k}={\bf z}^{k}+\mathbf{P}{\bf m}^{k}.$ (22)

Then the following inequality holds:

$\Psi^{k+1}\leq\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu_{H}}{L_{H}}}\right)\Psi^{k}.$
(23)

###### Proof.

Using (22) together with lines 5 and 6 of Algorithm 1, we get

$\displaystyle\hat{{\bf z}}^{k+1}$ $\displaystyle={\bf z}^{k+1}+\mathbf{P}{\bf
m}^{k+1}$ $\displaystyle={\bf z}^{k}+\eta\alpha({\bf z}_{g}^{k}-{\bf
z}^{k})+\Delta^{k}+\mathbf{P}({\bf m}^{k}-\eta\nabla H({\bf
z}_{g}^{k})-\Delta^{k})$ $\displaystyle={\bf z}^{k}+\mathbf{P}{\bf
m}^{k}+\eta\alpha({\bf z}_{g}^{k}-{\bf z}^{k})-\eta\mathbf{P}\nabla H({\bf
z}_{g}^{k})+\Delta^{k}-\mathbf{P}\Delta^{k}.$

From line 4 of Algorithm 1 and $\mathbf{P}\mathbf{W}(k)=\mathbf{W}(k)$ it
follows that $\mathbf{P}\Delta^{k}=\Delta^{k}$, which implies

$\displaystyle\hat{{\bf z}}^{k+1}$ $\displaystyle={\bf z}^{k}+\mathbf{P}{\bf
m}^{k}+\eta\alpha({\bf z}_{g}^{k}-{\bf z}^{k})-\eta\mathbf{P}\nabla H({\bf
z}_{g}^{k})$ $\displaystyle=\hat{{\bf z}}^{k}+\eta\alpha({\bf z}_{g}^{k}-{\bf
z}^{k})-\eta\mathbf{P}\nabla H({\bf z}_{g}^{k}).$

Hence,

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle=\left\|\hat{{\bf z}}^{k}-{\bf z}^{*}+\eta\alpha({\bf
z}_{g}^{k}-{\bf z}^{k})-\eta\mathbf{P}\nabla H({\bf z}_{g}^{k})\right\|^{2}$
$\displaystyle=\left\|(1-\eta\alpha)(\hat{{\bf z}}^{k}-{\bf
z}^{*})+\eta\alpha({\bf z}_{g}^{k}+\mathbf{P}{\bf m}^{k}-{\bf
z}^{*})\right\|^{2}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle-2\eta\langle\mathbf{P}\nabla H({\bf z}_{g}^{k}),{\bf
z}^{k}+\mathbf{P}{\bf m}^{k}-{\bf z}^{*}+\eta\alpha({\bf z}_{g}^{k}-{\bf
z}^{k})\rangle$ $\displaystyle\leq.$

Using inequality
$\left\|a+b\right\|^{2}\leq(1+\gamma)\left\|a\right\|^{2}+(1+\frac{1}{\gamma})\left\|b\right\|^{2},~{}\gamma>0$
with $\gamma=\frac{\eta\alpha}{1-\eta\alpha}$ we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle=(1-\eta\alpha)\left\|\hat{{\bf z}}^{k}-{\bf
z}^{*}\right\|^{2}+\eta\alpha\left\|{\bf z}_{g}^{k}+\mathbf{P}{\bf m}^{k}-{\bf
z}^{*}\right\|^{2}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$ $\displaystyle-2\eta\langle\nabla H({\bf
z}_{g}^{k}),\mathbf{P}({\bf z}_{g}^{k}-{\bf
z}^{*})\rangle+2\eta(1-\eta\alpha)\langle\nabla H({\bf
z}_{g}^{k}),\mathbf{P}({\bf z}_{g}^{k}-{\bf
z}^{k})\rangle-2\eta\langle\mathbf{P}\nabla H({\bf z}_{g}^{k}),{\bf
m}^{k}\rangle$ $\displaystyle\leq(1-\eta\alpha)\left\|\hat{{\bf z}}^{k}-{\bf
z}^{*}\right\|^{2}+2\eta\alpha\left\|{\bf z}_{g}^{k}-{\bf
z}^{*}\right\|^{2}+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$ $\displaystyle-2\eta\langle\nabla H({\bf
z}_{g}^{k}),\mathbf{P}({\bf z}_{g}^{k}-{\bf
z}^{*})\rangle+2\eta(1-\eta\alpha)\langle\nabla H({\bf
z}_{g}^{k}),\mathbf{P}({\bf z}_{g}^{k}-{\bf
z}^{k})\rangle-2\eta\langle\mathbf{P}\nabla H({\bf z}_{g}^{k}),{\bf
m}^{k}\rangle$

One can observe, that ${\bf z}^{k},{\bf z}_{g}^{k},{\bf
z}^{*}\in\operatorname{Im}\mathbf{P}\mathbf{B}$. Hence,

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq(1-\eta\alpha)\left\|\hat{{\bf z}}^{k}-{\bf
z}^{*}\right\|^{2}+2\eta\alpha\left\|{\bf z}_{g}^{k}-{\bf
z}^{*}\right\|^{2}+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$ $\displaystyle-2\eta\langle\nabla H({\bf
z}_{g}^{k}),{\bf z}_{g}^{k}-{\bf
z}^{*}\rangle+2\eta(1-\eta\alpha)\langle\nabla H({\bf z}_{g}^{k}),{\bf
z}_{g}^{k}-{\bf z}^{k}\rangle-2\eta\langle\mathbf{P}\nabla H({\bf
z}_{g}^{k}),{\bf m}^{k}\rangle.$

Using line 3 of Algorithm 1 we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq(1-\eta\alpha)\left\|\hat{{\bf z}}^{k}-{\bf
z}^{*}\right\|^{2}+2\eta\alpha\left\|{\bf z}_{g}^{k}-{\bf
z}^{*}\right\|^{2}+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$ $\displaystyle-2\eta\langle\nabla H({\bf
z}_{g}^{k}),{\bf z}_{g}^{k}-{\bf
z}^{*}\rangle+2\eta(1-\eta\alpha)\frac{(1-\tau)}{\tau}\langle\nabla H({\bf
z}_{g}^{k}),{\bf z}_{f}^{k}-{\bf
z}_{g}^{k}\rangle-2\eta\langle\mathbf{P}\nabla H({\bf z}_{g}^{k}),{\bf
m}^{k}\rangle.$

Using convexity and $\mu_{H}$-strong convexity of $H({\bf z})$ on
$\operatorname{Im}\mathbf{P}\mathbf{B}$ we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq(1-\eta\alpha)\left\|\hat{{\bf z}}^{k}-{\bf
z}^{*}\right\|^{2}+2\eta\alpha\left\|{\bf z}_{g}^{k}-{\bf
z}^{*}\right\|^{2}+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$ $\displaystyle-2\eta(H({\bf
z}_{g}^{k})-H({\bf z}^{*}))-\eta\mu_{H}\left\|{\bf z}_{g}^{k}-{\bf
z}^{*}\right\|^{2}$
$\displaystyle+2\eta(1-\eta\alpha)\frac{(1-\tau)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}_{g}^{k}))-2\eta\langle\mathbf{P}\nabla H({\bf
z}_{g}^{k}),{\bf m}^{k}\rangle$ $\displaystyle=(1-\eta\alpha)\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\left(2\eta\alpha-\eta\mu_{H}\right)\left\|{\bf
z}_{g}^{k}-{\bf z}^{*}\right\|^{2}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$ $\displaystyle-2\eta(H({\bf
z}_{g}^{k})-H({\bf z}^{*}))+2\eta(1-\eta\alpha)\frac{(1-\tau)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}_{g}^{k}))$ $\displaystyle-2\eta\langle\mathbf{P}\nabla
H({\bf z}_{g}^{k}),{\bf m}^{k}\rangle+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}.$

Using $\alpha$ defined by (16) we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf z}^{*}\right\|^{2}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$ $\displaystyle-2\eta(H({\bf
z}_{g}^{k})-H({\bf z}^{*}))+2\eta(1-\eta\alpha)\frac{(1-\tau)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}_{g}^{k}))$ $\displaystyle-2\eta\langle\mathbf{P}\nabla
H({\bf z}_{g}^{k}),{\bf m}^{k}\rangle+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}.$

Since $H({\bf z}_{g}^{k})\geq H({\bf z}^{*})$, we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf z}^{*}\right\|^{2}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle-2\eta(1-\eta\alpha)(H({\bf z}_{g}^{k})-H({\bf
z}^{*}))+2\eta(1-\eta\alpha)\frac{(1-\tau)}{\tau}(H({\bf z}_{f}^{k})-H({\bf
z}_{g}^{k}))$ $\displaystyle-2\eta\langle\mathbf{P}\nabla H({\bf
z}_{g}^{k}),{\bf m}^{k}\rangle+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}$
$\displaystyle=\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf z}^{*}\right\|^{2}+\eta^{2}\left\|\nabla H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle+2\eta(1-\eta\alpha)\left(\frac{(1-\tau)}{\tau}H({\bf
z}_{f}^{k})+H({\bf z}^{*})-\frac{1}{\tau}H({\bf z}_{g}^{k})\right)$
$\displaystyle-2\eta\langle\mathbf{P}\nabla H({\bf z}_{g}^{k}),{\bf
m}^{k}\rangle+2\eta\alpha\left\|{\bf m}^{k}\right\|^{2}_{\mathbf{P}}.$

Using (14) and $\theta$ defined by (18) we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\left(\eta^{2}-\frac{(1-\eta\alpha)\eta\lambda_{\min}^{+}}{\tau\lambda_{\max}L_{H}}\right)\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle+(1-\tau)\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))-\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k+1})-H({\bf z}^{*}))$ $\displaystyle-2\eta\langle\mathbf{P}\nabla
H({\bf z}_{g}^{k}),{\bf m}^{k}\rangle+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}.$

Using Young’s inequality we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\left(\eta^{2}-\frac{(1-\eta\alpha)\eta\lambda_{\min}^{+}}{\tau\lambda_{\max}L_{H}}\right)\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle+(1-\tau)\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))-\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k+1})-H({\bf z}^{*}))$
$\displaystyle+\frac{\eta^{2}\lambda_{\max}}{\lambda_{\min}^{+}}\left\|\nabla
H({\bf
z}_{g}^{k})\right\|^{2}_{\mathbf{P}}+\frac{\lambda_{\min}^{+}}{\lambda_{\max}}\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+2\eta\alpha\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}$
$\displaystyle=\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\left(\eta^{2}+\frac{\eta^{2}\lambda_{\max}}{\lambda_{\min}^{+}}-\frac{(1-\eta\alpha)\eta\lambda_{\min}^{+}}{\tau\lambda_{\max}L_{H}}\right)\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle+(1-\tau)\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))-\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k+1})-H({\bf z}^{*}))$
$\displaystyle+\left(\frac{\lambda_{\min}^{+}}{\lambda_{\max}}+2\eta\alpha\right)\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}.$

Using (17) and (16), that imply
$\eta\alpha\leq\frac{\lambda_{\min}^{+}}{4\lambda_{\max}}$, we obtain

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\left(\eta^{2}+\frac{\eta^{2}\lambda_{\max}}{\lambda_{\min}^{+}}-\frac{3\eta\lambda_{\min}^{+}}{4\tau\lambda_{\max}L_{H}}\right)\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle+(1-\tau)\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))-\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k+1})-H({\bf
z}^{*}))+\frac{3\lambda_{\min}^{+}}{2\lambda_{\max}}\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}.$

Using (15) and $\sigma$ defined by (19) we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\left(\eta^{2}+\frac{\eta^{2}\lambda_{\max}}{\lambda_{\min}^{+}}-\frac{3\eta\lambda_{\min}^{+}}{4\tau\lambda_{\max}L_{H}}\right)\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle+(1-\tau)\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))-\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k+1})-H({\bf z}^{*}))$
$\displaystyle+\left(1-\frac{\lambda_{\min}^{+}}{4\lambda_{\max}}\right)6\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}-6\left\|{\bf
m}^{k+1}\right\|^{2}_{\mathbf{P}}+\frac{12\eta^{2}\lambda_{\max}}{\lambda_{\min}^{+}}\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle\leq\left(1-\frac{\eta\mu_{H}}{2}\right)\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\left(\frac{14\eta^{2}\lambda_{\max}}{\lambda_{\min}^{+}}-\frac{3\eta\lambda_{\min}^{+}}{4\tau\lambda_{\max}L_{H}}\right)\left\|\nabla
H({\bf z}_{g}^{k})\right\|^{2}_{\mathbf{P}}$
$\displaystyle+(1-\tau)\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))-\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k+1})-H({\bf z}^{*}))$
$\displaystyle+\left(1-\frac{\lambda_{\min}^{+}}{4\lambda_{\max}}\right)6\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}-6\left\|{\bf
m}^{k+1}\right\|^{2}_{\mathbf{P}}.$

Using $\eta$ defined by (17) and $\tau$ defined by (20) we get

$\displaystyle\left\|\hat{{\bf z}}^{k+1}-{\bf z}^{*}\right\|^{2}$
$\displaystyle\leq\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu_{H}}{L_{H}}}\right)\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\left(1-\frac{\lambda_{\min}^{+}}{4\lambda_{\max}}\right)6\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}-6\left\|{\bf
m}^{k+1}\right\|^{2}_{\mathbf{P}}$
$\displaystyle+\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu_{H}}{L_{H}}}\right)\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))-\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k+1})-H({\bf z}^{*}))$
$\displaystyle\leq\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu_{H}}{L_{H}}}\right)\left(\left\|\hat{{\bf
z}}^{k}-{\bf z}^{*}\right\|^{2}+\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))+6\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}\right)$
$\displaystyle-\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf z}_{f}^{k+1})-H({\bf
z}^{*}))-6\left\|{\bf m}^{k+1}\right\|^{2}_{\mathbf{P}}.$

Rearranging and using (21) concludes the proof. ∎

###  9.4 Proof of Theorem 3

###### Proof.

From derivation of the reformulated problem and Demyanov-Danskin theorem it
follows that $\nabla F^{*}(\mathbf{B}^{\top}{\bf z}^{*})={\bf x}^{*}$.
Therefore $\nabla H({\bf z}^{*})=(0_{d},{\bf x}^{*})^{\top}$. Using
$L_{H}$-smoothness of $H$ on $\operatorname{Im}\mathbf{P}\mathbf{B}$ we get

$\displaystyle\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf z}_{g}^{k})-{\bf
x}^{*}\right\|^{2}$ $\displaystyle=\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf
z}_{g}^{k})-F^{*}(\mathbf{B}^{\top}{\bf z}^{*})\right\|^{2}\leq\left\|\nabla
H({\bf z}_{g}^{k})-\nabla H({\bf z}^{*})\right\|^{2}\leq L_{H}^{2}\left\|{\bf
z}_{g}^{k}-{\bf z}^{*}\right\|^{2}.$

Using line 3 of Algorithm 1 and inequality
$\left\|a+b\right\|^{2}\leq(1+\gamma)\left\|a\right\|^{2}+(1+\frac{1}{\gamma})\left\|b\right\|^{2},~{}\gamma>0$
with $\gamma=\frac{1}{\tau}-1$ we get we get

$\displaystyle\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf z}_{g}^{k})-{\bf
x}^{*}\right\|^{2}$ $\displaystyle\leq\tau L_{H}^{2}\left\|{\bf z}^{k}-{\bf
z}^{*}\right\|^{2}+(1-\tau)L_{H}^{2}\left\|{\bf z}_{f}^{k}-{\bf
z}^{*}\right\|^{2}.$

Using $\mu_{H}$-strong convexity of $H$ on
$\operatorname{Im}\mathbf{P}\mathbf{B}$ we get

$\displaystyle\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf z}_{g}^{k})-{\bf
x}^{*}\right\|^{2}$ $\displaystyle\leq\tau L_{H}^{2}\left\|{\bf z}^{k}-{\bf
z}^{*}\right\|^{2}+\frac{2(1-\tau)L_{H}^{2}}{\mu_{H}}(H({\bf
z}_{f}^{k})-H({\bf z}^{*})).$

Using (22) we get

$\displaystyle\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf z}_{g}^{k})-{\bf
x}^{*}\right\|^{2}$ $\displaystyle\leq 2\tau L_{H}^{2}\left\|\hat{{\bf
z}}^{k}-{\bf z}^{*}\right\|^{2}+2\tau L_{H}^{2}\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}+\frac{2(1-\tau)L_{H}^{2}}{\mu_{H}}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))$ $\displaystyle=2\tau L_{H}^{2}\left\|\hat{{\bf
z}}^{k}-{\bf
z}^{*}\right\|^{2}+\frac{\tau(1-\tau)L_{H}^{2}}{\eta(1-\eta\alpha)\mu_{H}}\frac{2\eta(1-\eta\alpha)}{\tau}(H({\bf
z}_{f}^{k})-H({\bf z}^{*}))+\frac{\tau L_{H}^{2}}{3}6\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}.$ $\displaystyle\leq\max\left\\{2\tau
L_{H}^{2},\frac{\tau(1-\tau)L_{H}^{2}}{\eta(1-\eta\alpha)\mu_{H}},\frac{\tau
L_{H}^{2}}{3}\right\\}\left(\left\|\hat{{\bf z}}^{k}-{\bf
z}^{*}\right\|^{2}+\frac{2\eta(1-\eta\alpha)}{\tau}(F^{*}({\bf
z}_{f}^{k})-F^{*}({\bf z}^{*}))+6\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}\right)$ $\displaystyle=\max\left\\{2\tau
L_{H}^{2},\frac{\tau(1-\tau)L_{H}^{2}}{\eta(1-\eta\alpha)\mu_{H}}\right\\}\left(\left\|\hat{{\bf
z}}^{k}-{\bf z}^{*}\right\|^{2}+\frac{2\eta(1-\eta\alpha)}{\tau}(F^{*}({\bf
z}_{f}^{k})-F^{*}({\bf z}^{*}))+6\left\|{\bf
m}^{k}\right\|^{2}_{\mathbf{P}}\right).$

Using the definition of $\Psi^{k}$ (21) and denoting
$C=\Psi^{0}\max\left\\{2\tau
L_{H}^{2},\frac{\tau(1-\tau)L_{H}^{2}}{\eta(1-\eta\alpha)\mu_{H}}\right\\}$ we
get

$\displaystyle\left\|\nabla F^{*}(\mathbf{B}^{\top}{\bf z}_{g}^{k})-{\bf
x}^{*}\right\|^{2}\leq\frac{C}{\Psi^{0}}\Psi^{k}.$

Applying Lemma 3 concludes the proof. ∎

## References

  *   * Nedic et al. [2010] Nedic, A., Ozdaglar, A., Parrilo, P.A.: Constrained consensus and optimization in multi-agent networks. IEEE Transactions on Automatic Control 55(4), 922–938 (2010)
  * Necoara et al. [2011] Necoara, I., Nedelcu, V., Dumitrache, I.: Parallel and distributed optimization methods for estimation and control in networks. Journal of Process Control 21(5), 756–766 (2011)
  * Scutari et al. [2016] Scutari, G., Facchinei, F., Lampariello, L.: Parallel and distributed methods for constrained nonconvex optimization—part i: Theory. IEEE Transactions on Signal Processing 65(8), 1929–1944 (2016)
  * Scutari and Sun [2019] Scutari, G., Sun, Y.: Distributed nonconvex constrained optimization over time-varying digraphs. Mathematical Programming 176(1), 497–544 (2019)
  * Wang and Hu [2022] Wang, J., Hu, G.: Distributed optimization with coupling constraints in multi-cluster networks based on dual proximal gradient method. arXiv preprint arXiv:2203.00956 (2022)
  * Wu et al. [2022] Wu, X., Wang, H., Lu, J.: Distributed optimization with coupling constraints. IEEE Transactions on Automatic Control (2022)
  * Zhu and Martinez [2011] Zhu, M., Martinez, S.: On distributed convex optimization under inequality and equality constraints. IEEE Transactions on Automatic Control 57(1), 151–164 (2011)
  * [8] Gong, K., Zhang, L.: Push-pull based distributed primal-dual algorithm for coupled constrained convex optimization in multi-agent networks. Available at SSRN 4109852
  * Liang et al. [2019] Liang, S., Yin, G., et al.: Distributed smooth convex optimization with coupled constraints. IEEE Transactions on Automatic Control 65(1), 347–353 (2019)
  * Alghunaim et al. [2018] Alghunaim, S.A., Yuan, K., Sayed, A.H.: Dual coupled diffusion for distributed optimization with affine constraints. In: 2018 IEEE Conference on Decision and Control (CDC), pp. 829–834 (2018). IEEE
  * Carli and Dotoli [2019] Carli, R., Dotoli, M.: Distributed alternating direction method of multipliers for linearly constrained optimization over a network. IEEE Control Systems Letters 4(1), 247–252 (2019)
  * Aybat and Hamedani [2019] Aybat, N.S., Hamedani, E.Y.: A distributed admm-like method for resource sharing over time-varying networks. SIAM Journal on Optimization 29(4), 3036–3068 (2019)
  * Chang [2016] Chang, T.-H.: A proximal dual consensus admm method for multi-agent constrained optimization. IEEE Transactions on Signal Processing 64(14), 3719–3734 (2016)
  * Kovalev et al. [2021a] Kovalev, D., Shulgin, E., Richtárik, P., Rogozin, A., Gasnikov, A.: Adom: Accelerated decentralized optimization method for time-varying networks. arXiv preprint arXiv:2102.09234 (2021)
  * Kovalev et al. [2021b] Kovalev, D., Gasanov, E., Gasnikov, A., Richtarik, P.: Lower bounds and optimal algorithms for smooth and strongly convex decentralized optimization over time-varying networks. Advances in Neural Information Processing Systems 34 (2021)
  * Rogozin et al. [2022] Rogozin, A., Yarmoshik, D., Kopylova, K., Gasnikov, A.: Decentralized strongly-convex optimization with affine constraints: Primal and dual approaches. arXiv preprint arXiv:2207.04555 (2022)
  * Zhou and Lange [2013] Zhou, H., Lange, K.: A path algorithm for constrained estimation. Journal of Computational and Graphical Statistics 22(2), 261–283 (2013)
  * Nesterov [2004] Nesterov, Y.: Introductory Lectures on Convex Optimization: a Basic Course. Kluwer Academic Publishers, Massachusetts, ??? (2004)
  * Salim et al. [2022] Salim, A., Condat, L., Kovalev, D., Richtárik, P.: An optimal algorithm for strongly convex minimization under affine constraints. In: International Conference on Artificial Intelligence and Statistics, pp. 4482–4498 (2022). PMLR
  * Scaman et al. [2017] Scaman, K., Bach, F., Bubeck, S., Lee, Y.T., Massoulié, L.: Optimal algorithms for smooth and strongly convex distributed optimization in networks. In: Proceedings of the 34th International Conference on Machine Learning-Volume 70, pp. 3027–3036 (2017). JMLR. org
  * Yarmoshik et al. [2022] Yarmoshik, D., Rogozin, A., Khamisov, O., Dvurechensky, P., Gasnikov, A., et al.: Decentralized convex optimization under affine constraints for power systems control. arXiv preprint arXiv:2203.16686 (2022)


# Learning Models for Query by Vocal Percussion: A Comparative Study

###### Abstract

The imitation of percussive sounds via the human voice is a natural and
effective tool for communicating rhythmic ideas on the fly. Thus, the
automatic retrieval of drum sounds using vocal percussion can help artists
prototype drum patterns in a comfortable and quick way, smoothing the creative
workflow as a result. Here we explore different strategies to perform this
type of query, making use of both traditional machine learning algorithms and
recent deep learning techniques. The main hyperparameters from the models
involved are carefully selected by feeding performance metrics to a grid
search algorithm. We also look into several audio data augmentation
techniques, which can potentially regularise deep learning models and improve
generalisation. We compare the final performances in terms of effectiveness
(classification accuracy), efficiency (computational speed), stability
(performance consistency), and interpretability (decision patterns), and
discuss the relevance of these results when it comes to the design of
successful query-by-vocal-percussion systems.

##  1 Introduction

One of the aims of Music Information Retrieval (MIR) as a field of research is
to provide musicians and music producers with new tools that can make their
workflow more pleasant. Thanks to this, mechanical and sometimes laborious
tasks can now be fully automated, allowing artists to better focus on the most
creative aspects of their craft. Subfields of MIR like musical key detection,
chord estimation, and automatic mixing are dedicated to the development of
these methods.

Another subfield that shares the same goal is Query by Example (QE). It
investigates how to quickly query a specific audio file from a sound library
by providing another file that sounds like it, like a sketch recording. The
searching process is usually fully audio content-based, and therefore does not
require any other kind of metadata like labels or text descriptions. A special
case of QE is that of Query by Vocal Imitation (QVI), in which the example
audio file fed to the system contains a vocal imitation of the desired sound.
The subject covered in this study, Query by Vocal Percussion (QVP), is a form
of QVI in which percussive sounds are retrieved from percussive vocal
imitations. In particular, we focus on retrieving four types of drum sounds:
kick drum, snare drum, closed hi-hat, and opened hi-hat.

Querying drum sounds by vocal percussion requires distinct approaches given
the context and resources available. For instance, a basic “naive” algorithm
could simply retrieve the most relevant drum samples by objective sound
content similarity with the vocal imitation (e.g. by spectral distance) [1].
Another more nuanced case-study is when one wants to retrieve the sound of a
snare drum from a particular brand or articulated in a specific manner [2].
Here, the algorithms must be aware of both the different sounds that snare
drums can make and how differently the user vocalises them. This is,
therefore, a problem of correspondence between two sounds, and the main
objective is to find the set of audio features that best link them. Lastly, if
one wants to trigger a fixed drum sample when a specific vocal imitation is
given, then the algorithms are just required to learn how to distinguish
between vocal imitations alone, meant to trigger different drum samples [3].
This is a problem of correspondence between a sound and a label
(classification), and the relevant set of audio features would now be the one
that best separates all classes. The focus of the present work is put on this
latter task.

A literature review of relevant past work on QVP is given in section 2. The
methods and sub-routines to study are laid out in section 3 and we give a
detailed summary of the final results in section 4. We conclude with a brief
recapitulation and future perspectives in section 5.

##  2 Problem Definition and Past Work

An important thing to note about our approach to performing QVP is that, when
considering a realistic scenario, it is not entirely reducible to a regular
sound event classification task. Here the users are not supposed to label the
utterances they record for the algorithm to train successfully, so the labels
have to be made implicit in some way. This is the reason why the final
training and testing datasets end up being of a slightly different nature.

For the training dataset, the users of the QVP system record several
utterances of each drum type in four separated audio files (kick drum, snare
drum, closed hi-hat, and opened hi-hat). This allows the algorithm to
automatically assign a different class to sounds depending on the file that
contains them. In this way, all the utterances are automatically labelled
without requiring the users to do it themselves. For the testing dataset, the
users record a beatbox-style improvisation with utterances of mixed classes in
a single file, without label annotations. The trained system would try to
correctly classify the vocalisations despite the expressive resources in the
improvised performance, like modulations in pitch, loudness, or duration. The
challenge is, then, to extract the most descriptive audio features possible in
the training process so that the algorithm is robust enough to handle these
expressive nuances in the testing data.

Most relevant past work on QVP is concerned with utterance classification
rather than querying. This means that only one of the two datasets (isolated
utterances [4] or full improvisation [5]) is used for both training and
testing, generally resulting in higher classification accuracies111We obtain
98% classification accuracy on this task using our dataset. The works most
often feature three classes of imitated sounds (kick drum, snare drum, and hi-
hat) instead of four. While automatic onset detection is explored as well in
some of these projects, our study is exclusively focused on the classification
of vocal percussion utterances into their correct drum types.

In [4], Hazan and Ramirez feed spectral and temporal audio features from the
vocal percussion utterances to a decision tree algorithm, reporting results of
up to 90% classification accuracy in testing sets. Sinyor et al. also explored
the AdaBoost algorithm in [6], achieving the best classification results (98%)
together with a decision tree. Nakano et al. used Mel Frequency Cepstral
Coefficients and a Hidden Markov Model to classify drum patterns containing
kick and snare drums in [7], hitting 93% accuracy. Kapur et al. [8] fed
different types of audio low-level descriptors to an artificial neural network
and investigated how well each of them could perform the task on its own. They
found that the zero-crossing rate descriptor was the most informative,
reaching an accuracy of 97%. Stowell et al. [5] observed that one could obtain
higher classification accuracies in beatbox utterance classification by
delaying the start of the first analysis frame 23 milliseconds from each
event’s onset, ignoring the information in the sound’s transient. Finally, in
[9], Ramires et al. used a set of engineered features and a K-Nearest
Neighbours algorithm to build a user-based vocal percussion classifier that
consistently performed better than a general-purpose one. They also introduce
a QVP plugin for Digital Audio Workstations in [3].

We draw inspiration from the studies above when selecting audio descriptors
and machine learning algorithms to classify the utterances. To our knowledge,
this is the first study that introduces data augmentation techniques and deep
learning strategies for QVP.

##  3 Methodology

Figure 1: Diagram illustrating the QVP workflow followed in this study. CNN
stands for Convolutional Neural Network, TCNN for Triplet Convolutional Neural
Network, MLP for Multi-Layer Perceptron, SLP for Single-Layer Perceptron and
ML for other Machine Learning models.

###  3.1 Dataset

The Amateur Vocal Percussion (AVP) [10] dataset contains vocal percussion
recordings from people with little or no experience in beatboxing. It includes
9780 utterances imitating four drum types (kick drum, snare drum, closed hi-
hat, and opened hi-hat) recorded by 28 participants and with annotated onsets
and labels. We use the personal subset, where the participants are asked to
imitate the drums in their own way. For every participant, this subset is
composed of five audio files: four files containing around twenty-five
isolated utterances from the four classes (train set) and one file of
arbitrary length with a beatbox-style improvisation using the chosen personal
imitations (test set). The data belonging to the three discarded participants
is included here for evaluation, as we also want to see how the models decide
between labels when given ambiguous utterances that are hard to classify.

###  3.2 Feature Extraction and Selection

Instead of feeding the raw waveforms directly to the machine learning models,
a good practice in MIR is to compress the information into lower-dimensional
feature vectors to make the analysis more manageable.

One of these arrays is the mel frequency spectrogram, a time-frequency
representation that compresses the frequency channels of the Fourier
spectrogram into a small number of bands by applying a transformation based on
the mel scale. We feed mel spectrograms from single utterances to the
convolutional models in sections 3.4.2 and 3.4.3. After fixing the length of
all audio waveforms to 0.68 seconds (30k samples at 44100 hz), we calculate
the mel spectrogram using a window size of 96 milliseconds and a hop size of
$floor(1000\cdot\frac{30}{N_{b}})$, where $N_{b}$ is the number of mel bands
(8, 12 or 16). By making the hop size dependent on the number of mel bands
that way, we obtain square mel spectrogram representations (8x8, 12x12 or
16x16). This squareness property allows us to explore different hop sizes and
frequency resolutions while making both data augmentation filtering and the
design of network architectures more comfortable.

We also study another type of feature set composed of engineered features.
These are made out of audio low-level descriptors (e.g. spectral centroid) and
statistical functionals (e.g. mean value), used to aggregate these descriptors
through time. We feed these engineered features to the multi-layer perceptron
in 3.4.1 and the rest of machine learning models in 3.4.4. For each vocal
percussion utterance, we compute 13 Mel Frequency Cepstral Coefficients
(MFCCs), the zero-crossing rate, and several spectral descriptors (centroid,
spread, skewness, kurtosis, flatness, and roll-off). We extract all of these
using Librosa [11]. In the case of statistical functionals, we compute the
mean, standard deviation, minimum, maximum, and interquartile range of each
descriptor through time. We have, therefore, 20 descriptors and 5 functionals,
making a set of 100 engineered features.

Finally, before feeding these engineered features to the final classifiers, we
use random forests to perform feature selection in order to minimise
overfitting and maximise accuracy. This way, random feature subsets are fed to
decision trees that perform the utterance classification task in the usual way
and then derive feature importances based on how many times they appear in
these trees. We implement the random forest algorithm using the Scikit-Learn
library [12]. We also study the effect of selecting different amounts of
features (1, 2, 4, 8, 16 and 32) to make the final set.

###  3.3 Data Augmentation

The scarcity of training data is usually a matter of concern if this is high-
dimensional. The curse of dimensionality specifies that when our feature space
is high-dimensional (e.g. waveform, spectrogram…) also is the solution space
that the machine learning model is operating in. This causes sparsity in the
data, making the model unable to obtain results with enough statistical
reliability.

In our case, our training sets have only around one hundred utterances each,
which is not an appropriate amount of data to train high-dimensional
parametric models like neural networks. In order to tackle this problem, we
use a regularisation technique known as data augmentation, which applies
realistic transformations to the original samples (e.g. a small amount of
pitch-shifting) so to create new files and include them in the final dataset.

Two different data augmentation techniques are explored in this study:
waveform data augmentation and spectrogram data augmentation.

Waveform Data Augmentation is applied to the raw audio waveform of each file.
We use three different methods: pitch-shifting (-1 and 1 semitones), as vocal
percussion utterances sometimes include vowels with harmonic content; time-
stretching (0.85 and 1.15), as utterances from the same class usually vary in
length; and noise addition (0.005 and 0.01 amplitude factor from pink noise),
which gently distorts the signal to let algorithms focus on the core timbral
characteristics of the utterances. We use Librosa [11] to make these
transformations.

Spectrogram Data Augmentation, which is based on the SpecAugment method for
speech recognition [13], is applied to the spectrogram representation of each
file and is composed by three different techniques: frequency channel zero-
masking, time step zero-masking and time warping. With $N_{b}$ being the
number of mel bands of the original spectrogram, we apply a random number of
masks to each file in the range $[1,floor(2\cdot\frac{N_{b}}{8})]$, each of
one of random width in the range $[1,ceil(\frac{N_{b}}{3})]$.

We experiment with serial combinations of these data augmentation techniques
so to create expanded datasets that could improve performance accuracy for
some models.

###  3.4 Models

####  3.4.1 Multi-Layer Perceptron

The Multi-Layer Perceptron (MLP) [14] is the most simple feedforward neural
network model, as well as the oldest. It has at least three layers of
artificial neurons: an input layer with the extracted features, a hidden layer
with learnt representations, and an output layer with the labels. These
neurons are fully connected by trainable parameters called weights in a
sequential way (input-hidden-output) and nonlinear activation functions are
applied to neurons in the hidden and output layers, increasing the complexity
of the learnt representation. Like most neural network models, MLPs use an
optimization technique called gradient descent that is made possible using the
backpropagation algorithm [15].

We feed engineered features to MLPs. The number of neurons in the hidden
layers that we explore for these network architectures are
[$\frac{n_{i}}{2}$], [$\frac{n_{i}}{2}$,$\frac{n_{i}}{4}$] and
[$\frac{n_{i}}{4}$], where $n_{i}$ is the number of input (selected) features.

####  3.4.2 End-To-End Convolutional Neural Network

Convolutional Neural Networks (CNN) [16] are another type of feed-forward
networks that include convolutional layers in their architecture. These
convolutional layers are composed of several small local filter kernels that
convolve with the input data, generating learnt feature maps as a result. The
usual method to build end-to-end CNNs is to stack various convolutional layers
to filter the data and create feature maps, pooling layers to subsample these
maps, and fully connected layers (MLPs) to classify the samples.

We feed the two-dimensional mel spectrograms to different CNN architectures.
The number of convolutional layers for these architectures are either 2, 4, or
6, and we also explore different number of filters for each of the layers (8,
16, 32, and 64). The kernel size of all convolutional filters is 3x3 with
stride 1, while the kernel size corresponding to the max-pooling layers is 2x2
with stride 2, which downsamples the feature maps by a factor of 2. We
experiment with different numbers of pooling layers (1, 2, and 3) and use
batch normalization as an extra regularisation step.

####  3.4.3 Triplet Convolutional Neural Network

A Triplet Convolutional Neural Network (TCNN) [17] is a CNN that uses a
special type of contrasting loss called the triplet loss to make data
embeddings (learnt features) that maximise separability between classes. It
is, therefore, a learning-based feature extractor and not a classifier. This
loss works by forming random data triplets with data samples: an anchor, a
positive and a negative sample. The anchor and the positive samples are of the
same class while the negative is not, and the goal of the algorithm is to
minimise the distance (usually Euclidean) between the embeddings from the
anchor and the positive samples while maximising the distance between the
embeddings from the anchor and the negative samples. In summary, the algorithm
tries to minimise the following function:

$\mathcal{L}=\max(d(a,p)-d(a,n)+\textit{margin},0)$ (1)

, where $a$, $p$ and $n$ are the embeddings corresponding to the anchor,
positive and negative samples respectively, $d(x,y)$ is a distance function
applied to two arrays of embeddings and margin is a fixed hyperparameter that
regulates the target closeness and separation between the embeddings. We apply
a Euclidean distance function and set the margin to 0 for the calculation of
this loss.

After calculating the losses for all triplets in a training batch using
equation 1, we need to make a final decision on which of these triplet losses
we want to pass to the network for backpropagation [18]. For instance, we may
only want to pass the greatest loss value, the one belonging to the hardest
triplet, or maybe pass the average value of all the triplet losses calculated
for that specific batch. Two triplet selectors are explored for this purpose:
random negative and hardest negative. We avoid averaging an elevated number of
triplet losses, as this is known to cause instability in the training process
[18].

Once the embeddings have been calculated, a Single-Layer Perceptron (SLP)
classifier, composed of only an input and an output layer, is trained
separately to extract the final performance accuracies. We build the networks
with the same architectures used for the CNNs (section 3.4.2) except for the
output layer, which has 16 or 32 neurons (embeddings) instead of 4 (classes).

####  3.4.4 Other Machine Learning Models

A final set of traditional machine learning methods is used for classification
using engineered features. It includes the following algorithms: K-Nearest
Neighbours, Support Vector Machine (with linear and radial basis function
kernels), Decision Tree, Random Forest, AdaBoost, Naive Bayes, Quadratic
Discriminant Analysis, and Single-Layer Perceptron. The values of the main
hyperparameters that give each algorithm their best performance (e.g. number
of tree estimators in the random forest) are found using grid search.

For this cluster of algorithms, the model that performs best is the one taken
for evaluation in section 4.

###  3.5 Training Methodology

After standardising the input features, the models described in the last
section are trained using different techniques. For the neural networks (MLP,
CNN, TCNN), several hyperparameters and routines must be chosen. These
parameters are the validation split, which determines the percentage of
samples used for network validation; the maximum number of epochs that the
network is allowed to reach during training; the optimiser type, which makes
the network gradually lower its performance error; the loss type, which is the
error metric used to guide the optimiser; the learning rate, which dictates
how fast the network is learning222If the learning rate is too high, it causes
instability on the network’s performance. If it is too low, it might take a
long time to train the network.; and the batch size, which indicates the
number of samples from the training set used to calculate a single loss value.

For all neural network models, we choose a validation split of 10%, a maximum
number of epochs of 200, and an Adam optimiser. We also investigate the effect
of different learning rates on the training of the models, with eleven
logarithmically spaced values in the range [0.1,0.000001]. We make sure that
the chosen batch size is not larger than the tenth part of the training set.
We assign a batch size of 8 to the original (non-augmented) datasets and
experiment with sizes of 128, 512, and 1024 for augmented datasets. We choose
cross-entropy as the loss function for both the MLP and the CNN (TCNN already
includes the triplet loss). We choose all these parameters and parameter
ranges in light of good practices for CNNs [18].

The networks also have two more sets of parameters that control the training
process through the validation loss. These are the learning rate scheduler and
the early-stopping engine. The former is concerned with reducing the learning
rate by a factor $f_{l}$ for subsequent epochs if the validation error has not
decreased in the last $N_{l}$ epochs. The latter automatically stops the
training process if the validation error has not decreased in the last $N_{e}$
epochs. For all networks, we set $N_{l}$ to 0.1. For the MLP, we set $f_{l}$
to 2 and $N_{e}$ to 3; For the CNN, $f_{l}$ is 3 and $N_{e}$ is 8. Finally,
for the TCNN, we have a $f_{l}$ of 2 and a $N_{e}$ of 3.

##  4 Results

###  4.1 Effectiveness

Final classification accuracies are gathered in table 1. Random Forests (RF)
with a hundred tree estimators were the best performing classifiers of all
machine learning models in 3.4.4. The values shown are the percentages of
correctly classified utterances in the test set (improvisation) for all 28
participants.

We observe that neural network models performed best using the augmented set.
They also show a preference for medium learning rates between 0.001 and
0.0001, which is in line with the observations in [18]. Best results for CNNs
are collected using spectrogram representations of sizes 12x12 and 16x16,
while the best results for TCNNs are obtained with representations of size
8x8. Data augmentation seems to not have any major effect on the performance
of machine learning models. Also, triplet models perform best when the loss
from the hardest negative triplet (see 3.4.3) is passed to the backpropagation
algorithm.

| MLP | CNN | TCNN | RF
---|---|---|---|---
Original Dataset | 74.3 | 77.7 | 77.5 | 77.6
Augmented Dataset | 76.9 | 82.2 | 80.9 | 78.3
Table 1: QVP classification accuracies (%) of Multi-Layer Perceptron (MLP),
Convolutional Neural Network (CNN), Triplet Convolutional Neural Network
(TCNN) and Random Forests (RF).

###  4.2 Efficiency

We choose the algorithms that gave the best results for each of the models to
evaluate performance efficiency (speed). In particular, we select the ones
whose accuracies lie between the top accuracy and one percentage point below
it. We include all preprocessing stages in the measurements (data
augmentation, feature selection…). Training and testing times are averaged
across the 28 participants, as model parameters and dataset lengths gave rise
to very different training and testing times depending on the participant. All
experiments are performed in the same machine (MacBook Pro 2015 laptop) and
under the same conditions.

| MLP | CNN | TCNN | RF
---|---|---|---|---
Training Speed | 14.08 | 217.42 | 28.83 | 7.79
Testing Speed | 1.19 | 1.46 | 1.80 | 1.44
Table 2: Performace speeds for each type of model in seconds. Values are
averaged across all participants.

Results are shown in table 2. Both the CNN and the TCNN models take augmented
datasets as input, while the MLP and the random forest (RF) use the original
dataset. We see that the CNN model, while being the most accurate one, it
trains very slowly. Its testing speed, however, is comparable with those of
the rest of the models. Hence, if waiting for training to finish is not an
issue, the user of the QVP system could use an end-to-end CNN algorithm to
classify utterances as fast as the other methods and with higher accuracy.

###  4.3 Stability

In order to evaluate performance stability, we take a look at how each model
performs when given different combinations of hyperparameter values. This
would be the first metric of performance stability. Then, once discovered the
optimal hyperparameter values for each model (giving the accuracies in table
1), we study how these specific algorithms perform when their internal
parameters (e.g. decision tree’s random state, neural networks’ weights… etc.)
are given by random initialisation. This would be the second metric of
performance stability.

The probability density functions of both measures of stability for all four
models are shown in figure 2. The shape of these distributions reveals whether
the results of a particular model are strongly tied to the choice of
hyperparameters or to the random initialisation of its internal parameters. In
this way, we see that the most stable models to hyperparameter values are the
two convolutional neural networks (CNN and TCNN), whose probability density
functions showcase the least amount of statistical variance. When it comes to
random initialisation of internal parameters, however, we see that random
forests are more stable than neural networks. This could be explained by
considering random forests as a type of ensemble learning (various decision
trees deciding the most probable label), making them more prone to stability
than the rest of the algorithms, which stand on their own.

Figure 2: Probability density functions to evaluate performance stability of
the four models: MLP in red, CNN in blue, TCNN in black and RF in green.

###  4.4 Interpretability

While it is challenging to interpret neural networks’ decisions, models that
take engineered features as inputs are generally more scrutable. For example,
we can observe the kind of audio descriptors that were derived from the
feature selection process (see 3.2) so to then relate them to psychoacoustic
phenomena that we can understand better. The most relevant selected descriptor
across participants and models was the spectral centroid, which is often
described as the “centre of mass” of the frequency axis of the spectrum and
happens to be one of the most important descriptors in sound timbre analysis.
It was closely followed by the first MFCC and the third MFCC, which are also
related to the utterance’s low to high frequency ratio. For some participants,
higher MFCCs were included in the top selected features, which encode more
complex spectral information. They are possibly used to separate vocal
imitations of closed hi-hat and opened hi-hat, which are often close in
timbral similarity.

Related to the last point, the most common mistake for all four models is the
confusion between closed hi-hat and opened hi-hat vocalisations. Many
participants imitate them in a similar manner and, while the utterances are
carefully pronounced when recording the isolated audio files for training, the
diction within the improvisation files is sometimes erratic, often mixing both
sounds in a hybrid utterance that the algorithm misclassifies in the testing
stage. A greater pre-emphasis in the higher frequency bands for feature
extraction or classification could improve accuracies for many participants.

##  5 Conclusions

We have explored here several approaches for effective and efficient Query By
Vocal Percussion. Different machine learning models and data augmentation
routines have been tested, and several combinations of their hyperparameters
have been investigated and later presented in the results section. We have
compared the approaches in terms of effectiveness, efficiency, stability and
interpretability. We found that the combination of data augmentation and
convolutional neural networks achieved the best classification accuracies in
our setting, at the expense of significantly longer training time and lower
performance stability to random initialisation. In the future, we intend to
carry out experiments in high-resolution percussive onset detection and build
a quick and accurate user-based QVP system with our results here in mind.

Acknowledgments

$\begin{array}[]{l}\includegraphics[scale={0.1}]{EU_Low.jpg}\end{array}$ This
project has received funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Sk$\l{}$odowska-Curie grant agreement
No. 765068

## References

  * [1] H. Marko and T. Virtanen, “Audio Query by Example Using Similarity Measures between Probability Density Functions of Features,” _EURASIP Journal on Audio, Speech, and Music Processing_ , vol. 2010, 01 2010.
  * [2] A. Mehrabi, K. Choi, S. Dixon, and M. Sandler, “Similarity measures for vocal-based drum sample retrieval using deep convolutional auto-encoders,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 356–360.
  * [3] A. Ramires, R. Penha, and M. E. P. Davies, “User Specific Adaptation in Automatic Transcription of Vocalised Percussion,” _arXiv:1811.02406_ , 2018.
  * [4] A. Hazan, “Towards automatic transcription of expressive oral percussive performances,” in _Proceedings of the 10th international conference on Intelligent user interfaces - IUI ’05_. San Diego, California, USA: ACM Press, 2005.
  * [5] D. Stowell and M. D. Plumbley, “Delayed decision-making in real-time beatbox percussion classification,” _Journal of New Music Research_ , 2010.
  * [6] E. Sinyor, C. M. Rebecca, D. Mcennis, and I. Fujinaga, “Beatbox classification using ACE,” in _Proceedings of the International Conference on Music Information Retrieval_ , 2005.
  * [7] T. Nakano, J. Ogata, M. Goto, and Y. Hiraga, “A Drum Pattern Retrieval Method by Voice Percussion,” in _5th International Conference on Music Information Retrieval (ISMIR), Barcelona, Spain_ , 2004.
  * [8] A. Kapur, M. Benning, and G. Tzanetakis, “Query-By-Beat-Boxing: Music Retrieval for the DJ,” in _5th International Conference on Music Information Retrieval (ISMIR), Barcelona, Spain_ , 2004.
  * [9] A. F. S. Ramires, “Automatic Transcription of Drums and Vocalised percussion,” _Universidade do Porto_ , 2017.
  * [10] A. Delgado, S. McDonald, N. Xu, and M. Sandler, “A New Dataset for Amateur Vocal Percussion Analysis,” in _Proceedings of the 14th International Audio Mostly Conference: A Journey in Sound_ , 2019, pp. 17–23.
  * [11] B. McFee, C. Raffel, D. Liang, D. P. Ellis, M. McVicar, E. Battenberg, and O. Nieto, “librosa: Audio and music signal analysis in python,” in _Proceedings of the 14th python in science conference_ , vol. 8, 2015.
  * [12] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg _et al._ , “Scikit-learn: Machine learning in Python,” _Journal of machine learning research_ , vol. 12, no. Oct, pp. 2825–2830, 2011.
  * [13] D. S. Park, W. Chan, Y. Zhang, C.-C. Chiu, B. Zoph, E. D. Cubuk, and Q. V. Le, “Specaugment: A simple data augmentation method for automatic speech recognition,” _arXiv preprint arXiv:1904.08779_ , 2019.
  * [14] F. Murtagh, “Multilayer perceptrons for classification and regression,” _Neurocomputing_ , vol. 2, no. 5-6, pp. 183–197, 1991.
  * [15] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning internal representations by error propagation,” California Univ San Diego La Jolla Inst for Cognitive Science, Tech. Rep., 1985.
  * [16] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [17] E. Hoffer and N. Ailon, “Deep metric learning using triplet network,” in _International Workshop on Similarity-Based Pattern Recognition_. Springer, 2015, pp. 84–92.
  * [18] A. Hermans, L. Beyer, and B. Leibe, “In defense of the triplet loss for person re-identification,” _arXiv preprint arXiv:1703.07737_ , 2017.