# Topological Estimation of Number of Sources
in Linear Monocomponent Mixtures

Sean Kennedy  Electrical and Computer Engineering
Naval Postgraduate School
Monterey, CA, USA
<EMAIL_ADDRESS>Murali Tummala  Electrical and Computer Engineering
Naval Postgraduate School
Monterey, CA, USA
<EMAIL_ADDRESS>John McEachen  Electrical and Computer Engineering
Naval Postgraduate School
Monterey, CA, USA
<EMAIL_ADDRESS>
###### Abstract

Estimation of the number of sources in a linear mixture is a critical
preprocessing step in the separation and analysis of the sources for many
applications. Historically, statistical methods, such as the minimum
description length and Akaike information criterion, have been used to
estimate the number of sources based on the autocorrelation matrix of the
received mixture. In this paper, we introduce an alternative, topology-based
method to compute the number of source signals present in a linear mixture for
the class of constant-amplitude, monocomponent source signals. As a proof-of-
concept, we include an example of three such source signals that overlap at
multiple points in time and frequency, which the method correctly identifies
from a set of eight redundant measurements. These preliminary results are
promising and encourage further investigation into applications of topological
data analysis to signal processing problems.

###### Index Terms:

Number of sources, persistent homology, embedding, monocomponent, array

##  I Introduction

The objective of Blind Source Separation (BSS), also called blind signal
separation or source separation, is to separate a group of source signals from
a mixture without detailed knowledge of the sources or mixing process itself.
BSS has applicability to many radar, communication, and imaging scenarios.
Numerous techniques have been developed for BSS of mixed signals [1]. However,
in most cases, the separation techniques assume that the number of sources
equals the number of observations or that the number of sources is known in
advance [1, 2, 3]. Thus, the ability to estimate or determine the number of
sources to be unmixed is a critical preprocessing step in practical
implementation.

The problem of estimating the number of sources in a linear mixture has been
studied in the literature [4, 5, 6, 7]. The two most popular methods used for
this estimation are information-theoric measures based on the statistics of
the mixture’s autocorrelation matrix: the Minimum Description Length (MDL) and
Akaike Information Criterion (AIC) estimators [5, 6]. Both have been shown to
work well under the assumption of temporally and spatially white noise and
Gaussian-distributed random sources. However, they have also been shown to be
non-robust when real-world data deviates from these source and noise models
[6]. While numerous modifications and enhancements to these methods have been
proposed, e.g., [5, 6, 7], the fundamental underpinning of these methods
remains rooted in the analysis of eigenvalues and eigenvectors of the sample
autocorrelation matrix with the assumption that there exist fewer sources than
mixed measurements.

In this paper, we introduce an approach that departs from the statistical
inference model for estimating the number of sources. Instead, we frame the
problem in terms of topology and show that existing Topological Data Analysis
(TDA) tools can be used to estimate the number of independent sources under
certain scenarios. The method is mathematically motivated and developed in
Section II. In Section III, we provide a simple validation of the method for
three nonstationary sources. Conclusions and a discussion of future research
goals are provided in Section IV.

##  II Method Development and Discussion

In practical terms, the method consists of only three steps: (1) embed an
observed signal as a manifold in higher dimensional space, (2) compute the
Betti number sequence [8] of the manifold using TDA, and (3) match the Betti
number sequence to a known reference sequence. As shown in Section III,
implementation of these steps is straightforward with appropriate software.
Thus, the primary contribution of this paper is the analysis of the
mathematical mechanisms by which the method achieves valid results. In
Subsection II-A we provide the basic topological theory behind the method. In
Subsection II-B we apply this theory to the case of a linear mixture as might
be encountered in a radar, sonar, or multi-receiver communication array. In
Subsection II-C, we discuss the primary constraint the measured data must meet
in order for the method to work, and provide a potential avenue to circumvent
the constraint in practice. In Subsection II-D, we briefly describe the tool
used to compute the Betti number sequence of a given data set, and how to use
this sequence to estimate the number of sources in the mixture.

###  II-A Topological Analysis of Monocomponent Mixtures

We begin by analyzing the mixing problem through the lens of topology. First,
let $x(t)=[x_{1}(t),x_{2}(t),...,x_{n}(t)]^{\top}$ be the vector of $n$
independent sources. If we consider each $x_{i}(t)$ as the motion of a point
along an orthogonal basis vector of $\mathbb{R}^{n}$, then $x(t)$ can be
interpreted as a parametric path through $\mathbb{R}^{n}$. We restrict each
$x_{i}(t)$ to be a continuous signal of the form
$x_{i}(t)=A_{i}\textrm{cos}(\alpha_{i}(t))$, where $A_{i}$ is the constant
amplitude of $x_{i}(t)$, and $\alpha_{i}(t)$ is a continuous function of time
encoding the instantaneous frequency and phase of $x_{i}(t)$. Since each
$\alpha_{i}$ could include a constant phase term in the interval $[-\pi,\pi]$,
we can say that these sinusoidal sources are “cosines” without any loss of
generality. Sources of this type are often referred to as constant-amplitude
“monocomponent signals” in the literature and are frequently encountered in
radar (e.g., chirps) and telecommunication (e.g., continuous phase frequency
modulated signals) applications.

Constant-amplitude monocomponent signals are of special interest since they
can be embedded as topological circles in $\mathbb{R}^{2}$ [9]. To see this,
consider that the Hilbert transform of $x_{i}(t)$ is immediately obtained as
$\widetilde{x}_{i}(t)=A_{i}\textrm{sin}(\alpha_{i}(t))$, where $A_{i}$ and
$\alpha_{i}$ are unchanged from their $x_{i}(t)$ counterparts [10] . Then, by
simple trigonometric identity, the expression
$x_{i}^{2}(t)+\widetilde{x}_{i}^{2}(t)$ equals the constant $A_{i}^{2}$ for
all times $t$, which is the definition of a circle in the plane [9].
Considering now each component $x_{i}(t)$ and its Hilbert transform
$\widetilde{x}_{i}(t)$ as the motion of a point along mutually orthogonal axes
in $\mathbb{R}^{2n}$, we find that the trajectory of this point (i.e., the
phase portrait of $x(t)$) forms a path on an $n$-torus since
$(\mathbb{S}^{1})^{n}=\mathbb{T}^{n}$ [11]. As discussed in [11], the path
itself may not actually become dense on the torus but form a torus knot
instead when the ratio of instantaneous frequencies of individual components
of $x(t)$ are rational throughout the observation window. This condition is
unlikely to occur for incoherent signals with nonstationary frequencies, so we
set aside this outlier case and assume that the path becomes dense on the
$\mathbb{T}^{n}$ manifold.

By the Künneth formula [12], the Betti number sequence for an $n$-torus is
given by the coefficients of the Poincaré polynomial $(1+q)^{n}$. This
sequence is a topological invariant of the $n$-torus, i.e., it is invariant
under homeomorphisms and embeddings [8]. The existence of this topological
invariant allows us to compute the number of sources in a mixture of
monocomponent signals as follows. Let $y(t)\in\mathbb{R}^{m}$ be the vector
composed of $m$ independent observations of mixtures of $x(t)$, and let
$z(t)\in\mathbb{R}^{2n}$ be defined as
$[x_{1}(t),\widetilde{x}_{1}(t),x_{2}(t),\widetilde{x}_{2}(t),...,x_{n}(t),\widetilde{x}_{n}(t)]^{\top}$.
Let $w(t)\in\mathbb{R}^{2m}$ be a vector defined by $w(t)=f(z(t))$, where $f$
is an arbitrary function mapping $\mathbb{R}^{2n}\to\mathbb{R}^{2m}$. When $f$
is a homeomorphism (when $m=n$) or an embedding (when $m>n$), the Betti number
sequences of the phase portraits of $w(t)$ and $z(t)$ will be equal. So, if we
can find a suitable method to embed the observed vector $y(t)$ into
$\mathbb{R}^{2m}$ as $w(t)$ such that $f$ exists and is a homeomorphism or
embedding, then we can simply compute the Betti number sequence of the phase
portrait of $w(t)$ to recover the Betti number sequence of the phase portrait
of $z(t)$. Subsection II-B provides such a method for the scenario of a
generic receiver array. Once computed, if we find that the Betti number
sequence of $z(t)$ matches the the coefficients of the polynomial $(1+q)^{n}$,
we know that there exist $n$ monocomponent sources in the mixture. In summary,
our estimation strategy is: (1) embed an $\mathbb{R}^{m}$ signal mixture into
$\mathbb{R}^{2m}$, (2) compute the Betti number sequence, and (3) compare the
sequence to the coefficients of $(1+q)^{n}$.

###  II-B Embedding Observations into $\mathbb{R}^{2m}$

In order to identify a suitable embedding of $y(t)$ into $\mathbb{R}^{2m}$ as
$w(t)$, we impose a constraint on the mixing process itself: that each
observation vector $y_{i}(t)$ is given by:

$y_{i}(t)=\sum_{j=1}^{n}B_{i,j}\,\textrm{cos}(\alpha_{j}(t)+\phi_{i,j})$ (1)

for some values $B_{i,j}>0$ and $\phi_{i,j}\in[-\pi,\pi]$. In other words,
each observation vector $y_{i}(t)$ is a sum of each of the sources of $x(t)$,
modified by a relative magnitude and phase at each measurement. This model was
chosen as a simplification of the case of a receiver array with multiple,
directional elements receiving mixtures of the source signals, as might be
encounted in radar, sonar, or communications applications. The directionality
of the receiver elements along with the physical spacing of the source signals
induces differences in magnitude among the components of each mixture.
Likewise, the physical spacing of the receiver’s antenna induces a small time
delay between the reception of different sources, which, if small, can be
approximated by a phase-shift [13].

An effective way to ensure that that the map
$f:\mathbb{R}^{2n}\to\mathbb{R}^{2m}$ ($m\geq n$) exists and is a
homeomorphism (or embedding) is to let $f$ be defined by a matrix $T$ such
that $w(t)=Tz(t)$. Then, $T$ is a linear map, which is guaranteed to be a
homeomorphism (or embedding) if $T$ has full column rank of $2n$. Letting
$B_{i,j}=R_{i,j}\cdot A_{i}$ for some $R_{i,j}$, and using the trigonometric
identity for angle sums, we rewrite $y_{i}(t)$ as:

$y_{i}(t)=\sum_{j=1}^{n}R_{i,j}\,\textrm{cos}(\phi_{i,j})\,x_{j}(t)-R_{i,j}\,\textrm{sin}(\phi_{i,j})\,\widetilde{x}_{j}(t)$
(2)

where $\widetilde{x}_{j}(t)$ is the Hilbert transform of $x_{j}(t)$. Based on
the form of Eq. 2, we can represent each component of $y(t)$ as a linear
combination of the components of $z(t)$: the coefficient terms
$R_{i,j}\,\textrm{cos}(\phi_{i,j})$ and $R_{i,j}\,\textrm{sin}(\phi_{i,j})$
form the entries of $T$. By taking the Hilbert transform of each $y_{i}(t)$,
we can immediately obtain a second set of $m$ observations given by:

$\widetilde{y}_{i}(t)=\sum_{j=1}^{n}R_{i,j}\,\textrm{sin}(\phi_{i,j})\,x_{j}(t)+R_{i,j}\,\textrm{cos}(\phi_{i,j})\,\widetilde{x}_{j}(t)$
(3)

where each $R_{i,j}$ and $\phi_{i,j}$ are unchanged between $y_{i}(t)$ and its
Hilbert transform $\widetilde{y}_{i}(t)$. This second set of $m$ observations
is likewise a linear combination of the components of $z(t)$, providing the
additional $m$ rows of $T$. Therefore, $w(t)$ is defined as
$[y_{i}(t),\widetilde{y}_{i}(t)]^{\top}$, and $w(t)=Tz(t)$. We note that in
practice the embedding observation vector $\widetilde{y}_{i}(t)$ can be
obtained either through direct measurement (e.g., analog phase shifter
circuitry), or through digital analysis via the discrete Hilbert transform.

A useful result of requiring that the map
$f:\mathbb{R}^{2n}\to\mathbb{R}^{2m}$ $(m\geq n)$ be a linear transformation
$T$ is that if $T$ exists and is full rank, then there also must exist a
pseudo-inverse of $T$ that recovers the source vectors from the observation
vectors. Therefore when the method is successful, in addition to the mere
assertion that there are $n$ sources in the mixture, we can further assert
that the $n$ sources are able to be unmixed by finding the appropriate
pseudoinverse, e.g., via a least squares fit of transformed data to the
$n$-dimensional square torus. This is a powerful insight that will be explored
in future work.

Another potential benefit of this method is that the only constraint on $T$ is
that $T$ be full rank, and so exact knowledge of the entries of $T$ (i.e.,
each $B_{i,j}$ and $\phi_{i,j}$) is unnecessary. In practical terms, this
means that prior knowledge of the array geometry, incidence angle, propagation
velocity, or wave frequencies are generally not required for the method to
work. This could allow the method to be used on mobile, distributed, and/or
dynamic arrays of receiver elements.

###  II-C Defining Independence in the Observations

As explained in Section II-B, $T$ must have full column rank under this
“Hilbert embedding” for our topological feature recovery strategy to work, so
we now discuss the conditions under which this is true. Label each row of $T$
corresponding to the equation $y_{i}(t)$ as $r_{i}(t)$ and each row of $T$
corresponding to the equation $\widetilde{y}_{i}(t)$ as
$\widetilde{r}_{i}(t)$. Then, each pair of rows
$[r_{i}(t),\widetilde{r}_{i}(t)]^{\top}$ forms $n$ blocks of size $2\times 2$
of the form:

$\begin{bmatrix}R_{i,j}\cos(\phi_{i,j})&-R_{i,j}\sin(\phi_{i,j})\\\
R_{i,j}\sin(\phi_{i,j})&R_{i,j}\cos(\phi_{i,j})\end{bmatrix}$ (4)

where the terms have the same meanings as in Eqs. 2 and 3. This is the
familiar form of a matrix in $\mathbb{C}^{m\times n}$ as implemented in
$\mathbb{R}^{2m\times 2n}$. As such, we construct a dual of the matrix $T$,
called $U$, where each complex term of $U$ replaces each $2\times 2$ block
from Eq. 4 with $U_{i,j}=R_{i,j}\angle{\phi_{i,j}}$ (in phasor notation).
Under this construction, we can say that $T$ has full column rank of $2n$ if
and only if $U$ has full column rank $n$. Labeling each row of $U$ as $u_{i}$
for $i\in[1,m]$, we note that each $u_{i}$ corresponds to the coefficients of
the complex observation vector
$v_{i}(t)=y_{i}(t)+\sqrt{-1}\,\widetilde{y}_{i}(t)$. Then, for $U$ to possess
rank $n$, there must exist a subset of $n$ rows chosen from $u_{i}$ that are
linearly independent. Since the effect of multiplying the complex vector
$v_{k}(t)$ by a complex constant $c_{k}=A_{k}\angle{\theta_{k}}$
($A_{k}\in\mathbb{R}$, $\theta_{k}\in[-\pi,\pi]$) is to multiply the
magnitudes of the associated $y_{i}(t)$ and $\widetilde{y}_{i}(t)$ vectors by
$A_{k}$ and add $\theta_{k}$ to the phases of $y_{i}(t)$ and
$\widetilde{y}_{i}(t)$, we can informally say that a set of observations
$y(t)$ are “independent” if no element $y_{i}(t)$ is a linear combination
(over $\mathbb{R}$) of any phase-shifted versions of the others. Consequently,
it is the relative ratios in the magnitudes and relative differences in the
phases of each component among the $y_{i}(t)$ observations that determine the
rank of $T$.

We suspect that the condition for a full rank $T$ will hold in most practical
instances where $m\gg n$ (e.g., phased array with many elements) as only $n$
of the $m$ observation vectors must be “independent” in the meaning given
above. However, in cases where there are not expected to be $n$ independent
observation vectors available (e.g., $m<n$), it may be possible to derive
additional independent observation vectors through the use of a Time Delay
Embedding (TDE), as discussed in [11]. A full discussion of the theory, usage
and drawbacks of TDEs as embedding functions is beyond the scope of this
paper; the interested reader is referred to references [11, 14, 9, 15]. With
respect to monocomponent sources, TDEs have the potential to induce relative
phase shifts among the source components in a given mixture [9], and are
therefore a potential source for additional independent measurement vectors
when there are fewer array elements than sources present in the mixture.

###  II-D Topological Computation via Persistent Homology

As discussed in Sections II-A and II-B, once we have obtained a properly
embedded analysis signal $w(t)$, our next step is to compute the Betti number
sequence corresponding to the topology of the phase portrait of $w(t)$. The
TDA tool we choose to perform this computation with is Persistent Homology
(PH) [16], which has been used to study the topology of datasets with many
real-world applications [17, 18, 19, 20]. We direct the unfamiliar reader to
reference [21] for an introduction to the theory and computation of PH. Many
PH computation packages have been developed in recent years [22]; in Section
III, we use the JavaPlex software to compute PH due to its simple integration
into MATLAB [23].

Treating the PH computation as a “black box” processing step, we can estimate
the Betti number sequence of the input manifold through simple analysis of the
output “barcode” plot [23]. For the purposes of this paper, we determine the
Betti number sequence by simply counting the number of features in the barcode
plot that persist for at least half of the PH computation interval in each
dimension. We use this sequence to compute the number of sources by finding an
integer value $n$ for which our sequence matches the coefficients of the
polynomial $(1+q)^{n}$, as discussed in Section II-A. This integer value, $n$,
is the output of the method, providing the estimate of the number of sources
in the mixture.

In practice, we need not actually search through all possible $n$ to find a
match. As the coefficients of $(1+q)^{n}$ follow from the binomial theorem,
and ${n\choose 1}=n$, we let our initial “guess” of $n$ equal the second Betti
number (corresponding to homology group $H_{1}$) determined from the barcode
plot. We then verify that the remainder of Betti numbers from the sequence
match the remaining coefficients of $(1+q)^{n}$. If they all match, we have
confidence in our computation of the topology of $w(t)$ and thus $z(t)$, so we
assert that there are $n$ sources to be unmixed. If they do not match, then
the Betti numbers do not conform to those of a torus $\mathbb{T}^{n}$, as
there are no other possible values of $n$ which could produce the obtained
Betti number sequence. We assume in the non-matching case that our computation
has failed and we cannot accurately estimate the number of sources; this can
occur, for example, when the mixing matrix $T$ is either singular or otherwise
badly conditioned.

##  III Example Computation with Three Monocomponent Sources

Figure 1: Short-time Fourier Transform (STFT) Spectrograms of clean source
signals, along with one example of a mixture of all three sources with AWGN at
20 dB SNR. As can be seen, all three signals are highly nonstationary causing
overlap in time and frequency at various points.

We now present a proof-of-concept demonstration of the method using three
synthetic digital signals in MATLAB. While unknown to the receiver, the
received signals are a mixture of three independent monocomponent sources with
nonstationary frequency characteristics as could be encountered in a radar
scenario. The first source is a “barrage jammer” waveform whose instantaneous
frequency continuously sweeps the Nyquist range, while the other two are
“chirp” signals that begin and end at disparate instantaneous frequencies.
Plots of the spectrograms (i.e., short-time Fourier transforms) of the three
source signals, as well as an example of a noisy mixture of the three sources,
are provided in Figure 1. These signals were chosen somewhat arbitrarily as an
example of wideband, nonstationary signals that overlap at irregular points in
the time and frequency domains and are therefore inseparable through classical
linear filtering. The specific parameters (e.g., frequency sweep ranges,
initial phase, etc.) for these chosen signals are irrelevant; our topological
method is general enough that virtually any form of monocomponent signals is
viable under this method. The authors have conducted many additional trials
using alternate monocomponent signal types with results similar to those
presented herein.

Assuming that the receiver array contains eight elements (eight being an
arbitrary number substantially larger than the expected three sources), we
generate eight random values of $R_{i,j}\in[0.75,1.25]$ and
$\phi_{i,j}\in[-\pi,\pi]$ and create the eight observation signals,
$y_{i}[n]$, according to a discretized version of Eq. 1. Each observation
vector contains 30000 samples, equivalent to a 1-MHz sampling frequency
observed over 30 milliseconds. Additive white Gaussian noise (AWGN) is added
separately to each observation signal with a signal-to-noise ratio (SNR)
randomly selected from 15 to 25 dB. A variable SNR on each measurement is used
to provide some evidence of robustness in the presence of nonspatially white
noise. As discussed in Section II-B, we compute the Hilbert transform of each
of our observed $y_{i}[n]$ to obtain $\widetilde{y}_{i}[n]$. As a final
preparatory step, we remove the first and last 10% of samples from each
$y_{i}[n]$ and $\widetilde{y}_{i}[n]$ to minimize the windowing effects of
computing the Hilbert transform of a finite-duration signal. These trimmed
signals are combined into the 16-dimensional input signal $w[n]$ via
concatenation.

We then provide the 16-dimensional signal $w[n]$ as the point-cloud input to
JavaPlex. We use the Witness Complex construction with 150 landmark points
chosen using the sequential min-max procedure, up to a maximum filtration
distance of 0.24 [23]. The results are provided in barcode format in Figure 2,
showing persistent Betti numbers of $\\{1,3,3,1,0,0,0,0,0\\}$. As up to eight
sources could be detected with this method, the PH was computed over
dimensions 0-8; however, there were no features found in dimensions 4-8, so
these empty results are omitted from Figure 2 to save space. It is trivial to
identify that the computed Betti number sequence of $\\{1,3,3,1\\}$ matches
the coefficients of the polynomial $(1+q)^{3}=1+3q+3q^{2}+q^{3}$, thus the
method computes that there are exactly 3 sources making up the mixed
observation vectors, as was desired.

For comparison purposes, we also computed the outputs of the MDL and AIC
estimators as defined in [6] on the same data set. Surprisingly, these
estimators improperly estimate that the mixtures consist of 7 independent
sources. Similar results are obtained when the parameters of the mixtures are
changed and different noise levels are generated. As discussed in [5], this
error in the MDL and AIC estimators is likely due to variations in AWGN power
among the observed signals. While far from conclusive, this difference in
outcomes hints that the topological method may provide additional robustness
over the statistical methods in certain scenarios, and thus the method
warrants additional research.

Figure 2: Barcodes corresponding to the PH computation in dimensions 0-3. The
Betti numbers of $\\{1,3,3,1\\}$ are found by counting the number of lines in
each dimension which persist over a majority of the filtration interval.

##  IV Conclusions and Future Work

In this paper, we introduced a topology-based method to estimate the number of
source signals present in a linear, monocomponent mixture. In particular, we
showed how the analytic form of a mixture of $n$ monocomponent signals
generates a phase space which is homeomorphic to the $n$-torus, so by using
TDA tools to recover the topological features of the mixture, we can determine
the number of signals making up the mixture. While this method is only
directly applicable to constant-amplitude monocomponent source signals in its
current form, these signals are nonetheless commonplace in wireless
communications and radar scenarios. Therefore, the method could have
substantial real-world applicability without significant modification.

Our provided demonstration of the method is merely a proof-of-concept, as
there are many aspects to the method that require further research in order to
refine. In particular, the TDA tool of PH is known to suffer from high
computational complexity and sensitivity to noisy outliers in the data set
[23]. These issues are being addressed by other researchers in the literature
[22], but the nuances are complex. Accordingly, we omitted such discussions in
this paper as they distract from the presentation of the underlying theory of
the technique. Since we consider the PH computation as a “black box” for
computing the Betti number sequence of the underlying mixture, any
improvements in the robustness, accuracy, or efficiency of the PH computation
can likewise be leveraged by our method. We suspect additional preprocessing
of the observation signals will also play a role in improving the method;
these techniques will be explored in future work.

As briefly demonstrated in Section III, the method can potentially outperform
the standard MDL and AIC methods in some specific scenarios, such as when the
received noise is not spatially white. While this outcome encourages further
research into the method, a more robust investigation of this method is needed
to adequately compare performance with the existing statistical methods. Such
a comparison requires additional control over the scenario model, any
preprocessing and optimization steps, and variations in the numbers of
sources, signal types, measurements, and noise types/levels. Due to the large
number of experimental parameters and alternate techniques to test, we have
opted to save such an investigation for future publication, and present only
the theory and proof-of-concept of our method here.

More generally, we believe the novelty of the approach is of sufficient
interest: eschewing the usual tools of statistical analysis to instead analyze
signals on the basis of shape. As topological data analysis is a rapidly
evolving field, we believe that additional practical algorithms and tools will
emerge to tackle current problems in many fields such as signal processing,
control theory, and communications. As such, we intend to not only improve on
this particular application, but also leverage the topological approach taken
to search for novel solutions to other types of problems. We strongly
encourage other researchers to join us in this endeavour.

## References

  * [1] Z. Luo, C. Li, and L. Zhu,  “A Comprehensive Survey on Blind Source Separation for Wireless Adaptive Processing: Principles, Perspectives, Challenges and New Research Directions,”  IEEE Access, vol. 6, pp. 66685–66708, 2018.
  * [2] R. Schmidt,  “Multiple emitter location and signal parameter estimation,”  IEEE Transactions on Antennas and Propagation, vol. 34, no. 3, pp. 276–280, Mar. 1986.
  * [3] R. Roy and T. Kailath,  “ESPRIT-estimation of signal parameters via rotational invariance techniques,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 7, pp. 984–995, July 1989.
  * [4] L. C. Zhao, P. R. Krishnaiah, and Z. D. Bai,  “On detection of the number of signals in presence of white noise,”  Journal of Multivariate Analysis, vol. 20, no. 1, pp. 1–25, Oct. 1986.
  * [5] E. Fishler and H. V. Poor,  “Estimation of the Number of Sources in Unbalanced Arrays via Information Theoretic Criteria,”  IEEE Transactions on Signal Processing, vol. 53, no. 9, pp. 3543–3553, Sept. 2005.
  * [6] T. Salman, A. Badawy, T. M. Elfouly, A. Mohamed, and T. Khattab,  “Estimating the Number of Sources: An Efficient Maximization Approach,”  in 2015 International Wireless Communications and Mobile Computing Conference (IWCMC), Aug. 2015, pp. 199–204.
  * [7] A. Weiss and A. Yeredor,  “Blind Determination of the Number of Sources Using Distance Correlation,”  IEEE Signal Processing Letters, vol. 26, no. 6, pp. 828–832, June 2019.
  * [8] “Betti number - Encyclopedia of Mathematics,” https://encyclopediaofmath.org/wiki/Betti_number.
  * [9] S. M. Kennedy, J. D. Roth, and J. W. Scrofani,  “A Novel Method for Topological Embedding of Time-Series Data,”  in 2018 26th European Signal Processing Conference (EUSIPCO), Rome, Sept. 2018, pp. 2350–2354, IEEE.
  * [10] E. Bedrosian,  “A Product Theorem for Hilbert Transforms,”  Tech. Rep., RAND Corporation, Dec. 1962.
  * [11] V. de Silva, P. Skraba, and M. Vejdemo-Johansson,  “Topological Analysis of Recurrent Systems,”  in Workshop of Algebraic Topology and Machine Learning, NIPS 2012, 2012, pp. 1–5.
  * [12] “Künneth formula - Encyclopedia of Mathematics,” https://encyclopediaofmath.org/wiki/K%C3%BCnneth_ formula.
  * [13] L. J. Ziomek,  An Introduction to Sonar Systems Enginering,  CRC Press, 2023.
  * [14] Floris Takens,  “Detecting strange attractors in turbulence,”  in Dynamical Systems and Turbulence, Warwick 1980, D. Rand and L. Young, Eds., vol. 898, pp. 366–381. Springer Berlin Heidelberg, Berlin, Heidelberg, 1981.
  * [15] M. Casdagli, S. Eubank, J. D. Farmer, and J. Gibson,  “A theory of state space reconstruction in the presence of noise,”  Tech. Rep. LA-UR-90-4342; CONF-9007195-1, Los Alamos National Lab., NM (USA), Jan. 1990.
  * [16] G. Carlsson,  “Topology and data,”  Bulletin of the American Mathematical Society, vol. 46, no. 2, pp. 255–308, Jan. 2009.
  * [17] F. Erden and A. E. Cetin,  “Period Estimation of an Almost Periodic Signal Using Persistent Homology With Application to Respiratory Rate Measurement,”  IEEE Signal Processing Letters, vol. 24, no. 7, pp. 958–962, July 2017.
  * [18] S. Emrani, T. Gentimis, and H. Krim,  “Persistent Homology of Delay Embeddings and its Application to Wheeze Detection,”  IEEE Signal Processing Letters, vol. 21, no. 4, pp. 459–463, Apr. 2014.
  * [19] J. A. Perea, A. Deckard, S. B. Haase, and J. Harer,  “SW1PerS: Sliding windows and 1-persistence scoring; discovering periodicity in gene expression time series data,”  BMC Bioinformatics, vol. 16, pp. 257, Aug. 2015.
  * [20] V. Selvaraju, P. Namadurai, and R. Swaminathan,  “Analysis of Electromyography Burst Signals using Topological Feature Extraction for Diagnosis of Preterm Birth,”  in 2020 IEEE Signal Processing in Medicine and Biology Symposium (SPMB), Dec. 2020, pp. 1–4.
  * [21] R. Ghrist,  “Barcodes: The persistent topology of data,”  Bulletin of the American Mathematical Society, vol. 45, no. 01, pp. 61–76, Oct. 2007.
  * [22] N. Otter, M. A Porter, U. Tillmann, P. Grindrod, and H. A Harrington,  “A roadmap for the computation of persistent homology,”  EPJ Data Science, vol. 6, no. 1, pp. 17, Dec. 2017.
  * [23] H. Adams, A. Tausz, and M. Vejdemo-Johansson,  “javaPlex: A Research Software Package for Persistent (Co)Homology,”  in Mathematical Software – ICMS 2014, Hoon Hong and Chee Yap, Eds., vol. 8592, pp. 129–136. Springer Berlin Heidelberg, Berlin, Heidelberg, 2014.


# Fusing Multiple Algorithms for Heterogeneous Online Learning

Darshan Gadginmath<EMAIL_ADDRESS>
Shivanshu Tripathi<EMAIL_ADDRESS>
Fabio Pasqualetti<EMAIL_ADDRESS>
University of California  Riverside This material is based upon work supported
in part by awards ONR-N00014-19-1-2264, ARO W911NF-20-2-0267, and
AFOSRFA9550-19-1-0235.

###### Abstract

This study addresses the challenge of online learning in contexts where agents
accumulate disparate data, face resource constraints, and use different local
algorithms. This paper introduces the Switched Online Learning Algorithm
(SOLA), designed to solve the heterogeneous online learning problem by
amalgamating updates from diverse agents through a dynamic switching mechanism
contingent upon their respective performance and available resources. We
theoretically analyze the design of the selecting mechanism to ensure that the
regret of SOLA is bounded. Our findings show that the number of changes in
selection needs to be bounded by a parameter dependent on the performance of
the different local algorithms. Additionally, two test cases are presented to
emphasize the effectiveness of SOLA, first on an online linear regression
problem and then on an online classification problem with the MNIST dataset.

##  1 Introduction

Multi-agent learning frequently involves scenarios in which agents gather
disparate data at varying rates, collectively seeking to address an online
optimization problem. Some instances include collaborative localization,
search-and-rescue operations, coverage control, etc. Compounding the
complexity of these scenarios is the constraint of limited data processing and
computational capabilities, and the need for time-sensitive decision-making.
Conventionally, either heterogeneous data or gradients are pooled [Zhang et
al.(2018)Zhang, Cisse, Dauphin, and Lopez-Paz, Esfandiari et
al.(2021)Esfandiari, Tan, Jiang, Balu, Herron, Hegde, and Sarkar, Aketi et
al.(2023)Aketi, Kodge, and Roy], or heterogeneity is ignored and a distributed
learning algorithm is deployed [Le Bars et al.(2023)Le Bars, Bellet, Tommasi,
Lavoie, and Kermarrec, Zhu et al.(2021)Zhu, Xu, Liu, and Jin]. The former
approach raises privacy concerns, while the latter proves suboptimal due to
the inherent heterogeneity of the data. Effectively addressing the
heterogeneity of data and resource constraints in an online learning problem
remains an unresolved challenge.

An alternative approach to mitigate the challenges posed by data heterogeneity
and resource constraints is to employ distinct algorithms tailored to the
specific characteristics of the data, computation, and communication
resources. Nevertheless, adopting distinct algorithms introduces the risk of
underutilizing the available data. In this study, we present a systematic
method to integrate updates provided by distinct algorithms to solve the
heterogeneous online learning problem.

###  1.1 Problem Statement

We seek to solve the following online minimization problem:

$\displaystyle\min_{x}F(x,\bigcup\limits_{i}^{M}\mathcal{D}^{i}(t)).$ (1)

Here, $x\in^{n}$ is a parameter that needs to be collectively estimated by
agents $i\in\\{1,2,\dots,M\\}$. The data gathered by each agent up to time $t$
is denoted by $\mathcal{D}^{i}(t)$. More precisely, let $T^{i}\in$ denote the
set of time instances when agent $i$ collects new data. Then
$T^{i}=\\{t^{i}_{1},t^{i}_{2},\dots\\}$, where $t^{i}_{j}$ denotes the
$j^{th}$ round of sampling by agent $i$. We use the set $T$ to denote all the
time instances when new data is available. That is,
$T=\bigcup_{i}~{}T^{i}=\\{t_{1},t_{2},\dots\\}$. We drop the superscript $i$
to denote a time instance when new data is acquired by any agent. The elements
of $T$ satisfy $t_{1}\leq t_{2}\leq\dots\bar{t}$, where
$\bar{t}=\max_{i,j}t^{i}_{j}$. Let the samples collected by agent $i$ at the
$j^{th}$ round be denoted by $s(t^{i}_{j})$. The data collected by agent $i$
until time $t$ is

$\displaystyle\mathcal{D}^{i}(t)=\bigcup\limits_{\begin{subarray}{c}t^{i}_{j}\in
T^{i},t^{i}_{j}\leq t\end{subarray}}s(t^{i}_{j}).$ (2)

Every agent $i$ is constrained by their computation power and communication
capabilities. Hence, they need to employ a local algorithm $\mathcal{A}^{i}$
with their local data to solve the online optimization problem. For instance,
a single agent can employ a centralized algorithm such as gradient descent
(GD), stochastic gradient descent (SGD), or batch-wise GD [Dixit et
al.(2019)Dixit, Bedi, Tripathi, and Rajawat]. An agent who is composed of a
system of smaller units can collectively employ distributed algorithms such as
Decentralized SGD or Federated Learning depending on the availability of
resources. As an example, consider the following situation: an e-commerce
entity equipped with data centers strategically situated across varied
geographical locations, adopts an asynchronous data acquisition methodology
from online users for the purpose of targeted advertisement display. In this
scenario, each center employs a localized algorithm to process its specific
dataset, owing to limitations in both computational and communicative
capabilities.

Agents may choose to delay their decision-making process in order to
accumulate an ample amount of data or computational resources necessary for
solving problem (1). However, this approach may be suboptimal, given that
decisions are frequently time-sensitive. A practical example is that of naval
vessels mapping the sea for adversarial entities. Vessels positioned at
varying distances from the shoreline collect sonar data, yet their
computational and communicative capabilities are constrained by their
respective locations. Specifically, vessels in closer proximity to the coast
benefit from superior computational resources, albeit with diminished data
quality, as elucidated by [Ferla and Jensen(2002)]. In this scenario, waiting
to gather sufficient resources can be extremely dangerous. Therefore, the
updates from different algorithms need to be fused in an online fashion as
soon as updates are available to solve the problem (1).

###  1.2 Related Work

Distributed online optimization has been extensively studied as detailed in
[McMahan(2017), Hoi et al.(2021)Hoi, Sahoo, Lu, and Zhao, Li et al.(2023b)Li,
Xie, and Li]. However, these studies do not consider cooperatively using
different algorithms in a constrained time-sensitive setting. Popular
algorithms such as decentralized SGD and Federated Learning in the presence of
asynchronous agents were studied in [Jiang et al.(2021)Jiang, Zhang, Gu, and
Zhu, Chen et al.(2020)Chen, Ning, Slawski, and Rangwala]. In [Jiang et
al.(2021)Jiang, Zhang, Gu, and Zhu], decentralized SGD is proposed with
asynchronous agents but it does not incorporate agents running different
algorithms and it also requires extensive communication between agents.
Asynchronous online Federated Learning [Chen et al.(2020)Chen, Ning, Slawski,
and Rangwala] requires the presence of a coordinator, and still does not fuse
different algorithms. Model fusion has received attention in supervised
learning [Li et al.(2023a)Li, Peng, Zhang, Ding, Hu, and Shen]. Model fusion
in online learning case has received significantly less interest due to its
complexities [Foster et al.(2015)Foster, Rakhlin, and Sridharan, Hoang et
al.(2019)Hoang, Hoang, Low, and How, Cutkosky(2019)]. These works typically
consider selecting models from several algorithms at every time step.
Particularly, [Cutkosky(2019)] shows that algorithms with bounded regrets can
be fused simply by averaging the parameters and still maintain bounded regret.
However, in our case, we seek to fuse the updates from different agents by
employing only a single agent at any given time. Therefore, it is unclear how
agents with different data and resources could cooperatively solve problem (1)
by running their own local algorithm.

###  1.3 Contributions

We provide an algorithm called the Switched Online Learning Algorithm (SOLA)
to solve Problem (1) by fusing the updates from agents running different
algorithms. We solve the considered problem by switching between the agents
and fusing their updates based on their performance. We provide a sufficient
condition to guarantee a bound on the regret of SOLA based on the rate at
which different algorithms are chosen. To this end, we model SOLA as a
switched dynamical system and ensure its contractivity. We numerically analyze
the performance of our algorithm for the online linear regression problem and
also the online classification problem with the MNIST dataset111Code
repository: https://github.com/Shivanshu-007/Heterogeneous-online-optimization
.

##  2 Switched Online Learning Algorithm (SOLA)

In this section, we describe the proposed Switched Online Learning Algorithm
(SOLA). The input to any local algorithm $\mathcal{A}^{i}$ at time $t$ is the
data $\mathcal{D}^{i}(t)\in^{m(t)\times p}$, and the parameter $x\in^{n}$.
Note that the dimension of $\mathcal{D}^{i}(t)$ is dependent on time as each
agent acquires new data over time. The number of samples is given by $m(t)$
and the number of features is $p$. Henceforth, we simply use the variable $k$
to denote the discrete instances $t_{k}\in T$, i.e. $k\in\\{1,2,\dots,|T|\\}$.
The update provided by algorithm $\mathcal{A}^{i}$ is given by the map
$\mathcal{A}^{i}(x,\mathcal{D}^{i}(k)):^{n}\times^{m(k)\times
p}\rightarrow^{n}$. At any instance $k$, the _selecting signal_
$\sigma(k):\\{1,2,...,|T|\\}\rightarrow\\{1,2,\dots,M\\}$ selects an agent $i$
if the agent has new data. Agent $\sigma(k)$ uses its local algorithm
$\mathcal{A}^{\sigma(k)}$ to update the parameter $x(k-1)$. The update
provided by $\mathcal{A}^{\sigma(k)}$ is used to solve the problem (1) as

$\displaystyle x(k)=\alpha(k)\
\mathcal{A}^{\sigma(k)}(x(k-1),\mathcal{D}^{\sigma(k)}(k))+(1-\alpha(k))\
x(k-1).$ (3)

We call $\alpha(k)\in[0,1]$ as the _fusing variable_. We introduce it to
smoothly incorporate new updates from the chosen local algorithms
$\mathcal{A}^{\sigma(k)}$. The fusing variable depends on the performance of
the algorithm of the local algorithms. Let us define a performance metric as
$P(x,\mathcal{D}):^{n}\times^{m\times p}\rightarrow_{\geq 0}$. Common
performance metrics for classification problems are precision, recall, true
positive rate, etc. [Hand(2012)]. In regression problems, performance is often
computed as the inverse of the norm of error, or the trace of the inverse of
error covariance [Kay(1993)]. We are particularly interested in metrics that
have higher values to signify better performance. Given the performance
metric, the fusing variable $\alpha(k)$ is defined as:

$\displaystyle\alpha(k)$
$\displaystyle=\frac{P(x(k^{+}),\mathcal{D}^{\sigma(k)}(k))}{P(x(k^{+}),\mathcal{D}^{\sigma(k)}(k))+P(x(k-1),\mathcal{D}^{\sigma(k-1)}(k-1))},$
(4)

where $x(k^{+})=\mathcal{A}^{\sigma(k)}(x(k-1),D^{\sigma(k)}(k))$. Note that
$\alpha(k)$ is using the performance of $\mathcal{A}^{\sigma(k)}$ to
incorporate its update. If the performance is poor on the local data
$\mathcal{D}^{\sigma(k)}(k)$, the fusing variable $\alpha(k)$ is closer to
$0$, whereas $1-\alpha(k)$ is closer to 1. This implies that the former
parameter $x(k-1)$ has a higher influence on the updated parameter $x(k)$.
Similarly, if the performance of $\mathcal{A}^{\sigma(k)}$ is superior to that
of $x(k)$, more weight is given to $x(k^{+})$. We describe the proposed
algorithm SOLA in detail in Algorithm 1.

######  Remark 2.1.

(Naive switching): The design of the fusing variable $\alpha(k)$ is crucial to
SOLA because the updates from $\mathcal{A}^{\sigma(k)}$ can be vastly
different from $x(k)$. For instance, by setting $\alpha(k)=1$, the update of
one agent acts as the input to the subsequent agents. This is a naive method
of switching between agents without any consideration for the performance of
each agent. Switching naively between different agents may cause abrupt and
large changes in the parameters, which may be undesirable. This is illustrated
in the online regression problem represented in Figure 1. SOLA chooses between
two agents running centralized GD and decentralized SGD with five sub-units.
Figure 1 compares the naive case when $\alpha(k)=1$, and the case when
$\alpha(k)$ is dependent on the performance of the local algorithms. We see
frequent jumps in the parameter $x(k)$ resulting in a chattering behavior with
an improper choice of $\alpha(k)$. However, when $0\leq\alpha(k)\leq 1$, we
see an improvement in the performance and the chattering behavior is absent.
The detailed simulation setting is provided in Section 4.1.

time ($t$)$\|x(t)-x^{*}\|$ Figure 1: Online linear regression with SOLA
comparing a naive choice of the fusing variable $\alpha_{t}=1$ and a fusing
variable that incorporates the performance of each agent. It is evident that a
naive choice of the fusing variable results in frequent jumps of the parameter
$x(t)$.

######  Remark 2.2.

(Choice of performance metric): SOLA uses updates from different agents for a
general online optimization problem. Therefore, the optimal choice of a
performance metric is dependent on the cost function, data and the constraints
of each agent. For instance, consider the online linear regression problem.
The recursive least squares algorithm in [Kay(1993)] incorporates data
arriving incrementally to solve an online linear regression problem. Although
the algorithm operates with only a single agent, it incorporates incoming
updates by weighting them with the inverse of the error covariance. It can be
shown that the choice of using the inverse of the error covariance in the
recursive least-squares algorithm produces the Best Linear Unbiased Estimator
(BLUE).

Algorithm 1 Switched Learning Algorithm

Initialize: $x(0),~{}\sigma(0),~{}\mathcal{D}^{\sigma(0)}(0)$.

1\. For $k=1$ to $|T|$:

2\. Select agent $i$ using switching signal, $\sigma(k)=i$.

3\. Agent $\sigma(k)$ receives $x(k-1)$ and
$P(x(k-1),\mathcal{D}^{\sigma}(k-1)(k-1))$ from $\sigma(k-1)$.

4\. Update parameter $x(k-1)$ using algorithm $\mathcal{A}^{\sigma(k)}$:
$x(k^{+})=\mathcal{A}^{\sigma(k)}(x(k-1),D_{\sigma(k)}(k))$.

5\. Compute performance metric $P(x(k^{+})),D_{\sigma(k)}(k))$ and fusing
variable $\alpha(k)$ (4).

6\. Incorporate update:
$x(k)=\alpha(k)\mathcal{A}^{\sigma(k)}(x(k-1),\mathcal{D}^{\sigma(k)}(k))+(1-\alpha(k))x(k-1)$.

##  3 Design of Selecting Signal

The selecting signal $\sigma(k)$ determines the choice of subsystem for the
update, and in turn, determines the regret of SOLA. In this section, we
analyze the effect of the switching signal on the regret of the Switched
Learning Algorithm. Let us denote any algorithm employed to solve problem (1)
as $\mathcal{A}$. The regret of algorithm $\mathcal{A}$ is denoted by
$\mathcal{R}_{\mathcal{A}}$ and is defined for $K\geq 1$ as

$\displaystyle\mathcal{R}_{\mathcal{A}}(K)$
$\displaystyle=\sum_{k=1}^{K}F(x(k),\mathcal{D}(k))-\sum_{k=1}^{T}F(x^{*},\mathcal{D}(k)).$
(5)

Here $x(k)$ refers to the parameters provided by algorithm $\mathcal{A}$ at
time $k$ using data $\mathcal{D}(k)$, and $x^{*}$ is the optimal parameter
given all the data a priori. It is essential for an algorithm to have bounded
regret as it compares the performance of the algorithm to an optimal choice.
We now show that SOLA has bounded regret under the right conditions of the
selecting signal $\sigma(k)$.

We analyze the regret achieved by SOLA by viewing the algorithm as a dynamical
system. Optimization algorithms have received extensive attention from the
view of dynamical systems [Elisseeff et al.(2005)Elisseeff, Evgeniou, Pontil,
and Kaelbing, Shikhman and Stein(2009), Ross and Bagnell(2011), Saha et
al.(2012)Saha, Jain, and Tewari, Sontag(2022), Cisneros-Velarde and
Bullo(2022), Kozachkov et al.(2022)Kozachkov, Wensing, and Slotine].
Particularly, algorithms that are stable have been shown to have bounded
regret [Ross and Bagnell(2011), Saha et al.(2012)Saha, Jain, and Tewari]. A
useful tool to study these optimization algorithms and stability is
contraction theory [Cisneros-Velarde and Bullo(2022), Kozachkov et
al.(2022)Kozachkov, Wensing, and Slotine]. A contracting algorithm is defined
as follows.

######  Definition 3.1.

(Contracting algorithm, [Lohmiller and Slotine(1998)]) Let the updates
provided by an algorithm $\mathcal{A}$ be given by the dynamical system

$\displaystyle x(k)=\mathcal{A}(x(k-1),\mathcal{D}(k))$ (6)

where $x(k)$ and $\mathcal{D}(k)$ are the parameter and data at time $k$,
respectively. The differential dynamics of the algorithm is then given by

$\displaystyle\delta_{x(k)}=\frac{\partial\mathcal{A}(x(k-1),\mathcal{D}(k))}{\partial
x(k-1)}\delta_{x(k-1)}.$

The associated distance for the differential dynamics is denoted by

$\displaystyle V(x(k))$
$\displaystyle=\delta_{x(k)}^{\top}M(x(k))\delta_{x(k)},$ $\displaystyle
M(x(k))$
$\displaystyle=\left(\frac{\partial\mathcal{A}(x(k),\mathcal{D}(t))}{\partial
x(k)}\right)^{\top}\frac{\partial\mathcal{A}(x(k),\mathcal{D}(t))}{\partial
x(k)}.$

Here $M(x(k))$ is a symmetric positive-definite matrix function, and it is
uniformly bounded as $\gamma_{1}I\leq M(x(k))\leq\gamma_{2}I$, for all $x(k)$.
The algorithm $\mathcal{A}$ is said to be contracting if

$\displaystyle V(x(k))\leq\beta V(x(k-1)),$ (7)

where $0<\beta\leq 1$, is the rate of contraction.

We show that any contracting algorithm achieves bounded regret in Appendix A.
Given that contracting algorithms achieve bounded regret, it is sufficient to
ensure that SOLA is a contracting algorithm by designing the switching signal
$\sigma(k)$. We make the following assumptions for our regret analysis.

  1. (A1)

$F(x,\mathcal{D})$ is $\ell$-convex in $x$: $\nabla^{2}F(x,\mathcal{D})>\ell
I$.

  2. (A2)

Every local algorithm $\mathcal{A}^{i}$ is contracting with a rate
$\beta^{i}$. Further, for each pair of algorithms $(i,j)$, there exists
$\mu^{ij}>1$ such that the distance of differential dynamics is bounded as
$V^{i}(x(k))=\mu^{ij}V^{j}(x(k))$.

Assumption (A1) is a commonly used in regret analysis for online problems
[McMahan(2017)]. However, (A2) needs more careful attention as it comes from
the perspective of dynamical systems. Several algorithms such as gradient
descent [Cisneros-Velarde and Bullo(2022)], SGD and decentralized SGD
algorithms have been proven to be contractive [Boffi and Slotine(2020)]. When
these algorithms are used to solve a common online learning problem,
assumption (A2) captures the relationship between the differential dynamics of
the local algorithms. Particularly, it provides an upper bound on the relative
difference between the performance of the various local algorithms. We
introduce $\bar{\mu}$ and $\bar{\beta}$ which will be useful for our regret
analysis.

$\displaystyle\bar{\mu}=\mathrm{arg}\max\limits_{ij}\mu^{ij},\quad\bar{\beta}=\mathrm{arg}\max\limits_{i}\beta^{i}$

The variables $\bar{\mu}$ corresponds to the largest difference in
differential dynamics and $\bar{\beta}$ corresponds to the slowest contracting
rate. The following theorem addresses the design of the selecting signal such
that SOLA achieves bounded regret.

######  Theorem 3.2.

(SOLA achieves bounded regret) For a selecting signal $\sigma(k)$, let the
number of switches between different agents be $N(k_{1},k_{2})$ over any
horizon $\left[k_{1},k_{2}\right]$, where $k_{2}>k_{1}$. That is,

$\displaystyle N(k_{1},k_{2})$
$\displaystyle=\sum\limits_{k=k_{1}+1}^{k_{2}}r(k),\quad\mathrm{where~{}}r(k)=\begin{cases}0,\quad\sigma(k)=\sigma(k-1),\\\
1,\quad\mathrm{otherwise}.\end{cases}$ (8)

Then, if the number of switches satisfies

$\displaystyle N(k_{1},k_{2})$ $\displaystyle\leq
N_{0}+\frac{k_{2}-k_{1}}{\tau},$ (9)

where $N_{0}>0$ is a constant, and
$\tau=-\frac{\ln(\bar{\mu})}{\ln(\bar{\beta})}$, then SOLA achieves bounded
regret: $\mathcal{R}(K)\leq\epsilon(K)$, where $\epsilon(K)$ is a decreasing
function in $K$.

We refer the reader to Appendix B for the proof of Theorem 3.2. In the
literature related to switched systems [Liberzon(2003)], $N_{0}$ is typically
referred to as the chatter bound and $\tau$ is referred to as the average
dwell time.

######  Remark 3.3.

(Effect of cost function, data, and algorithm on the number of switches) It is
important to note that $\tau$ captures the effect of the cost function, the
data and the characteristics of the different local algorithms on the
admissible switching signal in Theorem 3.2. For instance, poor data or a slow
learning rate for the agent $i$ results in high values of $\beta^{i}$. This
leads to a higher value of $\bar{\beta}$, which in turn means that the number
of switches $N(k_{1},k_{2})$ needs to be small. Conversely, algorithms that
learn fast allow for fast switching. Another aspect of $\tau$ is the
similarity of the behavior of local algorithms captured through Assumption
(A2). If the different local algorithms behave similarly, $\bar{\mu}\approx
1$. This leads to a small value of $\tau$. Conversely, if the algorithms
behave very differently $\bar{\mu}$ is larger which restricts the frequency of
switching local algorithms.

##  4 Numerical results

In this section, we show the effectiveness of SOLA for online linear
regression and online classification using the MNIST dataset.

###  4.1 Online linear regression

We conduct experiments for online linear regression with a synthetic dataset.
The agents acquire data $\big{(}A^{i}(k),B^{i}(k)\big{)}$, where
$A^{i}(k)\in^{1\times m(k)}$, $B^{i}(k)\in^{3\times m(k)}$, and $x\in^{3\times
1}$. The online linear regression problem is

$\displaystyle\min\limits_{x}\left|\left|A^{\sigma(k)}(k)-x^{\top}~{}B^{\sigma(k)}(k)\right|\right|^{2}_{2}.$

In this experiment, the number of agents $M=2$. One agent acquires data
sampled as $A^{1}(k)={x^{*}}^{\top}B^{1}(k)+\zeta^{1}(k)$, where
$\zeta^{1}(k)\sim\mathcal{N}(0,3I)$. The other agent collects data
$A^{2}(k)={x^{*}}^{\top}B^{2}(k)+\zeta^{2}(k)$, where
$\zeta^{2}(k)\sim\mathcal{N}(0,30I)$. Here,
$B^{1}(k),B^{2}(k)\sim\mathcal{N}(0,0.5I).$ Figure 2 compares the performance
of SOLA for different settings. The selecting signal $\sigma(k)$ periodically
chooses between the two agents every ten instances, i.e. $N(k,k+10)=1$. In
Figure 2, we compare the performance of SOLA in two different settings, (i)
Agent 1 has one sub-unit that uses centralized gradient descent, and agent 2
uses decentralized SGD with 5 sub-units, (ii) Agent 1 has 5 sub-units that
perform FedAvg [McMahan et al.(2017)McMahan, Moore, Ramage, Hampson, and y
Arcas] and agent 2 uses decentralized SGD. Lastly, we compare the performance
of only agent 2 using decentralized SGD without SOLA. From Figure 2, it is
evident SOLA with GD and decentralized SGD performs the best, whereas SOLA
with FedAvg with decentralized SGD converges slower. Decentralized SGD by
itself converges slower and has a higher error.

In Figure 2(b), we compare SOLA with $M=3$. The three agents perform
centralized GD, decentralized GD with five sub-units and FedAvg with five sub-
units. The data used by the agents is the same as mentioned above. It is
evident that a naive choice of the fusing variable not only causes more error
but also leads to chattering. SOLA has a higher error when $M=3$ as compared
to when $M=3$. This is because fusing more agents requires a good choice of
the fusing variable and data of good quality.

time(a)$\|x(t)-x^{*}\|$

time(b)

Figure 2: Online linear regression with SOLA. (a) The blue curve represents
SOLA choosing between agents using centralized GD and decentralized SGD, the
green curve represents SOLA choosing between agents performing FedAvg and
Decentralized SGD, whereas the red curve represents only decentralized SGD.
SOLA in both cases performs better than pure decentralized SGD. (b) SOLA with
3 agents performing gradient descent, decentralized SGD and FedAvg compared
with SOLA with a naive choice of the fusing variable, $\alpha=1$. SOLA with
good choice of $\alpha$ not only eliminate the chattering behavior, it also
has a better overall error.

###  4.2 Online classification

time(a)Accuracy

time(b)Loss

Figure 3: The testing accuracy and loss for SOLA on the MNIST data set. Here
the selecting signal chooses between an agent running decentralized SGD and
centralized GD, respectively using the MNIST dataset. (a) accuracy of SOLA,
(b) entropy loss for SOLA.

In this experiment, we consider the online classification problem with $M=2$
and the cost is the cross-entropy loss. Here, one agent samples data from the
MNIST dataset and employs centralized gradient descent. It uses a neural
network that has one hidden layer with 128 neurons. The second agent performs
decentralized SGD with sub-units having the same neural network architecture
as agent 1. We compare the performance of SOLA for different numbers of sub-
units for the second agent. In the first case, there are five sub-units and
each sub-unit acquires only two labels. For example, the first sub-unit
receives labels $0$ and $1$, the second sub-unit receives labels $2$ and $3$,
and so on. In the second case, there are ten sub-units and each sub-unit
receives only one distinct label. Further, every sub-unit is restricted to
possess only $128$ images (reducing computation). The images sampled by the
sub-units are from the MNIST dataset but marred with Gaussian noise of
variance $0.5$ and zero mean. Figures 3(a) and (b) show the accuracy and
entropy loss of SOLA on testing data over time. We denote the number of sub-
units as $N$ in Figure 3. We observed that the performance metric
$P(x,\mathcal{D})=\frac{1}{F(x,\mathcal{D})}$ gave the best performance for
both accuracy and loss. The selecting signal periodically chooses between
agents every five instances, i.e. $N(k,k+5)=1$. It can be seen that SOLA still
achieves an accuracy close to eighty-two percent for both cases of five and
ten sub-units. However, the naive fusing choice with five sub-units has lower
accuracy and more chattering.

##  5 Conclusions and Future Work

In this work, we considered the scenario where agents with different data and
resources use different local algorithms to solve an online learning problem.
Our proposed algorithm, SOLA, provides a way to systematically fuse the
updates from different algorithms and ensure that regret is bounded. We also
numerically analyze the performance of SOLA for different online learning
scenarios. Future directions include the case of dynamically changing data
distributions, tighter regret bounds and the adversarial case where byzantine
agents provide malicious updates.

##  Appendix A Contracting Optimizers Achieve Bounded Regret

The connection between the stability of learning algorithms and bounded regret
for online learning problems has been studied in [Poggio et al.(2011)Poggio,
Voinea, and Rosasco, Ross and Bagnell(2011), Saha et al.(2012)Saha, Jain, and
Tewari]. Particularly, in [Poggio et al.(2011)Poggio, Voinea, and Rosasco],
the notion of online stability is defined as follows:

######  Definition A.1.

(Online stability) An algorithm $\mathcal{A}$ is said to be online stable if

$\displaystyle\mathbb{E}\left|\left|F(x^{\mathcal{A}}(k),\mathcal{D}(k))-F(x^{\mathcal{A}}(k-1),\mathcal{D}(k-1))\right|\right|\leq\epsilon_{os}(k),\quad\forall
t,$ (10)

where $\epsilon_{os}(k)\rightarrow 0$ as $k\rightarrow\infty$.

The notion of online stability captures the fact that for an online stable
algorithm, the change in cost incurred between any time instances is bounded
by a non-increasing function $\epsilon_{os}(k)$. Importantly, Theorem 18 in
[Ross and Bagnell(2011)] shows that online-stable algorithms achieve bounded
regret where $\epsilon(K)\leq\sum_{k=1}^{K}\epsilon_{os}(k)$. Also, [Poggio et
al.(2011)Poggio, Voinea, and Rosasco] shows that iterative gradient-based
methods such as GD, and SGD achieve online stability. We first model any
iterative gradient descent algorithm $\mathcal{A}$ as a perturbated gradient
descent.

$\displaystyle x(t)$
$\displaystyle=\mathcal{A}(x(t-1),\mathcal{D})=x(t-1)-\eta~{}(\nabla
F(x(t-1),\mathcal{D})+\xi(k)).$ (11)

The perturbation to the true gradient is
$\xi(k)\sim\mathcal{N}(0,\Sigma(\nabla F(x(k-1)))$ and $\eta$ is a constant
learning rate. Further, the covariance of the perturbation $\Sigma(\nabla
F(x(k)))$ satisfies $\lim\limits_{t\rightarrow\infty}\Sigma(\nabla
F(x(k)))=\bm{0}$, where $\bm{0}$ denotes the zero matrix. This model is
commonly used in the analysis of algorithms such as SGD [Jastrzebski et
al.(2017)Jastrzebski, Kenton, Arpit, Ballas, Fischer, Bengio, and Storkey,
Mandt et al.(2017)Mandt, Hoffman, and Blei, Li and Orabona(2019)] and
decentralized SGD [Boffi and Slotine(2020)]. The dynamics of any single sub-
unit of decentralized SGD or FedAvg can be expressed with (11). Further, the
average parameter of all the sub-units also follows the same dynamics,
however, the noise characteristics of $\xi$ differ. Here, we show that
perturbed iterative gradient-based algorithms are contracting as given by
Definition 3.1, and are are online stable.

######  Theorem A.2.

(Gradient-based algorithms are contractive and achieve online stablility)
Consider any iterative gradient-based algorithm $\mathcal{A}$ which updates
the parameter $x$ as given by (11). If an algorithm $\mathcal{A}$ is
contracting by Definition 3.1 in expectation, i.e.,

$\displaystyle\mathbb{E}_{\xi}\left[V(x(k))\right]\leq\beta~{}\mathbb{E}_{\xi}\left[V(x(k-1))\right],$

then $\mathcal{A}$ achieves online stability.

######  Proof A.3.

The differential dynamics of the system (11) is given by

$\displaystyle\delta_{x(k)}$
$\displaystyle=P(x(k-1))\delta_{x(k-1)}-\eta\delta_{\xi(k-1)},$ (12)
$\displaystyle P(x(k))$ $\displaystyle=\frac{\partial\big{(}x(k)-\eta~{}\nabla
F(x(k))\big{)}}{\partial x(k)}=I-\eta\nabla^{2}F(x(k))$ (13)

Here $\xi$ is treated as an external signal. For the Euclidean metric
$M(x(k))=I$, the distance of the differential dynamics is given by

$\displaystyle
V(x(k))=\delta_{x(k-1)}P(x(k-1))^{2}\delta_{x(k-1)}-2\eta\delta_{\xi(k-1)}^{\top}P(x(k-1))\delta_{x(k-1)}+\eta^{2}||\delta_{\xi(k-1)}||^{2}.$
(14)

Note that $P(x(k))<I-\eta\ell I$ due to the $\ell$-convexity of $F$.
Therefore,

$\displaystyle
V(x(k))<(1-\eta\ell)^{2}V(x(k-1))-2\eta\delta_{\xi(k-1)}^{\top}P(x(k-1))\delta_{x(k-1)}+\eta^{2}||\delta_{\xi(k-1)}||^{2}.$
(15)

Taking the expectation with respect to the perturbation $\xi(t)$,

$\displaystyle\mathbb{E}[V(x(k))]$
$\displaystyle<(1-\eta\ell)^{2}\mathbb{E}[V(x(k-1))]+\eta^{2}\mathrm{Var}[\xi(k-1)].$
(16)

If the learning rate $\eta$ ensures that

$\displaystyle(1-\eta\ell)^{2}\mathbb{E}[V(x(k))]+\eta^{2}\mathrm{Var}[\xi(k-1)]\leq\mathbb{E}[V(x(k))],$
(17)

then $\mathcal{A}$ is contractive. The proof of online stability follows along
the same lines as Theorem 2 of [Poggio et al.(2011)Poggio, Voinea, and
Rosasco]. The cost $F(x(t+1))$ can be expressed using the first two
derivatives by Taylor expansion. The stability of the perturbed gradient can
be then used to bound the change in cost as given by Definition A.1.

##  Appendix B Proof of Theorem 3.2

First, consider SOLA with the fusing variable $\alpha(k)=1$. Then,

$\displaystyle x(k)$
$\displaystyle=\mathcal{A}^{\sigma(k-1)}(x(k-1),\mathcal{D}(k-1)).$ (18)

By assumption (A2), we have that between any two instances $k_{1}$ and $k_{2}$

$\displaystyle V^{\sigma(k_{2})}(x(k_{2}))$
$\displaystyle\leq\bar{\beta}V^{\sigma(k_{2})}(x(k_{2}-1))\leq\bar{\mu}\bar{\beta}V^{\sigma(k_{2}-1)}(x(k_{2}-1))$
$\displaystyle\leq{\bar{\mu}}^{N(k_{1},k_{2})}{\bar{\beta}}^{(k_{2}-k_{1})}V^{\sigma(k_{1})}x(k_{1})=e^{N(k_{1},k_{2})\ln\bar{\mu}}~{}e^{(k_{2}-k_{1})\ln\bar{\beta}}~{}V^{\sigma(k_{1})}x(k_{1})$

To ensure that SOLA is a contraction, we need that

$\displaystyle N(k_{1},k_{2})\ln u^{*}+(k_{2}-k_{1})\ln\bar{\beta}$
$\displaystyle\leq 0$

To admit at least one switch if $\frac{(k_{2}-k_{1})}{\tau}<1$, we introduce
$N_{0}>0$. Therefore, when the number of switches $N(k_{1},k_{2})$ satisfies

$\displaystyle N(k_{1},k_{2})$ $\displaystyle\leq
N_{0}+\frac{(k_{2}-k_{1})}{\tau},$ (19)

SOLA is a contracting optimizer that achieves online stability. Further, when
$0\leq\alpha(k)\leq 1$, SOLA uses a convex combination of the update by the
local algorithm $\mathcal{A}^{\sigma(k)}$ and the parameter at the previous
time instance $x(k-1)$. The convex combination of contracting systems also
results in a contracting system as shown in [Lohmiller and Slotine(1998)].
Hence, the overall SOLA algorithm achieves online stability.

## References

  * [Aketi et al.(2023)Aketi, Kodge, and Roy] S. A. Aketi, S. Kodge, and K. Roy.  Neighborhood gradient mean: An efficient decentralized learning method for non-iid data.  _Transactions on Machine Learning Research_ , 2023.
  * [Boffi and Slotine(2020)] N. M. Boffi and J.-J. E. Slotine.  A continuous-time analysis of distributed stochastic gradient.  _Neural computation_ , 32(1):36–96, 2020.
  * [Chen et al.(2020)Chen, Ning, Slawski, and Rangwala] Y. Chen, Y. Ning, M. Slawski, and H. Rangwala.  Asynchronous online federated learning for edge devices with non-iid data.  In _2020 IEEE International Conference on Big Data (Big Data)_ , pages 15–24. IEEE, 2020.
  * [Cisneros-Velarde and Bullo(2022)] P. Cisneros-Velarde and F. Bullo.  A contraction theory approach to optimization algorithms from acceleration flows.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1321–1335. PMLR, 2022.
  * [Cutkosky(2019)] A. Cutkosky.  Combining online learning guarantees.  In _Conference on Learning Theory_ , pages 895–913. PMLR, 2019.
  * [Dixit et al.(2019)Dixit, Bedi, Tripathi, and Rajawat] Rishabh Dixit, Amrit Singh Bedi, Ruchi Tripathi, and Ketan Rajawat.  Online learning with inexact proximal online gradient descent algorithms.  _IEEE Transactions on Signal Processing_ , 67(5):1338–1352, 2019.
  * [Elisseeff et al.(2005)Elisseeff, Evgeniou, Pontil, and Kaelbing] A. Elisseeff, T. Evgeniou, M. Pontil, and L. P. Kaelbing.  Stability of randomized learning algorithms.  _Journal of Machine Learning Research_ , 6(1), 2005.
  * [Esfandiari et al.(2021)Esfandiari, Tan, Jiang, Balu, Herron, Hegde, and Sarkar] Y. Esfandiari, S. Y. Tan, Z. Jiang, A. Balu, E. Herron, C. Hegde, and S. Sarkar.  Cross-gradient aggregation for decentralized learning from non-iid data.  In _International Conference on Machine Learning_ , pages 3036–3046. PMLR, 2021.
  * [Ferla and Jensen(2002)] C. M. Ferla and F. B. Jensen.  Are current environmental databases adequate for sonar predictions in shallow water?  In _Impact of Littoral Environmental Variability of Acoustic Predictions and Sonar Performance_ , pages 555–562. Springer, 2002.
  * [Foster et al.(2015)Foster, Rakhlin, and Sridharan] D. J. Foster, A. Rakhlin, and K. Sridharan.  Adaptive online learning.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * [Hand(2012)] D. J. Hand.  Assessing the performance of classification methods.  _International Statistical Review_ , 80(3):400–414, 2012.
  * [Hoang et al.(2019)Hoang, Hoang, Low, and How] T. N. Hoang, Q. M. Hoang, K. H. Low, and J. How.  Collective online learning of gaussian processes in massive multi-agent systems.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 7850–7857, 2019.
  * [Hoi et al.(2021)Hoi, Sahoo, Lu, and Zhao] S. C. H. Hoi, D. Sahoo, J. Lu, and P. Zhao.  Online learning: A comprehensive survey.  _Neurocomputing_ , 459:249–289, 2021.
  * [Jastrzebski et al.(2017)Jastrzebski, Kenton, Arpit, Ballas, Fischer, Bengio, and Storkey] S. Jastrzebski, Z. Kenton, D. Arpit, N. Ballas, A. Fischer, Y. Bengio, and A. Storkey.  Three factors influencing minima in sgd.  _arXiv preprint arXiv:1711.04623_ , 2017.
  * [Jiang et al.(2021)Jiang, Zhang, Gu, and Zhu] J. Jiang, W. Zhang, J. Gu, and W. Zhu.  Asynchronous decentralized online learning.  _Advances in Neural Information Processing Systems_ , 34:20185–20196, 2021.
  * [Kay(1993)] S. M. Kay.  _Fundamentals of statistical signal processing: estimation theory_.  Prentice-Hall, Inc., 1993.
  * [Kozachkov et al.(2022)Kozachkov, Wensing, and Slotine] L. Kozachkov, P. Wensing, and J. J. E. Slotine.  Generalization as dynamical robustness–the role of riemannian contraction in supervised learning.  _Transactions on Machine Learning Research_ , 2022.
  * [Le Bars et al.(2023)Le Bars, Bellet, Tommasi, Lavoie, and Kermarrec] B. Le Bars, A. Bellet, M. Tommasi, E. Lavoie, and A.M. Kermarrec.  Refined convergence and topology learning for decentralized sgd with heterogeneous data.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1672–1702. PMLR, 2023.
  * [Li et al.(2023a)Li, Peng, Zhang, Ding, Hu, and Shen] W. Li, Yong Peng, Miao Zhang, Liang Ding, Han Hu, and Li Shen.  Deep model fusion: A survey.  _arXiv preprint arXiv:2309.15698_ , 2023a.
  * [Li and Orabona(2019)] X. Li and F. Orabona.  On the convergence of stochastic gradient descent with adaptive stepsizes.  In _The 22nd international conference on artificial intelligence and statistics_ , pages 983–992. PMLR, 2019.
  * [Li et al.(2023b)Li, Xie, and Li] X. Li, L. Xie, and N. Li.  A survey on distributed online optimization and online games.  _Annual Reviews in Control_ , 56:100904, 2023b.
  * [Liberzon(2003)] D. Liberzon.  _Switching in systems and control_ , volume 190.  Springer, 2003.
  * [Lohmiller and Slotine(1998)] W. Lohmiller and J.-J. E. Slotine.  On contraction analysis for non-linear systems.  _Automatica_ , 34(6):683–696, 1998.
  * [Mandt et al.(2017)Mandt, Hoffman, and Blei] S. Mandt, M. D. Hoffman, and D. M. Blei.  Stochastic gradient descent as approximate bayesian inference.  _arXiv preprint arXiv:1704.04289_ , 2017.
  * [McMahan et al.(2017)McMahan, Moore, Ramage, Hampson, and y Arcas] B. McMahan, Eider Moore, Daniel Ramage, Seth Hampson, and Blaise Aguera y Arcas.  Communication-efficient learning of deep networks from decentralized data.  In _Artificial intelligence and statistics_ , pages 1273–1282. PMLR, 2017.
  * [McMahan(2017)] H. B. McMahan.  A survey of algorithms and analysis for adaptive online learning.  _The Journal of Machine Learning Research_ , 18(1):3117–3166, 2017.
  * [Poggio et al.(2011)Poggio, Voinea, and Rosasco] T. Poggio, S. Voinea, and L. Rosasco.  Online learning, stability, and stochastic gradient descent.  _arXiv preprint arXiv:1105.4701_ , 2011.
  * [Ross and Bagnell(2011)] S. Ross and J. A. Bagnell.  Stability conditions for online learnability.  _arXiv preprint arXiv:1108.3154_ , 2011.
  * [Saha et al.(2012)Saha, Jain, and Tewari] A. Saha, P. Jain, and A. Tewari.  The interplay between stability and regret in online learning.  _arXiv preprint arXiv:1211.6158_ , 2012.
  * [Shikhman and Stein(2009)] V. Shikhman and O. Stein.  Constrained optimization: projected gradient flows.  _Journal of optimization theory and applications_ , 140:117–130, 2009.
  * [Sontag(2022)] E. D. Sontag.  Remarks on input to state stability of perturbed gradient flows, motivated by model-free feedback control learning.  _Systems & Control Letters_, 161:105138, 2022.
  * [Zhang et al.(2018)Zhang, Cisse, Dauphin, and Lopez-Paz] H. Zhang, M. Cisse, Y. N. Dauphin, and D. Lopez-Paz.  mixup: Beyond empirical risk minimization.  In _International Conference on Learning Representations_ , 2018.
  * [Zhu et al.(2021)Zhu, Xu, Liu, and Jin] H. Zhu, J. Xu, S. Liu, and Y. Jin.  Federated learning on non-iid data: A survey.  _Neurocomputing_ , 465:371–390, 2021.