# Fast Decentralized Linear Functions Over
Edge Fluctuating Graphs

Siavash Mollaebrahim, , and Baltasar Beferull-Lozano S. Mollaebrahim and B.
Lozano are with the Intelligent Signal Processing and Wireless Networks
(WISENET) Center, University of Agder, Norway. e-mail:{siavash.mollaebrahim,
<EMAIL_ADDRESS>
###### Abstract

Implementing linear transformations is a key task in the decentralized signal
processing framework, which performs learning tasks on data sets distributed
over multi-node networks. That kind of network can be represented by a graph.
Recently, some decentralized methods have been proposed to compute linear
transformations by leveraging the notion of graph shift operator, which
captures the local structure of the graph. However, existing approaches have
some drawbacks such as considering some special instances of linear
transformations, or reducing the family of transformations by assuming that a
shift matrix is given such that a subset of its eigenvectors spans the
subspace of interest. In contrast, this paper develops a framework for
computing a wide class of linear transformations in a decentralized fashion by
relying on the notion of graph shift operator. The main goal of the proposed
method is to compute the desired linear transformation in a small number of
iterations. To this end, a set of successive graph shift operators is
employed, then, a new optimization problem is proposed whose goal is to
compute the desired transformation as fast as possible. In addition, usually,
the topology of the networks, especially the wireless sensor networks, change
randomly because of node failures or random links. In this paper, the effect
of edge fluctuations on the performance of the proposed method is studied. To
deal with the negative effect of edge fluctuations, an online kernel-based
method is proposed which enables nodes to estimate the missed values with
their at hand information. The proposed method can also be employed to
sparsify the network graph or reduce the number of local exchanges between
nodes, which saves sensors power in the wireless sensor networks.

###### Index Terms:

graph signal processing, decentralized linear transformations, graph filters.

##  I Introduction

Processing the data sets defined over non-regular domains is the main
application of the graph signal processing (GSP) field. Those data sets are
characterized by graphs, and traditional signal processing notions such as
frequency analysis and sampling have been extended to process signals defined
over graphs (graph signals), i.e. data associated with the nodes of the
network [1, 2, 3]. Graph signals exist in different fields such as social
networks and wireless sensor networks (WSN).

Recently, a number of decentralized methods have been proposed to compute some
special instances of linear transformations, as a common application of
decentralized signal processing. Projecting signals onto a low-dimensional
subspace [4, 5] or interpolating bandlimited signals [6] are examples of those
linear transformations. For instance, the authors in [7] propose a
decentralized approach to implement the subspace projection task, which
encompasses a number of learning tasks in WSN, such as estimation and
denoising [8]. In the aforementioned approach [7], each node at every
iteration linearly combines its iterate with the ones of its neighbors. Then,
by optimizing a criterion that quantifies asymptotic convergence, the
coefficients of those linear combinations are obtained. However, this
algorithm can only accommodate a limited set of topologies [9], its
convergence is asymptotic, and in general requires a relatively large number
of local information exchanges. Those issues are addressed in [10, 11, 4] by
proposing some decentralized GSP-based approaches for the average consensus
task, which is a special case of subspace projection.

Moreover, some decentralized GSP-based approaches employ Graph filters (GFs)
to compute the linear transformations. GFs are expressed by polynomials of the
so-called graph shift matrices, which capture the local structure of the graph
[2]. Applying GFs gives rise to decentralized implementations, which means
that the nodes of the network can implement the graph filter by exchanging
only local information with their neighbors. Graph filters have been used for
different applications such as graph signal analysis [12, 2] reconstruction
[13, 14, 6, 15] and denoising [16, 17, 5, 18].

In [17], GFs are designed to compute pre-specified linear transformations.
Nevertheless, the GFs design in [17] for rank-1 projections or for cases,
where projections share their eigenvectors with given graph shift matrices
such as the Laplacian or adjacency matrices, which limits the possible linear
functions that can be implemented. Moreover, the aforementioned method
considers the design of GFs only for symmetric topologies. On the other hand,
due to interference and background noise, WSNs are often characterized by
asymmetric packet losses, leading to network topologies that can be viewed as
asymmetric graphs [19, 20]. To approximate the linear transformations, a
decentralized approach has been proposed based on the so-called edge-variant
(EV) GFs in [21]. This method uses a different graph shift matrix in each
iteration of the graph filtering which enables nodes to weigh the signal of
their neighbors with different values. Furthermore, in [5] a decentralized
subspace projection method via GFs for symmetric topologies has been proposed
by optimizing a criterion that yields convergence to the subspace projection
in a nearly minimal number of iterations.

In this paper, we propose an approach to compute the linear transformations as
fast as possible i.e. in a small number of local exchanges. Our approach is
based on a sequence of different shift operators are employed to compute the
desired transformation matrix. In our method, the difference between the
sequence of shift operators and the transformation matrix at each round is
minimized. In order to decrease the required number of rounds which leads to
fast convergence, a weighted-sum method is proposed, where larger weights are
assigned to the last rounds to enforce the error of those rounds to be
smaller.

In addition, edge fluctuation is a critical issue in the decentralized GSP-
based approaches, and it negatively affects the performance of them. For
instance, in WSN, a number of network edges are missed when edge fluctuations
occur due to node failures. In this case, some nodes miss a portion of their
neighbors’ information, which creates a deviation in the nodes’ exchange
information procedure. The effect of graph perturbations on _finite impulse
response_ (FIR) graph filters has been investigated in [22, 23, 24, 25]. Those
works analyze the effect of graph perturbation based on the random edge
sampling (RES) graph model [23], where the activation probability of graph
edges is assumed to be known.

In summary, the contributions that we provide, can be summarized as follows:

  * •

In contrast with [26, 5], where only symmetric topologies have been
considered, the proposed method designs shift operators to compute the linear
transformations over different topologies, including both symmetric and
asymmetric graphs. Compared to [21], the proposed method needs a noticeably
smaller number of local exchanges, leading to a faster convergence and
computes the desired linear transformation exactly in scenarios where the
method in [21] approximates them.

  * •

Furthermore, this work considers the issue of edge fluctuations on the problem
of computing the linear transformations, which is not considered by existing
works [26, 5, 21]. In this paper, the effect of edge fluctuation on the
proposed approach is studied. Then, to deal with the negative effect of edge
fluctuations, an efficient online method is proposed, where each node
estimates the missing values (due to the fluctuation of the edges) with its
at-hand information, by employing an online kernel-based learning approach. In
the proposed method, each node estimates the missing values just by applying
some simple local calculations.

  * •

Finally, the proposed estimator enables us to randomly remove a number of
edges at each iteration to sparsify the graph, which gives rise to saving
power of sensor nodes in WSN, increasing their lifetime, since they are
usually empowered by batteries. By using our method, the sensor nodes spend
less power to relay their information to their neighbours by decreasing the
number of graph edges in each iteration. Indeed, for each node in each
iteration, we can randomly remove some edges, and then by applying the
proposed estimator, nodes can estimate their neighbours’ information
corresponded to the removed edges online. This advantage can be further
enhanced if the proposed estimators are trained during a certain number of
iterations, so that the local exchanges between nodes stop, and instead the
nodes update their information based on an estimate of their neighbours’
information. By this new approach, nodes power is saved since the number of
local exchanges is decreased.

Notation: Bold uppercase (lowercase) letters denote matrices (column vectors),
while $(.)^{\top}$ stands for matrix transposition. $({a})_{n}$ and
$(\mathbf{A})_{n,n^{\prime}}$ denote the $n$-th entry of $\mathbf{a}$ and the
entries of $\mathbf{A}$, respectively. Also, $\mathbf{1}$ is the all ones
vector and $\mathbf{e}_{n}$ is the basis vector whose all entries are zero
except the $n$-th one which equals one. The indicator function and the
Hadamard product are denoted by $\mathbbm{1}$ and $\circ$, respectively.
Moreover, $\Pi$ and $\mathrm{vec}$ represent the product of a collection of
terms and the vectorization operation. Finally, $\text{trace}(\mathbf{A})$ and
the Frobenius norm are denoted by $\mathrm{tr}(\mathbf{A})$ and $||.||_{F}$,
respectively.

##  II Preliminaries

###  II-A System Model

Consider $N$ networked sensor nodes which are able to exchange messages with
their neighbour nodes to which they are connected. The network is modeled as a
directed connected graph $\mathcal{G}(\mathcal{V},\mathcal{E})$, where the
vertex (node) set $\mathcal{V}:={1,2,\cdots,N}$ corresponds to the network
nodes, and $\mathcal{E}$ represents the set of edges. The $n^{\prime}$-th
vertex $v_{n^{\prime}}$ is connected to $v_{n}$ if there is a directed edge
from $v_{n^{\prime}}$ to $v_{n}$ $(v_{n^{\prime}},{v}_{n})\in\mathcal{E}$, but
this does not imply that $v_{n}$ is connected to $v_{n^{\prime}}$ unless
$(v_{n},{v}_{n^{\prime}})\in\mathcal{E}$. The in-neighborhood of the $n$-th
node is defined as the set of nodes connected to it, which can be denoted as
$\mathcal{N}_{n}=({v_{n^{\prime}}|(v_{n^{\prime}},{v}_{n})\in\mathcal{E}}$).
In the GSP context, vector $\mathbf{x}=[x_{1},x_{2},\cdots,x_{N}]^{\top}$ is
refereed to as a _graph signal_ where the $i$-th component of $\mathbf{x}$
represents the signal value at the $i$-th node of $\mathcal{G}$.

###  II-B Decentralized linear transformations via graph filters

In this section, we briefly review implementing linear transformations via
GFs. _Finite impulse response_ (FIR) graph filtering [2] involves two steps,
in the first step, a sequence of graph signals
$\mathbf{x}^{(l)}=[x_{1}^{(l)},x_{2}^{(l)},\cdots,x_{N}^{(l)}]^{\top},l=0,\cdots,L$
is obtained via $L$ number of local exchanges between nodes where
$\mathbf{x}^{(0)}=\mathbf{x}$. Particularly, the graph signal at the $l$-th
iteration i.e. $\mathbf{x}^{(l)}$ is computed based on the graph signal of the
previous iteration i.e. $\mathbf{x}^{(l-1)}$ such that
$x_{n}^{(l)}=\sum_{n^{\prime}\in\mathcal{N}_{n}}s_{n,n^{\prime}}x_{n^{\prime}}^{(l-1)}$
where $s_{n,n^{\prime}}$ is the coefficient used for aggregation of
information between the $n,n^{\prime}$ nodes. Thus, in matrix form, we have
$\mathbf{x}^{(l)}=\mathbf{S}\mathbf{x}^{(l-1)}$ where
$(\mathbf{S})_{n,n^{\prime}}=0$ if
$({v}_{n},{v}_{n^{\prime}})\notin\mathcal{E}$.

In the GSP context, $\mathbf{S}$ is called the shift operator [2, 27]. The
graph shift operator is a decentralized operator because each node updates its
information only via information of its neighbours. Examples of the graph
shift operator are the adjacency and Laplacian matrices of $\mathcal{G}$ [17].
In the second step, outputs of all iterations computed in the first step are
linearly aggregated as follows

$\displaystyle\mathbf{y}=\mathbf{H}\mathbf{x}\overset{\Delta}{=}\sum_{l=0}^{L-1}c_{l}\mathbf{S}^{l}\mathbf{x}$
(1)

where $\\{c_{l}\\}_{l=0}^{L-1}$ are the filter coefficients, and $L$ is the
order of the filter [2]. The graph filter in (1) has been used in [17] to
implement a pre-specified transformation matrix $\mathbf{T}$. Alternatively,
in [17] the graph filter (1) is designed only for symmetric topologies such
that $\sum_{l=0}^{L-1}c_{l}\mathbf{S}^{l}=\mathbf{T}$ for rank-1 projections
or for scenarios where $\mathbf{T}$ and the given shift matrix such as the
Laplacian or the Adjacency matrix are simultaneously diagonaizable i.e. they
share the same eigenvectors.

The concept of FIR graph filters is generalized in [21] by introducing the
edge-variant (EV) graph filters where a different set of weights is used in
each iteration. The edge weighting matrices share the same support with the
adjacency matrix of the network graph. In [21], the edge-variant graph filter
is used to approximate linear transformations.

##  III Proposed method: Successive shift operators

In this section, by using the notion of graph shift operator, we propose a new
approach to implement $\mathbf{T}\mathbf{x}$ in a decentralized manner as fast
as possible. To this end, a sequence of different shift operators is applied
so that
$\mathbf{S}_{l}\cdots\mathbf{S}_{2}\mathbf{S}_{1}\mathbf{x}=\mathbf{T}\mathbf{x}$,
where $\forall{\mathbf{S}_{i}},i=1,\cdots,l$, we have
$(\mathbf{S}_{i})_{nn^{\prime}}=0$ if
$(v_{n^{\prime}},v_{n})\not\in\mathcal{E}$.

As stated before, the proposed method leads to a decentralized approach. In
fact, after the first round of information exchange among nodes, the nodes
compute $\mathbf{x}^{(1)}:=\mathbf{S}_{1}\mathbf{x}$. Then, at the second
round,
$\mathbf{x}^{(2)}:=\mathbf{S}_{2}\mathbf{x}^{(1)}=\mathbf{S}_{2}\mathbf{S}_{1}\mathbf{x}$.
The procedure is repeated for $l$ iterations. By properly designing
$\\{\mathbf{S}_{i}\\}_{i=1}^{l}$, $\mathbf{T}\mathbf{z}$ is exactly computed
or $\Pi_{i=1}^{l}\mathbf{S}_{i}\mathbf{x}\approx\mathbf{T}\mathbf{x}$ is
obtained with a small error.

To compute $\mathbf{T}$ as efficient as possible i.e. with the smallest number
of iterations, we design the sequence of feasible graph shift matrices that
minimize the number of local exchanges, i.e. $l$ by solving the following
problem:

$\displaystyle\underset{\mathbf{S}_{1},\mathbf{S}_{2},\cdots,\mathbf{S}_{l}}{\text{min}}$
$\displaystyle l$ (2a) s. t.
$\displaystyle(\mathbf{S}_{i})_{n,n^{\prime}}=0\hskip 5.69054pt$
$\displaystyle\text{if}\hskip
5.69054pt(v_{n^{\prime}},v_{n})\not\in\mathcal{E},n,n^{\prime}=1,....,N,i=1,...,l$
(2b) $\displaystyle\mathbf{T}=\Pi_{i=1}^{l}\mathbf{S}_{i}$ (2c)

where constraint (2b) is added to the problem because the shift operators must
satisfy the topological constraints of the graph connectivity.

Optimization problem (2) is intractable. In order to address this issue, an
upper limit for $l$ is considered such that $l\in\\{1,\cdots,L\\}$. Then, to
design the shift operators that compute $\mathbf{T}$ via
$\Pi_{l}\mathbf{S}_{l}$ as fast as possible, the error at each round can be
minimized i.e. the error between $\mathbf{T}$ and $\Pi_{l}\mathbf{S}_{l}$,
which can be attained by solving the following optimization problem:

$\displaystyle\underset{\mathbf{S}_{1},\mathbf{S}_{2},\cdots,\mathbf{S}_{L}}{\text{min}}$
$\displaystyle\sum_{l=1}^{L}{{\left\|\mathbf{T}-\Pi_{i=1}^{l}\mathbf{S}_{i}\right\|}_{F}^{2}}$
(3a) s. t. $\displaystyle(\mathbf{S}_{i})_{n,n^{\prime}}=0\hskip 5.69054pt$
$\displaystyle\text{if}\hskip
5.69054pt(v_{n^{\prime}},v_{n})\not\in\mathcal{E},n,n^{\prime}=1,....,N,i=1,...,L$
(3b)

In optimization problem (3), we put the same emphasis on all rounds, which
means that all of them are treated equally. However, if we put more emphasise
on the later rounds, we expect that a faster approach is obtained. By this
approach, we accept larger error on the earlier rounds to obtain smaller error
at the later ones, which can be implemented by assigning different weights to
the terms in cost function (3a). To end this, larger weights are assigned to
the later rounds. Thus, we have:

$\displaystyle\underset{\mathbf{S}_{1},\mathbf{S}_{2},\cdots,\mathbf{S}_{L}}{\text{min.}}$
$\displaystyle\sum_{l=1}^{L}\omega_{l}{{\left\|\mathbf{T}-\Pi_{i=1}^{l}\mathbf{S}_{i}\right\|}_{F}^{2}}$
(4a) s. t. $\displaystyle(\mathbf{S}_{i})_{n,n^{\prime}}=0\hskip 5.69054pt$
$\displaystyle\text{if}\hskip
5.69054pt(v_{n^{\prime}},v_{n})\not\in\mathcal{E},n,n^{\prime}=1,....,N,i=1,...,L$
(4b)

where $\boldsymbol{\omega}=[\omega_{1},\omega_{2},\cdots,\omega_{L}]^{\top}$
is the weight vector whose entries are non-negative increasing, and such that
$\sum_{l=1}^{L}\omega_{l}=1$.

Optimization problem (4) is non-convex with respect to all its variables;
however, the cost function (4a) is a strongly convex function with respect to
each $\mathbf{S}_{i}$ separately, i.e. when the other blocks of variables are
fixed. Therefore, we can apply the block coordinate descent (BCD) algorithm to
solve the optimization problem (4). The convergence of the BCD algorithm is
guaranteed for the optimization problem in (4) [28]. In this approach, all
blocks of variables are fixed except one of them, and the optimization problem
is solved with respect to the free block. This procedure continues until all
the matrices $\mathbf{S}_{1},\mathbf{S}_{2},\cdots,\mathbf{S}_{L}$ are
considered as the block variable of the optimization problem. The BCD
algorithm is repeated until a stopping criterion is satisfied.

Here, we try to obtain the closed form solution of each round of the BCD
algorithm, where one of shift matrices is the variable of optimization problem
(4) and the rest of shift matrices are fixed, the following optimization
problem is solved:

$\displaystyle{J^{*}_{j}}\overset{\Delta}{=}\underset{\mathbf{S}_{j}=\sum_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}S_{n,n^{\prime}}^{(j)}\mathbf{e}_{n}\mathbf{e}_{n^{\prime}}^{\top}}{\text{Inf}}$
$\displaystyle{J_{j}}(\mathbf{S}_{j})$ (5)

where
${J_{j}}(\mathbf{S}_{j})=\sum_{l=j}^{L}\omega_{l}{\left\|\mathbf{T}-\mathbf{S}_{l}\mathbf{S}_{l-1}\cdots\mathbf{S}_{j+1}\mathbf{S}_{j}\mathbf{S}_{j-1}\cdots\mathbf{S}_{1}\right\|}_{F}^{2}$,
$\mathbf{S}_{j}=\sum_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}S_{n,n^{\prime}}^{(j)}\mathbf{e}_{n}\mathbf{e}_{n^{\prime}}^{\top}$
satisfies constraint (4b), and $S_{n,n^{\prime}}^{(j)}$ is the entry that lies
in the $n$-th row and the $n^{\prime}$-th column of $\mathbf{S}_{j}$. Thus, we
have:

$\displaystyle{J^{*}_{j}}=\underset{\\{{S}^{(j)}_{n,n^{\prime}}\\}_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}}{\text{Inf}}$
$\displaystyle{J_{j}}(\sum_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}S_{n,n^{\prime}}^{(j)}\mathbf{e}_{n}\mathbf{e}_{n^{\prime}}^{\top})$
(6)

The term inside of (6) can be expressed in vector form by applying the
vectorization operator as follows:

$\displaystyle\mathrm{vec}(\sum_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}S_{n,n^{\prime}}^{(j)}\mathbf{e}_{n}\mathbf{e}_{n^{\prime}}^{\top})=\sum_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}S_{n,n^{\prime}}^{(j)}(\mathbf{e}_{n^{\prime}}\otimes\mathbf{e}_{n})$
(7)

where in (7), we use the fact that
$\mathrm{vec}(\mathbf{e}_{n}\mathbf{e}_{n^{\prime}}^{\top})=\mathbf{e}_{n^{\prime}}\otimes\mathbf{e}_{n}$.
Moreover, we have:

$\displaystyle\sum_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}S_{n,n^{\prime}}^{(j)}(\mathbf{e}_{n^{\prime}}\otimes\mathbf{e}_{n})=\underbrace{[\mathbf{e}_{n^{\prime}_{1}}\otimes\mathbf{e}_{n_{1}},\cdots,\mathbf{e}_{n^{\prime}_{E}}\otimes\mathbf{e}_{n_{E}}]}_{\mathbf{E}}{\mathbf{s}^{(j)}}$
(8)

where
$\mathbf{s}^{(j)}=\left[\begin{array}[]{ccc}S_{n_{1},n^{\prime}_{1}}^{(j)}&\cdots&S_{n_{E},n^{\prime}_{E}}^{(j)}\end{array}\right]^{\top}$.
Therefore, optimization problem (6) can be represented as follows:

$\displaystyle\underset{\\{S^{(j)}_{n,n^{\prime}}\\}_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}}{\text{inf}}J(\sum_{(v_{n^{\prime}},v_{n})\in\mathcal{E}}S_{n,n^{\prime}}^{(j)}\mathbf{e}_{n}\mathbf{e}_{n^{\prime}}^{\top})=$
$\displaystyle\underset{\mathbf{s}^{(j)}\in\mathbb{R}^{E}}{\text{inf}}J(\mathrm{vec}^{-1}(\mathbf{E}\mathbf{s}^{(j)}))$
(9)

Consequently, from (III), optimization problem (4) can be re-written as
follows:

$\displaystyle\underset{\mathbf{s}^{(j)}}{\text{inf}}$
$\displaystyle\sum_{l=j}^{L}\omega_{l}{{\left\|\mathbf{T}-\mathbf{S}_{l:j+1}\mathrm{vec}^{-1}(\mathbf{E}\mathbf{s}^{(j)}))\mathbf{S}_{j-1:1}\right\|}_{F}^{2}}$
(10)

where
$\mathbf{S}_{l:j}\overset{\Delta}{=}\mathbf{S}_{l}\mathbf{S}_{l-1}\cdots\mathbf{S}_{j}$.

Finally, the vectorized form of (10) equals:

$\displaystyle\underset{\mathbf{s}^{(j)}}{\text{inf}}$
$\displaystyle\sum_{l=j}^{L}\omega_{l}{{\left\|\mathrm{vec}(\mathbf{T})-(\mathbf{S}_{j-1:1}^{\top}\otimes\mathbf{S}_{l:j+1})\mathbf{E}\mathbf{s}^{(j)}\right\|}_{2}^{2}}$
(11)

The following algorithm represents the approach for solving (11). For the
stopping criterion
$|f(\mathbf{S}_{L:1}^{(i+1)})-f(\mathbf{S}_{L:1}^{(i)})|/|f(\mathbf{S}_{L:1}^{(i)})|<\varepsilon$
is used for a sufficiently small $\varepsilon$, where $f(\mathbf{S}_{L:1})$
equals cost function (3a), and $i$ denotes the $i$-th iteration of the BCD
algorithm.

input : $\mathcal{E},L,\omega_{1},\cdots,\omega_{L},\varepsilon$.

initialize : $\mathbf{S}_{2},\mathbf{S}_{3},\cdots,\mathbf{S}_{L}$.

1 while  _stopping criterion not satisfied_ do

2 for  _$j=1$ to $L$_ do

3 fix $\mathbf{S}_{m},m=1,\cdots,L,m\neq{j}$, obtain $\mathbf{s}^{(j)}$ via
(11)

4 end for

5

6 end while

output : $\mathbf{S}_{1},\mathbf{S}_{2},\cdots,\mathbf{S}_{L}$

Algorithm 1 Proposed solver

##  IV Edge fluctuation

###  IV-A Edge Fluctuation Model

In the previous section, we have assumed that all edges are fixed; however, in
practical scenarios, there are fluctuations on edges due to different reasons
such as random links or node failures in WSN, which leads to a random graph
topology. In this subsection, the effect of edge fluctuations on the proposed
method i.e. the successive method is studied. We consider that edge
fluctuations occur randomly and independently across edges. Please note that
in this section, to provide an analysis of the effect of edge fluctuations, we
assume that the estimation of link activation probabilities is available.
However, knowing the link activation probabilities is not required for our
proposed method to deal with the edge fluctuations effect. Moreover, we assume
that the spectral norm of all graph shift operators is upper-bounded i.e.
$||\mathbf{S}_{i}||_{2}\leq\rho<\infty\quad\forall i=1,\cdots{L}$. This
assumption is practical because it implies that the graph of interest has
finite edge weights [24].

When an edge is missed in a certain of iteration of the successive method, the
corresponded destination node to the missing edge loses the information of its
neighbor, which generates a deviation. The shift operator corresponding to the
graph with perturbations can be modeled as
$\hat{\mathbf{S}}=\mathbf{S}+\tilde{\mathbf{S}}$ where $\tilde{\mathbf{S}}$ is
created because of edge fluctuations, and its non-zero entries equal the
corresponding entries of $\mathbf{S}$ with a negative sign. Thus, the signal
output for the first iteration of the successive method can be expressed

$\displaystyle\mathbf{y}^{(1)}=(\mathbf{S}_{1}+\tilde{\mathbf{S}}_{1})\mathbf{x}=\mathbf{S}_{1}\mathbf{x}+\mathbf{z}_{1}$
(12)

where $\mathbf{z}_{1}\overset{\Delta}{=}\tilde{\mathbf{S}}_{1}\mathbf{x}$. We
assume that $\mathbf{z}_{1}$ is a realization of a random process, with mean
$\mathbf{m}_{z_{1}}$ and covariance matrix $\mathbf{\Sigma}_{z_{1}}$. Note
that $\tilde{\mathbf{S}}_{1}$ can be modeled by using the activation link
probabilities i.e.
$\mathbb{E}[\tilde{\mathbf{S}}_{1}]=(\mathbf{1}-\mathbf{P}_{\text{ac}})\circ(-\mathbf{S}_{1})$,
where $\mathbf{P}_{\text{ac}}$ is the link activation probability matrix.
Thus, since $\mathbf{x}$ is fixed we have
$\mathbf{m}_{z_{1}}=\mathbf{m}_{\tilde{\mathbf{S}}_{1}}\mathbf{x}$, where
$\mathbf{m}_{\tilde{\mathbf{S}}_{1}}=(\mathbf{P}_{\text{ac}}-\mathbf{1})\circ(\mathbf{S}_{1})$
is the mean of ${\tilde{\mathbf{S}}_{1}}$. In addition, for the covariance of
$\mathbf{z}_{1}$, we have

$\displaystyle\mathbf{\Sigma}_{z_{1}}=\mathbb{E}[\mathbf{z}_{1}\mathbf{z}_{1}^{\top}]-\mathbb{E}[\mathbf{z}_{1}]\mathbb{E}[\mathbf{z}_{1}^{\top}]$
$\displaystyle\quad\quad=\mathbb{E}[\tilde{\mathbf{S}}_{1}\mathbf{x}\mathbf{x}^{\top}\tilde{\mathbf{S}}_{1}^{\top}]-\mathbf{m}_{\mathbf{z}_{1}}\mathbf{m}_{\mathbf{z}_{1}}^{\top}$
(13)

The $ij$-th element of
$\mathbb{E}[\tilde{\mathbf{S}}_{1}\mathbf{x}\mathbf{x}^{\top}\tilde{\mathbf{S}}_{1}^{\top}]$
equals
$\mathbb{E}(\tilde{\mathbf{s}}_{1_{i}}^{\top}\mathbf{X}\tilde{\mathbf{s}}_{1_{j}})$,
where $\mathbf{X}=\mathbf{x}\mathbf{x}^{\top}$ and
$\tilde{\mathbf{S}}_{1}=[\tilde{\mathbf{s}}_{1_{1}}|\tilde{\mathbf{s}}_{1_{2}}|\cdots|\tilde{\mathbf{s}}_{1_{N}}]^{\top}$.
From the cyclic property of trace, and the fact that trace is a linear
operator, we have

$\displaystyle(\mathbb{E}[\tilde{\mathbf{S}}_{1}\mathbf{x}\mathbf{x}^{\top}\tilde{\mathbf{S}}_{1}^{\top}])_{i,j}=\mathbb{E}[\tilde{\mathbf{s}}_{1_{i}}^{\top}\mathbf{X}\tilde{\mathbf{s}}_{1_{j}}]=\mathrm{tr}(\mathbf{X}\mathbb{E}(\tilde{\mathbf{s}}_{1_{j}}\tilde{\mathbf{s}}_{1_{i}}^{\top}))$
$\displaystyle=\mathrm{tr}(\mathbf{X}\mathbf{C}_{1_{{ji}}})$ (14)

where
$\mathbf{C}_{1_{{ji}}}=\mathbb{E}[\mathbf{s}_{1_{j}}\mathbf{s}_{1_{i}}^{\top}]$
is the cross-correlation between
$\tilde{\mathbf{s}}_{1_{j}},\tilde{\mathbf{s}}_{1_{i}}$. Consequently:

$\displaystyle(\mathbf{\Sigma}_{z_{1}})_{i,j}=\mathrm{tr}(\mathbf{X}\mathbf{C}_{1_{{ji}}})-(\mathbf{M}_{\mathbf{z}_{1}})_{i,j}\quad\forall{i,j=1,\cdots,N}$
(15)

where
$\mathbf{M}_{\mathbf{z}_{1}}=\mathbf{m}_{\mathbf{z}_{1}}\mathbf{m}_{\mathbf{z}_{1}}^{\top}$.
In the next iteration, we have:

$\mathbf{y}^{(2)}=(\mathbf{S}_{2}+\tilde{\mathbf{S}}_{2})\mathbf{y}^{(1)}=\mathbf{S}_{2}\mathbf{S}_{1}\mathbf{x}+\mathbf{S}_{2}\mathbf{z}_{1}+\mathbf{z}_{2}$
(16)

where
$\mathbf{z}_{2}\overset{\Delta}{=}\tilde{\mathbf{S}}_{2}\mathbf{S}_{1}\mathbf{x}+\tilde{\mathbf{S}}_{2}\mathbf{z}_{1}$,
which implies that $\mathbf{z}_{2}$ is correlated with $\mathbf{z}_{1}$.
Indeed,
$\mathbf{m}_{\mathbf{z}_{2}}=\mathbf{m}_{\tilde{\mathbf{S}}_{2}}(\mathbf{S}_{1}+\mathbf{m}_{\tilde{\mathbf{S}}_{1}})\mathbf{x}$
and

$\displaystyle\mathbf{\Sigma}_{\mathbf{z}_{2}}=\mathbb{E}[\mathbf{z}_{2}\mathbf{z}_{2}^{\top}]=\mathbb{E}[(\tilde{\mathbf{S}}_{2}\mathbf{S}_{1}\mathbf{x}+\tilde{\mathbf{S}}_{2}\mathbf{z}_{1})$
$\displaystyle\times(\mathbf{x}^{\top}\mathbf{S}_{1}^{\top}\tilde{\mathbf{S}}_{2}^{\top}+\mathbf{z}_{1}^{\top}\tilde{\mathbf{S}}_{2}^{\top})]-\mathbf{m}_{\mathbf{z}_{2}}\mathbf{m}_{\mathbf{z}_{2}}^{\top}$
(17)

From $\mathbf{z}_{1}\overset{\Delta}{=}\tilde{\mathbf{S}}_{1}\mathbf{x}$ and
(IV-A), we have that

$\displaystyle\mathbf{\Sigma}_{\mathbf{z}_{2}}=\mathbb{E}[\tilde{\mathbf{S}}_{2}\mathbf{S}_{1}\mathbf{x}\mathbf{x}^{\top}\mathbf{S}_{1}^{\top}\tilde{\mathbf{S}}_{2}^{\top}+\tilde{\mathbf{S}}_{2}\mathbf{S}_{1}\mathbf{x}\mathbf{x}^{\top}\tilde{\mathbf{S}}_{1}^{\top}\tilde{\mathbf{S}}_{2}^{\top}+$
$\displaystyle\tilde{\mathbf{S}}_{2}\tilde{\mathbf{S}}_{1}\mathbf{x}\mathbf{x}^{\top}\tilde{\mathbf{S}_{1}}^{\top}\tilde{\mathbf{S}}_{2}^{\top}+\tilde{\mathbf{S}}_{2}\tilde{\mathbf{S}}_{1}\mathbf{x}\mathbf{x}^{\top}\mathbf{S}_{1}^{\top}\tilde{\mathbf{S}}_{2}^{\top}]-\mathbf{M}_{\mathbf{z}_{2}}$
(18)

where
$\mathbf{M}_{\mathbf{z}_{2}}=\mathbf{m}_{\mathbf{z}_{2}}\mathbf{m}_{\mathbf{z}_{2}}^{\top}$.
By taking expectations and the linearity property of the expectation operator,
we have:

$\displaystyle\mathbf{\Sigma}_{\mathbf{z}_{2}}=\mathbb{E}[\tilde{\mathbf{S}}_{2}\mathbf{S}_{1}\mathbf{X}\mathbf{S}_{1}^{\top}\tilde{\mathbf{S}}_{2}^{\top}]+\mathbb{E}[\tilde{\mathbf{S}}_{2}\mathbf{S}_{1}\mathbf{X}\mathbf{m}_{\tilde{\mathbf{S}}_{1}}^{\top}\tilde{\mathbf{S}}_{2}^{\top}]+$
$\displaystyle\mathbb{E}[\tilde{\mathbf{S}}_{2}\tilde{\mathbf{S}}_{1}\mathbf{X}\tilde{\mathbf{S}_{1}}^{\top}\tilde{\mathbf{S}}_{2}^{\top}]+\mathbb{E}[\tilde{\mathbf{S}}_{2}\mathbf{m}_{\tilde{\mathbf{S}}_{1}}\mathbf{X}\mathbf{S}_{1}^{\top}\tilde{\mathbf{S}}_{2}^{\top}]-\mathbf{M}_{\mathbf{z}_{2}}$
(19)

Similar to (IV-A) and (15), for each entry of
$\mathbf{\Sigma}_{\mathbf{z}_{2}}$, we have

$\displaystyle(\mathbf{\Sigma}_{\mathbf{z}_{2}})_{i,j}=\mathrm{tr}((\mathbf{S}_{1}\mathbf{X}\mathbf{S}_{1}^{\top}+\mathbf{S}_{1}\mathbf{X}\mathbf{m}_{\tilde{\mathbf{S}}_{1}}^{\top}+\boldsymbol{\chi}_{1}+\mathbf{m}_{\tilde{\mathbf{S}}_{1}}\mathbf{X}\mathbf{S}_{1}^{\top})\mathbf{C}_{2_{{ji}}})$
$\displaystyle-(\mathbf{M}_{\mathbf{z}_{2}})_{i,j}\quad\forall{i,j=1,\cdots,N}$
(20)

where
$\boldsymbol{\chi}_{1}=\mathbb{E}[\tilde{\mathbf{S}}_{1}\mathbf{X}\tilde{\mathbf{S}}_{1}^{\top}]$
and $\mathbf{C}_{2_{{ji}}}$ is the cross-correlation between
$\tilde{\mathbf{s}}_{2_{i}},\tilde{\mathbf{s}}_{2_{j}}$. By using a similar
approach to the one used in (IV-A), we have that
$(\boldsymbol{\chi}_{1})_{i,j}=\mathrm{tr}(\mathbf{X}\mathbf{C}_{1_{{ji}}}),\forall{i,j=1,\cdots,N}$.

After $L$ iterations, the output of the successive method under the edge
fluctuation model equals:

$\mathbf{y}=\mathbf{H}\mathbf{x}+\sum_{i=1}^{L-1}(\Pi_{j=i}^{L-1}\mathbf{S}_{j+1})\mathbf{z}_{i}+\mathbf{z}_{L}$
(21)

where $\mathbf{H}=\Pi_{i=1}^{L}\mathbf{S}_{i}$. The first and second order
statistics of $\mathbf{z}_{i},\forall{i=3,\cdots,L}$ can be obtained by the
similar approach as the one used for $\mathbf{z}_{1},\mathbf{z}_{2}$.

Next, we show that, under the edge fluctuation model, the mean square error
(MSE) of the the successive approach output after $L$ iterations expressed as
follows is upper bounded.

$\displaystyle\mathbb{E}[||\boldsymbol{\Omega}||_{2}^{2}]=\mathbb{E}[||\mathbf{z}_{L}+\sum_{i=1}^{L-1}(\Pi_{j=i}^{L}\mathbf{S}_{j+1})\mathbf{z}_{i}||_{2}^{2}]$
(22)

where $\boldsymbol{\Omega}=\mathbf{y}-\mathbf{H}\mathbf{x}$.

Proposition 1: The MSE of the output of the successive approach in (21) under
the edge fluctuation model is upper-bounded by:

$\displaystyle\mathbb{E}[]||\boldsymbol{\Omega}||_{2}^{2}]\leq\mathbb{E}[\mathrm{tr}(\boldsymbol{\Psi})]+2\sum_{j=2}^{L}\rho^{L-j+1}\mathrm{tr}(\boldsymbol{\Sigma}_{z_{j-1}}+\mathbf{m}_{z_{j-1}})+$
$\displaystyle\mathrm{tr}(\boldsymbol{\Sigma}_{z_{L}}+\mathbf{m}_{z_{L}})$
(23)

where
$\boldsymbol{\Psi}=\sum_{i=1}^{L-1}(\Pi_{j=i}^{L-1}\mathbf{S}_{j+1})\mathbf{z}_{i}\sum_{k=1,k\neq{i}}^{L-1}\mathbf{z}_{k}^{\top}(\Pi_{m=k}^{L-1}\mathbf{S}_{m+1})^{\top}\\\
+\mathbf{z}_{L}\sum_{k=1}^{L-1}\mathbf{z}_{k}^{\top}(\Pi_{j=k}^{L-1}\mathbf{S}_{j+1})^{\top}+\sum_{i=1}^{L-1}(\Pi_{j=i}^{L-1}\mathbf{S}_{j+1})\mathbf{z}_{i}\mathbf{z}_{L}^{\top}$.

###### Proof.

Please see Appendix A. ∎

From Proposition 1, the MSE of the successive approach output under the edge
fluctuation model is upper-bounded by $\rho$ and some other parameters i.e.
$\boldsymbol{\Psi},\mathbf{m}_{z_{i}}$ and
$\mathbf{\Sigma}_{z_{i}},\forall{i=1,\cdots,N}$, which are dependent on the
statistics of ${\tilde{\mathbf{S}}_{i}},\forall{i=1,\cdots,N}$ and the shift
matrices based on Proposition 1.

###  IV-B Proposed method: the missing values estimator

In this section, to deal with the edge fluctuations effect, we focus to
estimate the missing values via a kernel-based estimator. For this, the main
task can be expressed as follows: given $\mathbf{x},\mathcal{G},\mathbf{H}$,
the goal is to estimate the missing values due to the edge fluctuations. To
this end, an online kernel-based method is proposed such that the node whose
neighbors’ information is missed, estimates the missing value based on its own
available information.

Let us consider that in the $i$-th iteration of the successive approach, the
$n$-th node misses the information from the $n^{\prime}$-th node which is
given by

$\displaystyle
x_{n^{\prime}}^{(i-1)}=\sum_{n^{\prime\prime}\in\mathcal{N}_{n^{\prime}}}({\mathbf{S}_{i})_{n^{\prime},n^{\prime\prime}}}x_{n^{\prime\prime}}^{(i-2)}+({\mathbf{S}_{i})_{n^{\prime},n^{\prime\prime}}}x_{n^{\prime}}^{(i-2)}$
(24)

From (24), we can conclude that $x_{n^{\prime}}^{(i-2)}$ and $x_{n}^{(i-2)}$
are the most relevant information that is available for the $n$-th node to
estimate the missing value. Rest information of the $n$-th node could be also
informative since the $n$ and $n^{\prime}$ nodes can have one-hop or two-hops
common neighbors. Moreover, due to the fact that the network is connected, all
nodes are connected to each other with a number of hops. Consequently, in our
method, the $n$-th node utilizes
$\mathbf{u}_{n}=[x_{n^{\prime}}^{(i-2)},x_{n}^{(i-2)},\mathbf{c}_{n}^{(i-2)}]^{\top}\in\mathbb{R}^{|\mathcal{N}_{n}|+1\times{1}}$
to estimate the missing value (information of the $n^{\prime}$ node in the
$i-1$-th iteration) where $\mathbf{c}_{n}$ is a vector that contains the
neighbours of the $n$-th node except the $n^{\prime}$-th one.

To estimate the missing values, in supervised learning, an estimator function
is sought via the training samples. To end this, flexible non-parametric
models are needed for good results. A general class of models is based on
functions of the form [29]

$\displaystyle
f(\mathbf{q})=\sum_{k=1}^{K}\alpha_{k}\mathcal{F}(\mathbf{q},\mathbf{q}_{k})$
(25)

where $\mathcal{F}$ is a non-linear function,
$\alpha_{1},\cdots,\alpha_{K}\in\mathbb{R}$ are coefficients,
$\mathbf{q}_{1},\cdots,\mathbf{q}_{K}\in\mathbb{R}^{d}$ are the training
samples and $K$ is the number of samples. Kernel methods are a famous example
of an approach using functions of the form (25) [30].

In kernel-based learning, it is assumed that the estimator function belongs to
a reproducing kernel Hilbert space (RKHS) expressed as [30]

$\displaystyle\mathcal{H}:=\\{f:f(\mathbf{q})=\sum_{k=1}^{K}\alpha_{k}\kappa(\mathbf{q},\mathbf{q}_{k}),\alpha_{k}\in\mathbb{R}\\}$
(26)

where $\kappa:\mathbb{R}^{N}\times\mathbb{R}^{N}\to\mathbb{R}$ is a pre-
selected symmetric and positive definite function. Given $\mathbf{y}$ and
matrix samples $\mathbf{Q}=[\mathbf{q}_{1},\cdots.\mathbf{q}_{K}]$, the
kernel-based estimate is usually obtained by solving the following problem:

$\displaystyle\hat{f}=\underset{f\in\mathcal{H}}{\text{arg
min}}\\{1/K\sum_{k=1}^{K}\mathcal{L}(f(\mathbf{q}_{k})-y_{k})+\lambda||f||^{2}_{\mathcal{H}}\\}$
(27)

where $\mathcal{L}$ is a pre-selected cost function, for regression the least-
squares is utilized. Moreover,
$||f||^{2}_{\mathcal{H}}=\sum_{m=1}\sum_{m^{\prime}=1}\alpha_{m}\alpha_{m^{\prime}}\kappa(\mathbf{q}_{m}.\mathbf{q}_{m^{\prime}})$
denotes the RKHS norm used as a regularizer to deal with overfitting and
$\lambda$ is a regularization parameter.

Due to the fact that $\mathcal{H}$ is infinite dimensional, $\hat{f}$ cannot
be obtained directly by solving (27). However, based on the _representer
theorem_ [31], (27) admits a solution in form of

$\displaystyle\hat{f}(\mathbf{q})=\sum_{k=1}^{K}\alpha_{k}\kappa(\mathbf{q},\mathbf{q}_{k}):=\boldsymbol{\alpha}^{\top}\mathbf{t}(\mathbf{q})$
(28)

where
$\mathbf{t}(\mathbf{q}):=[\kappa(\mathbf{q},\mathbf{q}_{1}),\cdots,\kappa(\mathbf{q},\mathbf{q}_{K})]^{\top}$.
By substituting (28) into (27), $\boldsymbol{\alpha}$ is obtained via solving
the following problem:

$\displaystyle\underset{\boldsymbol{\alpha}\in\mathbb{R}^{K}}{\text{
min}}\\{1/K\sum_{k=1}^{K}\mathcal{L}(\boldsymbol{\alpha}^{\top}\mathbf{t}(\mathbf{q}_{k})-y_{k})+\lambda\boldsymbol{\alpha}^{\top}\boldsymbol{\Gamma}\boldsymbol{\alpha}\\}$
(29)

where $\boldsymbol{\Gamma}\in\mathbb{R}^{K\times{K}}$ whose
$(i,i^{\prime})$-th entry is $\kappa(\mathbf{q}_{i},\mathbf{q}_{i^{\prime}})$.

In our proposed method, each node for each of its neighbours employs a
function estimator to estimate the missing value. This can be obtained by
using estimators in form of (28) such that the $i$-th node employs the
following estimator for each of its neighbours:

$\displaystyle{f}_{ij}(\mathbf{u}_{i})=\sum_{m=1}^{K}\alpha_{m}^{(ij)}\kappa(\mathbf{u}_{i},\mathbf{u}_{m})$
(30)

where $j\in\mathcal{N}_{i}$, $\\{\mathbf{u}_{m}\\}_{m=1}^{K}$ is the training
samples. The estimators are trained by solving optimization problem (29) that
admits a closed-form solution as follows for kernel regression:

$\displaystyle\hat{\boldsymbol{\alpha}}=(\boldsymbol{\Gamma}+\lambda{K}\mathbf{I})^{-1}\mathbf{y}$
(31)

From (31), to obtain $\hat{\boldsymbol{\alpha}}$, each node for each of its
neighbors should invert a $K\times{K}$ matrix which is computationally heavy
especially when the number of samples is high. Indeed, each node, when an edge
fluctuation occurs, has to inverse a $K\times{K}$ matrix whose computation
complexity equals $\mathcal{O}(K^{3})$. Also, they have to save a high number
of scalars, for instance, the $i$-th node should save
${K}\times(|\mathcal{N}_{i}|+1)$ scalars.

To tackle this issue, we employ the random feature approach [32], which is an
approximation of $\kappa(\mathbf{u},\mathbf{u}_{m})$ obtained from the shift-
invariant kernels i.e.
$\kappa(\mathbf{u}_{n},\mathbf{u}_{m})=\kappa(\mathbf{u}_{n}-\mathbf{u}_{m})$.
In this paper, we employ shift-invariant and bounded kernels. The examples of
that kind of kernels are Gaussian, Laplacian and Cauchy kernels [32]. The
Fourier transforms of the shift-invariant kernels denoted by
$\Theta(\mathbf{w})$ exist in closed form as follows:

$\displaystyle\kappa(\mathbf{u}_{n},\mathbf{u}_{m})=\int\Theta(\mathbf{w})e^{j\mathbf{w}(\mathbf{u}_{n}-\mathbf{u}_{m})}d\mathbf{w}$
(32)

From the definition of the expected value, (32) equals
$\mathbb{E}_{\mathbf{w}}[e^{j\mathbf{w}(\mathbf{u}_{n}-\mathbf{u}_{m})}]$. The
main idea of the random feature approach is that by drawing a sufficient
number of samples $\\{\mathbf{w}_{i=1}^{D}\\}$ from $\Theta(\mathbf{w})$,
$\kappa(\mathbf{u}_{n},\mathbf{u}_{m})$ can be approximated by

$\displaystyle\hat{\kappa}(\mathbf{u}_{n},\mathbf{u}_{m})=\boldsymbol{\Delta}_{\mathbf{W}}^{\top}(\mathbf{u}_{n})\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u}_{m})$
(33)

where $\mathbf{W}:=[\mathbf{w}_{1},\cdots,\mathbf{w}_{D}]$ and

$\displaystyle\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u})=\frac{1}{\sqrt{D}}[\mathrm{sin}(\mathbf{w}_{1}^{\top}\mathbf{u}),\cdots,\mathrm{sin}(\mathbf{w}_{D}^{\top}\mathbf{u}),$
$\displaystyle\mathrm{cos}(\mathbf{w}_{1}^{\top}\mathbf{u}),\cdots,\mathrm{cos}(\mathbf{w}_{D}^{\top}\mathbf{u})]^{\top}$
(34)

From (32) and taking expectation from both sides of (33) with respect to
$\mathbf{W}$, we can conclude that $\hat{\kappa}(.)$ is unbiased. Moreover,
the variance of $\hat{\kappa}(.)$ is proportional with $1/D$, which implies
that it tends to zero when the number of random features tends to infinity
[33].

By substituting (33) to (30), we have the following estimator

$\displaystyle\hat{f}_{ij}(\mathbf{u}_{i})=\sum_{m=1}^{K}\alpha_{m}^{(ij)}\boldsymbol{\Delta}_{\mathbf{W}}^{\top}(\mathbf{u}_{m})\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u}_{i}):=(\boldsymbol{\beta}^{(ij)})^{\top}\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u}_{i})$
$\displaystyle\quad\forall{i=1,\cdots{N}},\forall{j\in\mathcal{N}_{i}}$ (35)

From (IV-B), $(\boldsymbol{\beta}^{(ij)})$ is obtained by solving the
following problem [33]:

$\displaystyle\underset{\boldsymbol{\beta}^{(ij)}\in\mathbb{R}^{2D}}{\text{
min}}\frac{1}{K}\sum_{m=1}^{K}\mathcal{C}((\boldsymbol{\beta}^{(ij)})^{\top}\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u}_{i}^{(m)}),y_{m})$
$\displaystyle\forall{i=1,\cdots{N}},\forall{j\in\mathcal{N}_{i}}$ (36)

where $\mathbf{u}_{i}^{(m)}$ denotes the $m$-th training sample. Also, we
have:

$\displaystyle\mathcal{C}((\boldsymbol{\beta}^{(ij)})^{\top}\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u}_{i}^{(m)}),y_{m})=\mathcal{L}((\boldsymbol{\beta}^{(ij)})^{\top}\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u}_{i}^{(m)})-y_{m})$
$\displaystyle+\lambda||\boldsymbol{\beta}^{(ij)}||_{2}^{2}$ (37)

where
$||\boldsymbol{\beta}^{(ij)}||_{2}^{2}:=\sum_{n}\sum_{k}\alpha_{n}^{(ij)}\alpha_{k}^{(ij)}\boldsymbol{\Delta}_{\mathbf{W}}^{\top}(\mathbf{u}_{n})\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u}_{k})\simeq||f||_{\mathcal{H}}^{2}$.
The random feature approach can be employed in an online manner where the
nodes estimate the missing values online. In this approach, in each iteration
of the successive approach denoted by $l$, each node solves optimization
problem (IV-B) via online gradient descent [34] to update its parameters which
is expressed as follows:

$\displaystyle\boldsymbol{\beta}^{(ij)}_{l+1}=\boldsymbol{\beta}^{(ij)}_{l}-\eta_{l}\nabla\mathcal{C}((\boldsymbol{\beta}^{(ij)}_{l})^{\top}\boldsymbol{\Delta}_{\mathbf{W}}(\mathbf{u}_{i}),\mathbf{x}_{j}^{l-1}))$
$\displaystyle\forall{i=1,\cdots{N}},\forall{j\in\mathcal{N}_{i}}$ (38)

where $\eta_{l}$ is the step-size of gradient descent. In (IV-B),
$y_{m}=x_{j}^{l-1}$ is used. Then, if the $i$-th node misses the $j$-th node
information ($j\in\mathcal{N}_{i}$), the $i$-th node estimates the missing
value via (IV-B). At the end of $L$ iteration, the complexity of (IV-B) is
$\mathcal{O}({L}D)$.

Please note that for the scenario, where edge fluctuations happen in the first
iteration, the estimators are trained offline for just finding
$\boldsymbol{\beta}^{(ij)}_{1},\forall{j\in\mathcal{N}_{i}}$. Then, the
offline trained parameters are used for initialization. For the rest of
iterations, the parameters of the estimators are updated online via (IV-B).
Finally, the proposed approach in this section can be used for also different
types of graph filters such as FIR [2] or edge-variant graph filters [21].

##  V Conclusion

This paper proposes a new approach to compute the linear transformations as
fast as possible by leveraging a set of successive different shift operators.
To compute the desired transformation by a small number of iterations, a new
optimization problem is proposed to design the shift operators. An exhaustive
simulation study demonstrates that the proposed approach can effectively
compute the desired transformation markedly faster than existing algorithms.
Furthermore, this work studies the effect of edge fluctuations, which is
overlooked by the existing works and a common issue in the wireless sensor
networks due to the random links or node failures. To deal with that issue, an
online learning method is proposed which enables nodes to estimate the missing
values due to the edge fluctuations.

##  Appendix A Proof of Proposition 1

###### Proof.

From (22) and the fact that
$||\mathbf{a}||_{2}^{2}=\mathrm{tr}(\mathbf{a}\mathbf{a}^{\top})$, we have:

$\displaystyle\mathbb{E}[||\boldsymbol{\Omega}||_{2}^{2}]=\mathbb{E}[\mathrm{tr}((\mathbf{z}_{L}+\sum_{i=1}^{L-1}(\Pi_{j=i}^{L-1}\mathbf{S}_{j+1})\mathbf{z}_{i})$
$\displaystyle(\mathbf{z}^{\top}_{L}+\sum_{k=1}^{L-1}\mathbf{z}^{\top}_{k}(\Pi_{m=k}^{L-1}\mathbf{S}_{m+1})^{\top})]$
(39)

where $\boldsymbol{\Omega}=\mathbf{y}-\mathbf{H}\mathbf{x}$. By expanding (A),
we have:

$\displaystyle\mathbb{E}[||\boldsymbol{\Omega}||_{2}^{2}]=\mathbb{E}[\mathrm{tr}(\boldsymbol{\Psi}+\sum_{i=1}^{L-1}(\Pi_{j=i}^{L-1}\mathbf{S}_{j+1})\mathbf{z}_{i}\mathbf{z}_{i}^{\top}(\Pi_{m=i}^{L-1}\mathbf{S}_{m+1})^{\top}+$
$\displaystyle\mathbf{z}_{L}\mathbf{z}_{L}^{\top})]$ (40)

where
$\boldsymbol{\Psi}=\sum_{i=1}^{L-1}(\Pi_{j=i}^{L-1}\mathbf{S}_{j+1})\mathbf{z}_{i}\sum_{k=1,k\neq{i}}^{L-1}\mathbf{z}_{k}^{\top}(\Pi_{m=k}^{L-1}\mathbf{S}_{m+1})^{\top}\\\
+\mathbf{z}_{L}\sum_{k=1}^{L-1}\mathbf{z}_{k}^{\top}(\Pi_{j=k}^{L-1}\mathbf{S}_{j+1})^{\top}+\sum_{i=1}^{L-1}(\Pi_{j=i}^{L-1}\mathbf{S}_{j+1})\mathbf{z}_{i}\mathbf{z}_{L}^{\top}$.

By using the cyclic property of the trace and the fact that trace is a linear
operator, we have:

$\displaystyle\mathbb{E}[||\boldsymbol{\Omega}||_{2}^{2}]=\mathrm{tr}(\mathbb{E}[\mathbf{z}_{L}\mathbf{z}^{\top}_{L}]+\mathbb{E}[\boldsymbol{\Psi}]+\sum_{k=1}^{L-1}\boldsymbol{\Upsilon}_{k}^{\top}\boldsymbol{\Upsilon}_{k}\mathbb{E}[\mathbf{z}_{k}\mathbf{z}_{k}^{\top}])$
(41)

where $\boldsymbol{\Upsilon}_{i}=\Pi_{m=i}^{L-1}\mathbf{S}_{m+1}$. Finally, by
applying the inequality
$\mathrm{tr}(\mathbf{A}\mathbf{B})\leq||\mathbf{A}||_{2}\mathrm{tr}(\mathbf{B})$
(which holds for any square matrix $\mathbf{A}$ and positive semi-definite
matrix $\mathbf{B}$ [35]) and the following properties of the spectral norm
$||\mathbf{A}+\mathbf{B}||_{2}=||\mathbf{A}||_{2}+||\mathbf{B}||_{2}$,
$||\mathbf{A}\mathbf{B}||_{2}\leq||\mathbf{A}||_{2}||\mathbf{B}||_{2}$, we
have

$\displaystyle\mathbb{E}[||\boldsymbol{\Omega}||_{2}^{2}]\leq\mathrm{tr}(\mathbb{E}[\boldsymbol{\Psi}])+2\sum_{j=2}^{L}\rho^{L-j+1}\mathrm{tr}(\mathbb{E}[\mathbf{z}_{j-1}\mathbf{z}_{j-1}^{\top}])+$
$\displaystyle\mathrm{tr}(\mathbb{E}[\mathbf{z}_{L}\mathbf{z}_{L}^{\top}])$
(42)

From
$\mathbb{E}(\mathbf{z}_{j-1}\mathbf{z}_{j-1}^{\top})=\boldsymbol{\Sigma}_{z_{j-1}}+\mathbf{m}_{z_{j-1}}$
and (A), we have:

$\displaystyle\mathbb{E}[||\boldsymbol{\Omega}||_{2}^{2}]\leq\mathbb{E}[\mathrm{tr}(\boldsymbol{\Psi})]+2\sum_{j=2}^{L}\rho^{L-j+1}\mathrm{tr}(\boldsymbol{\Sigma}_{z_{j-1}}+\mathbf{m}_{z_{j-1}})+$
$\displaystyle\mathrm{tr}(\boldsymbol{\Sigma}_{z_{L}}+\mathbf{m}_{z_{L}})$
(43)

∎

## References

  * [1] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” IEEE Signal Processing Magazine, vol. 30, pp. 83–98, May 2013\.
  * [2] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs:frequency analysis,” IEEE Trans. on Signal Processing, vol. 61, pp. 1644–1656, Apr 2013.
  * [3] S. Chen, R. Varma, A. Sandryhaila, and J. Kovačević, “Discrete signal processing on graphs: Sampling theory,” IEEE Transactions on Signal Processing, vol. 63, no. 24, pp. 6510–6523, 2015.
  * [4] S. Safavi and U. A. Khan, “Revisiting finite-time distributed algorithms via successive nulling of eigenvalues,” vol. 22, pp. 54–57, Jan. 2015.
  * [5] T. Weerasinghe, D. Romero, C. Asensio-Marco, and B. Beferull-Lozano, “Fast decentralized subspace projection via graph filtering,” in International Conference on Acoustics, Speech, and Signal Processing (ICASSP), (Calgary,Canada), pp. 4639–4643, Apr. 2018.
  * [6] S. Segarra, A. G. Marques, G. Leus, and A. Ribeiro, “Reconstruction of graph signals through percolation from seeding nodes,” IEEE Transactions on Signal Processing, vol. 64, no. 16, pp. 4363–4378, 2016.
  * [7] S. Barbarossa, G. Scutari, and T. Battisti, “Distributed signal subspace projection algorithms with maximum convergence rate for sensor networks with topological constraints,” in International Conference on Acoustics, Speech, and Signal Processing (ICASSP), (Taipei,Taiwan), pp. 2893–2896, Apr. 2009.
  * [8] P. Di Lorenzo, S. Barbarossa, and S. Sardellitti, “Distributed signal processing and optimization based on in-network subspace projections,” IEEE Transactions on Signal Processing, vol. 68, pp. 2061–2076, 2020.
  * [9] C. A.-M. F. Camaro-Nogues, D. Alonso-Roman and B. Beferull-Lozano, “Reducing the observation error in a wsn through a consensus-based subspace projection,” in WCNC, (Shanghai,China), pp. 3643–3648, Apr. 2013.
  * [10] A. Y. Kibangou, “Graph laplacian based matrix design for finite-time distributed average consensus,” in 2012 American Control Conference (ACC), pp. 1901–1906, 2012.
  * [11] A. Sandryhaila, S. Kar, and J. M. F. Moura, “Finite-time distributed consensus through graph filters,” in International Conference on Acoustics, Speech, and Signal Processing (ICASSP), (Florence, Italy), pp. 1080–1084, May 2014.
  * [12] D. I. Shuman, B. Ricaud, and P. Vandergheynst, “Vertex-frequency analysis on graphs,” Applied and Computational Harmonic Analysis, vol. 40, no. 2, pp. 260–291, 2016.
  * [13] S. K. Narang and A. Ortega, “Perfect reconstruction two-channel wavelet filter banks for graph structured data,” IEEE Transactions on Signal Processing, vol. 60, no. 6, pp. 2786–2799, 2012.
  * [14] B. Girault, Signal Processing on Graphs - Contributions to an Emerging Field.  PhD thesis, École Normale Supérieure de Lyon, 2015.
  * [15] D. Romero, M. Ma, and G. B. Giannakis, “Kernel-based reconstruction of graph signals,” IEEE Transactions on Signal Processing, vol. 65, no. 3, pp. 764–778, 2017.
  * [16] M. Onuki, S. Ono, M. Yamagishi, and Y. Tanaka, “Graph signal denoising via trilateral filter on graph spectral domain,” IEEE Transactions on Signal and Information Processing over Networks, vol. 2, no. 2, pp. 137–148, 2016.
  * [17] S. Segarra, A. G. Marques, and A. Ribeiro, “Optimal graph-filter design and applications to distributed linear network operators,” IEEE Trans. Sig. Process., vol. 65, pp. 4117–4131, May 2017.
  * [18] S. Mollaebrahim, C. Asensio-Marco, D. Romero, and B. Beferull-Lozano, “Decentralized subspace projection in large networks,” in Global Conference on Signal and Information Processing (GlobalSIP), (Anaheim,USA), pp. 788–792, Feb. 2018.
  * [19] M. Zamalloa and B. Krishnamachari, “An analysis of unreliability and asymmetry in low-power wireless links,” ACM Trans. Sen. Netw., vol. 3, p. 1–34, June 2007.
  * [20] L. Sang, A. Arora, and H. Zhang, “On link asymmetry and one-way estimation in wireless sensor networks,” ACM Trans. Sen. Netw., vol. 6, p. 12–25, Mar. 2010.
  * [21] M. Coutino, E. Isufi, and G. Leus, “Advances in distributed graph filtering,” IEEE Transactions on Signal Processing, vol. 67, pp. 2320–2333, May. 2019.
  * [22] E. Isufi and G. Leus, “Distributed sparsified graph filters for denoising and diffusion tasks,” in 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 5865–5869, 2017.
  * [23] E. Isufi, A. Loukas, A. Simonetto, and G. Leus, “Filtering random graph processes over random time-varying graphs,” IEEE Transactions on Signal Processing, vol. 65, no. 16, pp. 4406–4421, 2017.
  * [24] F. Gama, E. Isufi, A. Ribeiro, and G. Leus, “Controllability of bandlimited graph processes over random time varying graphs,” IEEE Transactions on Signal Processing, vol. 67, no. 24, pp. 6440–6454, 2019.
  * [25] L. B. Saad and B. Beferull-Lozano, “Accurate graph filtering in wireless sensor networks,” IEEE Internet of Things Journal, pp. 1–1, 2020.
  * [26] S. Segarra, A. G. Marques, G. Mateos, and A. Ribeiro, “Network topology inference from spectral templates,” IEEE Trans. on Signal and Information Processing over Networks,, vol. 3, p. 467–483, Sep 2017.
  * [27] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs: Frequency analysis,” IEEE Trans. on Signal Processing, vol. 62, pp. 3042–3054, Jun 2014.
  * [28] Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” SIAM Journal on Imaging Sciences, vol. 6, no. 3, pp. 1758–1789, 2013.
  * [29] V. N. Vapnik, Statistical Learning Theory.  Wiley-Interscience, 1998.
  * [30] B. Scholkopf and A. J. Smola, Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond.  MIT Press, 2001.
  * [31] B. Schölkopf, R. Herbrich, and A. J. Smola, “A generalized representer theorem,” in Proceedings of the 14th Annual Conference on Computational Learning Theory, p. 416–426, 2001.
  * [32] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” in Advances in Neural Information Processing Systems 20, pp. 1177–1184, 2008\.
  * [33] Y. Shen, T. Chen, and G. B. Giannakis, “Random feature-based online multi-kernel learning in environments with unknown dynamics,” J. Mach. Learn. Res., vol. 20, p. 773–808, Jan. 2019.
  * [34] E. Hazan, “Introduction to online convex optimization,” Found. Trends Optim., vol. 2, p. 157–325, Aug. 2016.
  * [35] Sheng-De Wang, Te-Son Kuo, and Chen-Fa Hsu, “Trace bounds on the solution of the algebraic matrix riccati and lyapunov equation,” IEEE Transactions on Automatic Control, vol. 31, no. 7, pp. 654–656, 1986.


[datatype=bibtex, overwrite=true] [fieldsource=booktitle,
match=.*Interspeech.*, replace=Proc. INTERSPEECH] [fieldsource=journal,
match=.*INTERSPEECH.*, replace=Proc. INTERSPEECH] [fieldsource=booktitle,
match=.*ICASSP.*, replace=Proc. ICASSP] [fieldsource=booktitle,
match=.*icassp_inpress.*, replace=ICASSP (in press)] [fieldsource=booktitle,
match=.*International.*Conference.*on.*Acoustics,.*Speech.*and.*Signal.*Processing.*,
replace=ICASSP] [fieldsource=booktitle,
match=.*International.*Conference.*on.*Learning.*Representations.*,
replace=ICLR] [fieldsource=booktitle,
match=.*International.*Conference.*on.*Machine.*Learning.*, replace=Proc.
ICML] [fieldsource=booktitle,
match=.*Automatic.*Speech.*Recognition.*and.*Understanding.*, replace=Proc.
ASRU] [fieldsource=booktitle, match=.*Spoken.*Language.*Technology.*,
replace=Proc. SLT] [fieldsource=booktitle,
match=.*Speech.*Synthesis.*Workshop.*, replace=Proc. SSW]
[fieldsource=booktitle, match=.*workshop.*on.*speech.*synthesis.*,
replace=Proc. SSW] [fieldsource=booktitle,
match=.*Advances.*in.*neural.*information.*processing.*, replace=Proc. NIPS]
[fieldsource=booktitle,
match=.*Advances.*in.*Neural.*Information.*Processing.*, replace=Proc. NIPS]
[fieldsource=booktitle, match=.*AAAI.*, replace=Proc. AAAI]
[fieldsource=booktitle, match=.*Workshop.*on.* Applications.* of.*
Signal.*Processing.*to.*Audio.*and.*Acoustics.*, replace=Proc. WASPAA]
[fieldsource=booktitle,
match=.*International.*Conference.*on.*Language.*Resources.*and.*Evaluation.*,
replace=Proc. LREC] [fieldsource=journal,
match=.*Spontaneous.*Speech.*Processing.*and.*Recognition, replace=Proc. SSPR]
[fieldsource=publisher, match=.+, replace=] [fieldsource=month, match=.+,
replace=] [fieldsource=location, match=.+, replace=] [fieldsource=address,
match=.+, replace=] [fieldsource=organization, match=.+, replace=]
[fieldsource=doi, match=.+, replace=] [fieldsource=url, match=.+, replace=]
[fieldsource=editor, match=.+, replace=] The University of Tokyo, Japan.

# Diversity-based core-set selection for text-to-speech
with linguistic and acoustic features

###### Abstract

This paper proposes a method for extracting a lightweight subset from a text-
to-speech (TTS) corpus ensuring synthetic speech quality. In recent years,
methods have been proposed for constructing large-scale TTS corpora by
collecting diverse data from massive sources such as audiobooks and YouTube.
Although these methods have gained significant attention for enhancing the
expressive capabilities of TTS systems, they often prioritize collecting vast
amounts of data without considering practical constraints like storage
capacity and computation time in training, which limits the available data
quantity. Consequently, the need arises to efficiently collect data within
these volume constraints. To address this, we propose a method for selecting
the core subset (known as core-set) from a TTS corpus on the basis of a
diversity metric, which measures the degree to which a subset encompasses a
wide range. Experimental results demonstrate that our proposed method performs
significantly better than the baseline phoneme-balanced data selection across
language and corpus size.

Index Terms—  text-to-speech synthesis, data selection, core-set selection,
corpus construction, diversification

##  1 Introduction

Although text-to-speech (TTS) has achieved human-level naturalness in
transforming text into speech waveform [1], its expressive capabilities do not
yet match those of humans. Previous studies have aimed to enhance TTS
expressiveness by addressing aspects such as speaker identity control [2],
emotional expression [3], and prosody representation [4]. These studies
predominantly adopt data-driven methods based on machine learning, and corpora
with wider speaker or style variation are increasingly being anticipated.

To construct TTS corpora containing diverse data, previous studies have
gathered data from vast sources such as audiobooks [5] and YouTube [6]. These
methodologies are predicated on the belief that collecting data from large-
scale sources inherently results in a diverse dataset, and they utilize all
available data. Consequently, recently released speech corpora [7, 8] comprise
several thousand or even tens of thousands of hours of data. In practice,
however, limitations on storage capacity and the time required for learning
impose constraints on the amount of available data, requiring datasets to be
efficiently collected within these constraints [9, 10]. Given this context,
TTS training in practical environments is expected to be achieved if data size
can be reduced without compromising the quality of synthetic speech.

As a relevant machine learning technique, core-set selection has been proposed
[11]. This method aims to extract a subset (core-set) that achieves an
equivalent learning effect as the entire dataset, as shown in Fig. 1. Unlike
point-wise data selection, which considers each data point independently, this
method is a kind of subset selection and considers a subset as a whole when
selecting data. Furthermore, as shown in Fig. 1, it seeks to cover the entire
range, rather than just extracting a specific region.

This paper proposes a core-set selection method for multi-speaker TTS. We
define a diversity metric based on language and speech features derived from
self-supervised learning (SSL) models to assess the coverage area of a subset,
formulating the core-set selection as a diversity maximization task under
constraints on subset size. Our proposed method is computationally
lightweight. Specifically, it does not involve any model training with the
entire dataset. We conduct experiments on TTS corpora in Japanese, Chinese,
and English to demonstrate that our core-set selection method mitigates the
degradation in the naturalness and intelligibility of synthetic speech
compared with phoneme balance based subset selection [12]. The contributions
of this work are as follows:

  * •

This paper is the first to introduce core-set selection in TTS tasks.

  * •

Our proposed method is computationally efficient, mitigating degradation of
naturalness and intelligibility in synthetic speech.

  * •

We conduct experiments with multiple languages and dataset sizes to
demonstrate the validity of our proposed method.

(a) The core-set ideally has a learning effect equivalent to the entire
dataset.

(b) The core-set aims to cover an entire range, not a specific region.

Fig. 1: The purpose and characteristic of core-set selection.

##  2 Related work

###  2.1 Designing phoneme-balanced speech corpora

A classical method for constructing a balanced speech dataset is to construct
a phoneme-balanced sentence set [12] by maximizing the following function
defined to evaluate the phoneme balance:

$H(\bm{p})=-\sum_{i=1}^{n}p_{i}\log p_{i}$ (1)

where $p_{i}$ is the occurrence probabilities of the $i$-th phoneme and $n$ is
the number of phonemes. This method aims to enhance the effectiveness of model
training by avoiding situations where a specific phoneme occurs extremely
infrequently.

###  2.2 Data selection for multi-speaker TTS

Previous studies [13, 6] proposed data selection methods for multi-speaker
TTS, but they were point-wise data selection methods whereas this study
addresses a subset selection method. A previous study [14] proposed a data
subset selection method based on speaker selection, but that method extracts a
specific region from the speaker distribution, whereas the subset selection
method in this study aims to cover the entire range of the original dataset.
Previous studies [15, 16] proposed text set construction methods by clustering
in the linguistic embedding space, and another [17] demonstrated a selection
method based on phoneme and prosody entropy that is effective for statistical
parametric speech synthesis. However, they do not consider other aspects such
as speaker identity, important factors for multi-speaker TTS. In other words,
core-set selection methods with the goal of efficiently training large-scale
multi-speaker TTS models have not been explored.

###  2.3 Diversity-based subset selection algorithms

A previous study [11] proposed a diversity-based core-set selection using the
$k$-center method [18], and that study uses a similar algorithm. Diversity-
based subset selection is applied in various domains including recommendation
systems [19], and several diversity evaluation metrics have been proposed. One
such metric [20] is used to drive the $k$-center algorithm. We adopt a
modified version of this metric [21], where $\mathrm{max}$ is replaced with
$\mathrm{sum}$.

##  3 Proposed method

We first extract utterance-level feature vectors that encompass linguistic,
speaker, and acoustic features for each data point. Using these feature
vectors, we define a diversity metric and use the proposed core-set selection
to maximize the diversity metric.

###  3.1 Feature extraction

For diversity-based core-set selection, we use feature vectors where similar
data appear close together and dissimilar data appear far apart in feature
space. Each data point in multi-speaker TTS corpora consists of a pair of
text, speaker identity, and speech, and we utilize the joint features of each
aspect.

Linguistic features $\bm{x}_{\text{linguistic}}$: We use sentence embeddings
of texts, which are expected to separate data from different text domains.
Specifically, we average the output vector sequences from BERT [22] to obtain
fixed-dimensional vectors, as detailed in [23], and then normalize them to
have a norm equal to $1$. Although this method may lead to the loss of
specific features of individual words, it is still effective to place similar
data points in close proximity in the feature space.

Speaker features $\bm{x}_{\text{speaker}}$: We use continuous speaker
representations like $x$-vectors [24] to incorporate speaker similarity into
similarity calculations for the data. We normalize $x$-vectors to have a norm
equal to $1$. For multi-style TTS, our proposed method can be applied by using
a pre-trained style encoder to extract style features.

Acoustic features $\bm{x}_{\text{acoustic}}$: We average the output vector
sequences of wav2vec 2.0 [25] and normalize their norms. SSL features have
demonstrated their effectiveness in speech recognition [25], and phonetic
information is believed to be represented in their frame-level features.
Averaging these features is analogous to using dense vectors instead of one-
hot vectors in time-frequency representations and is expected to contribute to
expanding the phonetic coverage of the core-set.

Since linguistic and speaker features correspond to the input to the TTS
model, we define input features $\bm{x}_{\text{ipnut}}$ by concatenating them,
and define output features $\bm{x}_{\text{output}}$ as the acoustic feature.
Finally, we concatenate input and output features to obtain the joint features
$\bm{x}_{\text{joint}}$ used for calculating diversity. The relationship
between the features is described in Eq. (2).

$\bm{x}_{\text{input}}=\begin{bmatrix}\bm{x}_{\text{linguistic}}\\\
\bm{x}_{\text{speaker}}\end{bmatrix},\;\bm{x}_{\text{output}}=\bm{x}_{\text{acoustic}},\;\bm{x}_{\text{joint}}=\begin{bmatrix}\bm{x}_{\text{input}}\\\
\bm{x}_{\text{output}}\end{bmatrix}.$ (2)

Algorithm 1 Select a core-set $S$ from dataset $D$

1: $S\leftarrow\varnothing$

2: $\bm{x}\sim U(D)$

3: while $T(S\cup\\{\bm{x}\\})\leq t_{\mathrm{max}}$ do

4: $S\leftarrow S\cup\\{\bm{x}\\}$

5: $\bm{x}\leftarrow\textrm{argmax}_{\bm{x}\in D\setminus S}\sum_{\bm{y}\in
S}\|\bm{x}-\bm{y}\|^{2}$

6: end while

###  3.2 Diversity evaluation metric

Let $S$ and $D$ respectively denote a data subset and the entire dataset,
$\bm{x},\bm{y}\in\mathbb{R}^{d}$ denote feature vectors, and
$\|\cdot\|,\mathrm{cossim}$ represent $L^{2}$ norm and cosine similarity. To
assess the diversity of $S$, previous studies [21] proposed calculating
$\sum_{\bm{x},\bm{y}\in S}d(\bm{x},\bm{y})$ where $d(\cdot,\cdot)$ is
dissimilarity (e.g., $\|\cdot\|$ or $-\mathrm{cossim}$). Since $V(S)$ takes a
higher value when $S$ contains many dissimilar data pairs, $V(S)$ is
anticipated to represent the spread of $S$ in the feature space. When the
norms of $\bm{x}$ and $\bm{y}$ are equal to $1$, the relationship
$\|\bm{x}-\bm{y}\|^{2}=2-2\mathrm{cossim}(\bm{x},\bm{y})$ holds and squared
distance has a close relationship with cosine similarity. Therefore, we adopt
squared Euclidean distance as a dissimilarity metric and evaluate diversity
using the following function $V(S)$, the same as a previous study [26]

$V(S):=\sum_{x,y\in S}\|\bm{x}-\bm{y}\|^{2}.$ (3)

###  3.3 Core-set selection algorithm

We conduct core-set selection by solving the optimization problem that
involves maximizing the diversity score $V(S)$ for subset $S$ subject to size
constraints about $S$. However, this optimization problem is a combinatorial
optimization and can lead to explosive computational complexity. Since this
study focuses on scenarios with a large corpus and requires algorithms with
low computational resources, we utilize a greedy algorithm. Specifically, we
execute the core-set selection procedure by sequentially adding the data that
maximizes the diversity score until we reach the desired core-set size.

When adding each data, we select $\bm{x}$ to maximize $V(S\cup\\{\bm{x}\\})$,
which can be expressed as the sum of $V(S)$ and $2\times\sum_{\bm{y}\in
S}\|\bm{x}-\bm{y}\|^{2}$. Since $V(S)$ does not depend on $\bm{x}$, the core-
set selection procedure follows the algorithm outlined in Algorithm 1, where
$U(D)$ represents a uniform distribution on $D$, $T(S)$ represents total
speech duration included in $S$, and $t_{\mathrm{max}}$ represents a
constraint on $T(S)$. Notably, this algorithm can be executed with feature
vectors in instead of the actual data, leading to reduced storage
requirements.

##  4 Experimental evaluation

###  4.1 Experimental conditions

####  4.1.1 Dataset

We trained a monolingual multi-speaker TTS model using a multi-speaker TTS
corpus for each of Japanese, Chinese, and English. We used parallel 100 and
nonparallel 30 subsets from the JVS [27] corpus for Japanese, the AISHELL-3
[28] corpus for Chinese, and the training sets from the LibriTTS-R [29] corpus
(train_clean_100 and train_clean_360) for English. The corpus sizes are
$25$-hour, $63$-hour, and $243$-hour, respectively, from which the core-sets
are selected. Each corpus includes $100$, $174$, and $1151$ speakers.

For text datasets for evaluating Japanese TTS models, we used $324$ sentences
from the ITA corpus [30]. For the other languages, we randomly selected $100$
sentences from the test set of the corpora.

Table 1: Model-wise pseudo-MOS for TTS models trained with each dataset. The values in parentheses represent the core-set in hours. Bold indicates the highest value among the subsets. Dataset | All | Phoneme | Input | Our
---|---|---|---|---
Data | Balance | Balance | method
JVS (3h) | $3.020$ | $2.971$ | $2.942$ | $\mathbf{2.996}$
JVS (6h) | $3.020$ | $2.966$ | $3.006$ | $\mathbf{3.080}$
JVS (12h) | $3.020$ | $3.022$ | $2.996$ | $\mathbf{3.080}$
AISHELL-3 (6h) | $2.604$ | $2.717$ | $2.724$ | $\mathbf{2.796}$
LibriTTS-R (25h) | $3.943$ | $3.875$ | $3.891$ | $\mathbf{3.902}$
Table 2: ASR error rate $(\%)$ in transcribing synthesized speech. A lower value implies higher intelligibility of synthesized speech. Bold indicates the lowest value among the subsets. Dataset | All | Phoneme | Input | Our
---|---|---|---|---
Data | Balanced | Balanced | method
JVS (3h) | $19.92$ | $24.53$ | $24.25$ | $\mathbf{21.26}$
JVS (6h) | $19.92$ | $23.37$ | $22.42$ | $\mathbf{20.25}$
JVS (12h) | $19.92$ | $21.19$ | $21.91$ | $\mathbf{20.37}$
AISHELL-3 (6h) | $15.55$ | $18.09$ | $18.93$ | $\mathbf{17.97}$
LibriTTS-R (25h) | $16.61$ | $18.22$ | $\mathbf{17.54}$ | $17.84$

####  4.1.2 Model and training

The multi-speaker TTS models included FastSpeech 2 [31] and the pre-trained
HiFi-GAN vocoder [32] UNIVERSAL_V1 [33]. We followed hyperparameters of the
open-source implementations [34, 35]. For speaker representation, we opted for
$512$-dimensional $x$-vector, using a pre-trained model [36]. Each unique
speaker corresponded to one $x$-vector, with all $x$-vectors having an L2 norm
equal to $1$. The $x$-vector was added to the output of the FastSpeech 2
encoder via a $512$-by-$256$ linear layer. The number of training steps was
set in accordance with the size of each corpus: $50k$ steps for JVS and
AISHELL-3, and $300k$ steps for LibriTTS-R.

####  4.1.3 Feature extractors for core-set selection

We used language-specific BERT models [37, 38, 39] for linguistic features,
$x$-vectors for speaker features, and wav2vec 2.0 [25] base model [40] pre-
trained on Librispeech [41] for acoustic features.

####  4.1.4 Compared methods

We compared the following data selection methods.

All Data: To assess the quality degradation in training with subsets, we
conducted training using the entire dataset.

Phoneme Balance: As a conventional method, a subset was selected to maximize
phoneme entropy as described in Sec. 2.1.

Input Balance: Expanding phoneme balance to multi-speaker scenarios, we used
subset selection by maximizing the sum of phoneme entropy and speaker ID
entropy, under the expectation of enhancing the learning effect by reducing
speaker imbalance.

Our method: Our diversity-based core-set selection. We calculated the joint
feature vector for each data point and then incrementally added data to the
core-set, maximizing the diversity score.

####  4.1.5 Comparison conditions

We conducted experiments to answer the following questions:

Does our method work?: validation with varying core-set size. To validate
whether our method is more effective than traditional balance-based methods,
experiments were conducted with multiple core-set size. Core-sets of
approximately $10,20,40\%$ of JVS were evaluated, corresponding to $3,6,12$
hours, respectively.

Does our method work across language and corpus size?: validation on multiple
corpora. We compared the balance-based methods and our method across varying
languages and corpus sizes, specifically using AISHELL-3 and LibriTTS-R. We
selected $6$-hour and $25$-hour core sets, which correspond to about $10\%$ of
the corpus.

Are joint features effective?: ablation study about features. To assess the
effectiveness of combination of input and output features, we conducted core-
set selection with each feature. Core-set sizes were set to $3,6,12$ hours.

####  4.1.6 Evaluation criteria

We synthesized speech for all speakers included in each corpus with test
sentences prepared in Sec. 4.1.1 and evaluate them using both automated and
human subjective assessments.

For automatic evaluation, we used pseudo-MOS, an automatically predicted mean
opinion score (MOS) of synthetic speech. Specifically, we used the UTMOS [42]
strong learner model [43], which has high accuracy in English, Chinese [44],
and Japanese [6].

Model-wise pseudo-MOS evaluation: To quantitatively compare the overall
performance of multi-speaker TTS models, pseudo-MOSs were averaged over all
speakers for each model.

Speaker-wise pseudo-MOS evaluation: To assess speaker-wise performance, the
pseudo-MOSs were averaged per speaker.

We also calculated recognition error rates by using an automatic speech
recognition (ASR) model, Whisper [45] large model. We evaluated character
error rate for Japanese, phoneme error rate in Chinese, and word error rate in
English, which are referred to as ASR error rate.

As a subjective evaluation experiment, we conducted MOS tests on speech
naturalness and aggregated speaker-wise and model-wise MOS. The evaluation for
JVS encompassed all data, both balance methods for the $3$-hour core-set, and
our method for the $3$, $6$, and $12$-hour core-set. For the other corpora,
all methods were included. There were $800$ evaluators for JVS and $400$ for
every other corpus. Each listener assessed $24$ samples using a 5-point scale.
Since LibriTTS-R has a large number of speakers, we sampled speakers at
intervals of $10$ in pseudo-MOS order and evaluated $116$ speakers. We
evaluated all speakers for the other corpora.

Furthermore, in an ablation study about features, we conducted subjective
preference tests comparing each feature with joint features, using $3$-hour
core-sets. For each comparison, $100$ evaluators assessed the naturalness of
synthesized speech for 10 randomly selected combinations of speakers and
sentences.

Fig. 2: Cumulative histogram of speaker-wise pseudo-MOS for each model. The
values on the $y$-axis represent the number of speakers with pseudo-MOS scores
higher than the values on the $x$-axis.

###  4.2 Results

####  4.2.1 Validation with varying core-set size

The upper half of Table 1 presents the model-wise pseudo-MOSs for the models
trained on core-sets from JVS. Our method consistently outperforms the other
balance methods across all core-set sizes, with an average improvement of
$0.068$. Considering a previous study [6] that demonstrated that pseudo-MOS
values exhibit approximately half the range of MOS variation in Japanese, we
can expect a wider range of improvement in MOS. Additionally, in the speaker-
wise pseudo-MOS evaluations shown in the upper part of Figure 2, our method
consistently achieves higher scores than balance methods across most ranges.
Furthermore, our method has better ASR error rates than the other balance
methods among all subset sizes, as shown in the upper half of Table 2. These
results demonstrate that our method works to mitigate the decrease in
naturalness and intelligibility of synthetic speech regardless of core-set
size.

(a) Comparison with the JVS corpus.

(b) Comparison with the AISHELL-3 corpus and the LibriTTS-R corpus.

Fig. 3: Cumulative histogram of speaker-wise MOS for each model.

Figure 3 illustrates the results of the subjective evaluation. Within the
$3$-hour core-sets, our method consistently produces higher curves than the
balance methods, indicating that our method improves naturalness of synthetic
speech for all speakers. The model-wise MOSs for phoneme balance, input
balance, and our method were $3.069$, $2.999$, and $3.217$, respectively.
Consequently, our method achieved an average improvement of $0.183$ in model-
wise MOS, demonstrating its effectiveness in terms of naturalness in
subjective evaluations.

Notably, the model-wise MOS with the $12$-hour core-set selected by our method
is $3.614$, which closely matches the $3.587$ achieved by all data. Also, this
core-set achieved nearly equal values to all data in terms of pseudo-MOS and
ASR error rates (see Tables 1,2). These results indicate that the core-set
attained an equivalent learning effect to that of the entire dataset in terms
of naturalness and intelligibility. In other words, it highlights the validity
of using core-set selection as advantageous in terms of data volume.

####  4.2.2 Validation on multiple corpora

The lower half of Table 1 displays model-wise pseudo-MOSs for Chinese and
English. In both AISHELL-3 and LibriTTS-R, our method outperforms the balance-
based methods, demonstrating an average improvement of $0.101$ and $0.018$,
respectively. Furthermore, in terms of the increment of ASR error rates (shown
in the lower half of Table 2) compared with all data, the input balance
exhibits $3.38\%$ in AISHELL-3 while the proposed method reduces to $2.42\%$.
In the case of LibriTTS-R, the phoneme balance method results in $1.61\%$,
whereas the proposed method reduces to $1.23\%$. From these results, we can
say the TTS model trained with the core-set selected by the proposed method
can mitigate degradation in naturalness and intelligibility better than the
balance methods.

The lower half of Figure 2 shows the results of speaker-wise pseudo-MOS in
Chinese and English. For AISHELL-3, our method’s curve is shifted more to the
right than the others, clearly indicating its superiority over the other
balanced methods. The lower performance of all data is attributed to its
imbalance, e.g., phonemic imbalance; the phonemic entropy (Eq. (1)) dropped
from $6.8$ in the phoneme balance method to $6.5$. For LibriTTS-R, although
the difference between our method and the other balance-based methods is
marginal, the zoomed-in figure (bottom left) reveals that our method has fewer
speakers with low pseudo-MOSs. This implies that, while other methods suffer a
decrease in quality for speakers with less data, our method effectively
addresses and corrects this issue.

The lower half of Figure 3 shows the speaker-wise MOS in Chinese and English.
For AISHELL-3, the model-wise MOS for all data, phoneme balance, input
balance, and our method were $3.627$, $3.679$, $3.634$, and $3.657$,
respectively. In the case of LibriTTS-R, the corresponding scores were
$3.642$, $3.568$, $3.632$ and $3.606$. MOSs for our method are nearly
equivalent to those of all data for all speakers in both corpora, indicating
that the proposed method performs at a level similar to that of the entire
dataset in terms of MOS.

From these results, we conclude our method is applicable and effective
irrespective of language or corpus size.

Table 3: Model-wise pseudo-MOS for each feature. Values in parentheses represent difference from all data. Core-set | All | Input | Output | Joint
---|---|---|---|---
size | Data | features | features | features
$3$-hour | $3.020$ | $2.943({-0.077})$ | $2.931({-0.089})$ | $2.996({-0.024})$
$6$-hour | $3.020$ | $3.032({+0.012})$ | $2.984({-0.036})$ | $3.080({+0.060})$
$12$-hour | $3.020$ | $3.040({+0.020})$ | $3.047({+0.027})$ | $3.080({+0.060})$
Table 4:  Subjective evaluation on the naturalness of synthesized speech. Comparison between individual features and a joint feature. Compared | Score | $p$ value
---|---|---
Compared | vs. | Joint features
Input features | $0.441$ | vs. |  $\mathbf{0.559}$ | $9.82\texttimes 10^{-5}$
Output features | $0.475$ | vs. |  $0.525$ | $1.03\texttimes 10^{-1}$

####  4.2.3 Ablation study about features

Table 3 presents model-wise pseudo-MOSs. Joint features exhibit higher values
than the other features. Particularly within the 3-hour core-sets, where the
other features showed a decrease of $0.077$ and $0.089$, joint features
reduced it to $0.024$. Additionally, Table 4 presents the results of the
subjective evaluation. Although the $p$-value against output features is not
very small, the results suggest that joint features work to reduce the
degradation of synthetic speech naturalness. From these results, we can say
that combining input and output features is suitable for measuring similarity
in the training data.

##  5 Conclusion

We proposed a core-set selection method for multi-speaker text-to-speech
(TTS), which extracts a diverse subset on the basis of language and acoustic
features. Experimental results demonstrated that our proposed method improves
the learning effect compared with phoneme balance based subset selection
across multiple languages, corpora, and core-set sizes. We anticipate that our
proposed method will remain applicable even for larger corpora, thanks to its
low computational and storage requirements. Our future work includes
conducting empirical experiments with even larger corpora.

## References

  * [1] Jonathan Shen et al. “Natural TTS synthesis by conditioning WaveNet on mel spectrogram predictions” In _Proc. ICASSP_ , 2018, pp. 4779–4783
  * [2] Sercan Arik et al. “Neural voice cloning with a few samples” In _Advances in neural information processing systems_ 31, 2018
  * [3] Tao Li et al. “Controllable emotion transfer for end-to-end speech synthesis” In _2021 12th International Symposium on Chinese Spoken Language Processing (ISCSLP)_ , 2021, pp. 1–5
  * [4] RJ Skerry-Ryan et al. “Towards end-to-end prosody transfer for expressive speech synthesis with tacotron” In _international conference on machine learning_ , 2018, pp. 4693–4702
  * [5] Heiga Zen et al. “LibriTTS: A Corpus Derived from LibriSpeech for Text-to-Speech” In _Proc. INTERSPEECH_ , 2019, pp. 1526–1530
  * [6] Kentaro Seki et al. “Text-to-speech synthesis from dark data with evaluation-in-the-loop data selection” In _Proc. ICASSP_ , 2023, pp. 1–5
  * [7] Guoguo Chen et al. “GigaSpeech: An Evolving, Multi-domain ASR Corpus with 10,000 Hours of Transcribed Audio” In _arXiv:2106.06909_ , 2021
  * [8] Shinnosuke Takamichi et al. “JTubeSpeech: corpus of Japanese speech collected from YouTube for speech recognition and speaker verification” In _arXiv:2112.09323_ , 2021
  * [9] Sujeendran Menon et al. “Development of a neural network library for resource constrained speech synthesis” In _2020 5th IEEE International Conference on Recent Advances and Innovations in Engineering (ICRAIE)_ , 2020, pp. 1–8
  * [10] Renqian Luo et al. “Lightspeech: Lightweight and fast text to speech with neural architecture search” In _Proc. ICASSP_ , 2021, pp. 5699–5703
  * [11] Ozan Sener and Silvio Savarese “Active learning for convolutional neural networks: A core-set approach” In _arXiv preprint arXiv:1708.00489_ , 2017
  * [12] Akira Kurematsu et al. “ATR Japanese speech database as a tool of speech recognition and synthesis” In _Speech communication_ 9.4, 1990, pp. 357–363
  * [13] Kai-Zhan Lee and Erica Cooper “A comparison of speaker-based and utterance-based data selection for text-to-speech synthesis” In _Interspeech 2018_ 12873, 2018
  * [14] Pilar Oplustil Gallegos et al. “An Unsupervised Method to Select a Speaker Subset from Large Multi-Speaker Speech Synthesis Datasets.” In _INTERSPEECH_ , 2020, pp. 1758–1762
  * [15] Meysam Shamsi et al. “Corpus Design Using Convolutional Auto-Encoder Embeddings for Audio-Book Synthesis.” In _Proc. INTERSPEECH_ , 2019, pp. 1531–1535
  * [16] Stefan Taubert et al. “A Comparison of Text Selection Algorithms for Sequence-to-Sequence Neural TTS” In _2022 IEEE International Conference on Signal Processing, Communications and Computing (ICSPCC)_ , 2022, pp. 1–6
  * [17] Takashi Nose et al. “Sentence selection based on extended entropy using phonetic and prosodic contexts for statistical parametric speech synthesis” In _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ 25.5, 2017, pp. 1107–1116
  * [18] Teofilo F Gonzalez “Clustering to minimize the maximum intercluster distance” In _Theoretical computer science_ 38, 1985, pp. 293–306
  * [19] Laming Chen, Guoxin Zhang and Eric Zhou “Fast greedy map inference for determinantal point process to improve recommendation diversity” In _Advances in Neural Information Processing Systems_ 31, 2018
  * [20] Jade Goldstein and Jaime Carbonell “The use of mmr and diversity-based reranking in document reranking and summarization” In _Proceedings of the 14th Twente Workshop on Language Technology in Multimedia Information Retrieval_ , 1998, pp. 152–166
  * [21] Sreenivas Gollapudi and Aneesh Sharma “An axiomatic approach for result diversification” In _Proceedings of the 18th international conference on World wide web_ , 2009, pp. 381–390
  * [22] Jacob Devlin Ming-Wei Chang Kenton and Lee Kristina Toutanova “Bert: Pre-training of deep bidirectional transformers for language understanding” In _Proceedings of naacL-HLT_ 1, 2019, pp. 2
  * [23] Nils Reimers and Iryna Gurevych “Sentence-bert: Sentence embeddings using siamese bert-networks” In _arXiv preprint arXiv:1908.10084_ , 2019
  * [24] David Snyder et al. “X-vectors: Robust DNN embeddings for speaker recognition” In _Proc. ICASSP_ , 2018, pp. 5329–5333
  * [25] Alexei Baevski et al. “wav2vec 2.0: A framework for self-supervised learning of speech representations” In _Proc. NeurIPS_ 33, 2020, pp. 12449–12460
  * [26] Zaiqiao Meng and Hong Shen “Scalable aspects learning for intent-aware diversified search on social networks” In _IEEE Access_ 6, 2018, pp. 37124–37137
  * [27] Shinnosuke Takamichi et al. “JVS corpus: free Japanese multi-speaker voice corpus” In _arXiv:1908.06248_ , 2019
  * [28] Yao Shi et al. “Aishell-3: A multi-speaker mandarin tts corpus and the baselines” In _arXiv preprint arXiv:2010.11567_ , 2020
  * [29] Yuma Koizumi et al. “LibriTTS-R: A Restored Multi-Speaker Text-to-Speech Corpus” In _Proc. INTERSPEECH 2023_ , 2023, pp. 5496–5500
  * [30] “ITA corpus”, https://github.com/mmorise/ita-corpus
  * [31] Yi Ren et al. “Fastspeech 2: Fast and high-quality end-to-end text to speech” In _Proc. ICLR_ , 2021
  * [32] Jungil Kong, Jaehyeon Kim and Jaekyoung Bae “HiFi-GAN: Generative adversarial networks for efficient and high fidelity speech synthesis” In _Proc. NeurIPS_ 33, 2020, pp. 17022–17033
  * [33] “HiFi-GAN”, https://github.com/jik876/hifi-gan
  * [34] “FastSpeech2 JSUT implementation”, https://github.com/Wataru-Nakata/FastSpeech2-JSUT
  * [35] “FastSpeech 2 - PyTorch Implementation”, https://github.com/ming024/FastSpeech2
  * [36] “x-vector extractor for Japanese speech”, https://github.com/sarulab-speech/xvector_jtubespeech
  * [37] Kobayashi ”Koga et al. “LINE DistilBERT Japanese”, http://github.com/line/LINE-DistilBERT-Japanese, 2023
  * [38] “bert-base-chinese”, https://huggingface.co/bert-base-chinese
  * [39] “distilbert-base-uncased”, https://huggingface.co/distilbert-base-uncased
  * [40] “wav2vec 2.0”, https://github.com/facebookresearch/fairseq/tree/main/examples/wav2vec
  * [41] Vassil Panayotov et al. “Librispeech: an asr corpus based on public domain audio books” In _Proc. ICASSP_ , 2015, pp. 5206–5210
  * [42] Takaaki Saeki et al. “UTMOS: UTokyo-SaruLab System for VoiceMOS Challenge 2022” In _Proc. INTERSPEECH_ , 2022, pp. 4521–4525
  * [43] “UTMOS: UTokyo-SaruLab MOS Prediction System”, https://github.com/sarulab-speech/UTMOS22
  * [44] Wen Chin Huang et al. “The VoiceMOS Challenge 2022” In _Proc. INTERSPEECH_ , 2022, pp. 4536–4540
  * [45] Alec Radford et al. “Robust speech recognition via large-scale weak supervision” In _Proc. ICML_ , 2023, pp. 28492–28518