# Computational fluid dynamics approach
for understanding oscillating and interacting convective flows

Attila Gergely and Zoltán Néda<EMAIL_ADDRESS>Physics Department,
Babeş-Bolyai University,
Kogalniceanu street 1, Cluj-Napoca, Romania

###### Abstract

A 2D numerical hydrodynamics approach is considered for modelling recent
experimental results on the oscillation and collective behavior of convective
flows. Our simulations consider the rising dynamics of heated fluid columns in
a gravitational field. Simulations are done on two entirely different length-
scales, showing also the generality of the investigated phenomena. For the
flow of a single heated fluid column, the effect of the inflow yield and the
nozzle diameter is studied. In agreement with the experiments, for a constant
nozzle diameter the oscillation frequency increases approximately linearly as
a function of the the flow yield and for a constant flow yield the frequency
decreases as a power law with the increasing nozzle diameter. Concerning the
collective behavior of two nearby flow we find a counter-phase synchronization
of the oscillations and an increasing trend in the common oscillation
frequency when the distance between the flows is decreased. These results are
again in agreement with the experimental findings.

numerical fluid dynamics, convective flows, oscillations, synchronization

##  I Introduction

Recent experimental results have proven that rising gas columns can perform
oscillations and their interaction leads to fascinating collective behavior
Gergely  _et al._ (2021). These oscillations have been reported much earlier
by Yuan et. al Yuan, Durox, and Villermaux (1994) and can be discussed in
analogy with the very similar phenomena known for diffusive flames Chamberlin
and Rose (1948); Durox  _et al._ (1995); Durox, Yuan, and Villermaux (1997);
Huang  _et al._ (1999); Kitahata  _et al._ (2009); Ghosh  _et al._ (2010);
Okamoto  _et al._ (2016); Chen  _et al._ (2019); Gergely  _et al._ (2020). The
toy-model presented in Gergely  _et al._ (2021) is clearly inadequate to
explain the fine details observed in our recent experiments, therefore a more
sophisticated theoretical approach is needed. The present study contributes in
this sense, considering a 2D numerical hydrodynamics approach to this
phenomenon.

For validating the numerical results, the present work considers as reference
the experiments realized using a controlled flow of Helium into air Gergely
_et al._ (2021) . In these experiments the Schlieren technique Settles (2001);
Leptuch and Agrawal (2006) was used to visualize the flow, allowing also a
digital processing of the oscillations. From the images processed by the Otsu
method Otsu (1979); FreddyCree (2020) the characteristic frequency and the
relevant synchronization order parameter was derived. For a better
understanding of the phenomena, some sample movies with original recordings
and the ones processed with the Otsu method are given in our YouTube channel
Gergely (2020a). For a single flow column the experiments investigated the
effect of the nozzle diameter and flow rate on the observed oscillation
frequency. For the collective behavior of two nearby flows the experiments
described in Gergely  _et al._ (2021) investigated the phase difference
between the oscillations, their frequency and a proper synchronization order
parameter as a function of the separation distance between the flows. It was
concluded that at a constant Helium flow yield the oscillation frequency of
the rising gas column decreases in from of a power law as a function of the
nozzle diameter. This finding is similar with the observed oscillation
frequency of the flames of candle bundles as a function of the number of
candles in the bundle Gergely  _et al._ (2020). For a constant nozzle diameter
it was found that the oscillation frequency of the flow increases linearly
with the flow yield. For the collective behavior of two nearby and clearly
separated flow columns with similar flow parameters only counter-phase
synchronization was observed. This is somehow different from the phenomena
observed for candle bundle flames, where both in-phase and counter phase
synchronization is present depending on the distance between the flames
Gergely  _et al._ (2020). The experiments concluded that for short distances
the oscillation frequency of the flow column becomes significantly higher than
the frequency observed for non-interacting Helium columns with the same
parameters (flow rate and nozzle diameter). All the above summarized results,
should be a test for any numerical simulation approach on this interesting
phenomenon.

Due to the complexity of the problems related to flows in different spatial
configurations, computation approaches are many times the only theoretical
possibilities to realistically model such phenomena (see for example Joel H.
Ferziger and Street (2020); Wendt (2009)). Even with such a modelling
methodology, imposing the right boundary conditions and offering a proper
discretization of space and time raises many technical challenges V.
Venkatakrishnan and R. Chakravarthy (1998). The incredible revolution we
experience nowadays in computational resources and methods, helped us to
overcome much of these difficulties and computational fluid dynamics became
the primary tool to investigate theoretically problems related to fluid
dynamics. However, even with the presently available computational power, many
times we are forced to investigate a simpler flow topology and reduce the
dimensionality of the problem Correa  _et al._ (2019). This is nowadays a
standard procedure for those cases where the problem becomes computationally
difficult in 3D. Usually a two-dimensional simplification is considered when
the periodicity and symmetry of the considered flow allows this. Assuming in
the followings a cylindrical symmetry for the flow, we consider a two-
dimensional numerical fluid dynamics approach to the above described
phenomenon. First, we discuss the theoretical background on which our approach
is build and the details of the applied numerical method. Using simple and
straightforward examples we test the simulation environment and gain
confidence in the method. After this methodological part, we approach the
proposed problem and compare critically the results of the simulations with
the experimental data from Gergely  _et al._ (2021). Finally, conclusions are
drown and universality features of this intriguing phenomenon is discussed.

##  II The numerical approach

We present here a 2D numerical approach, that is suitable for modelling the
oscillations and collective behavior observed in rising gas columns. In order
to further simplify the problem, instead of a Helium column injected from the
bottom we consider the flow of the same incompressible fluid as the
surrounding, heated in a restricted region at the bottom of the simulated
area. In such manner we get a rising gas column which is also realizable in
experiments.

Using the same Schlieren technique as before, on Figure 1 and on the movies
presented in Gergely (2020b) we show that very similar instabilities and
oscillations occur. In these kitchen-experiments the heating is realized by a
simple heating coil in which one controls the dissipated electric power.
Unfortunately in such experiments there is no good control over the flow
debit, therefore one cannot conduct such carefully monitored experiments as
were done for Helium. This is the main reason why the numerical results are
compared with our earlier experiments presented in Gergely  _et al._ (2021).

The advantage of the proposed setup is that we do not have to apply the
numerical fluid dynamics method for two component gases. We pay however for
this simplification by the non-homogeneous temperature field, therefore extra
transports and gradients has to be taken into account.

Figure 1: Visualization of a rising hot-air column by the Schlieren technique.
Similar instabilities and oscillations appear as in the case of the rising
Helium gas column.

In our approach the fluid is considered to be an ideal one described by the
Navier-Stokes equation. For an incompressible fluid in a gravitational field,
the Navier-Stokes equation writes in the following form:

$\begin{split}\rho\cdot\frac{\partial\mathbf{u}}{\partial
t}+\rho\cdot\left(\mathbf{u}\cdot\nabla\right)\cdot\mathbf{u}&=-\nabla
p+\mathbf{g}\cdot\rho+\mu\cdot\Delta\mathbf{u}\\\
\nabla\mathbf{u}&=0\end{split}$ (1)

Here $\rho$ denotes the density, $p$ is the pressure, $\mathbf{g}$ the
gravitational acceleration, $\mathbf{u}$ is the velocity of the fluid, and
$\mu$ denotes the fluid’s viscosity. The quantities $\mathbf{u}$ and $p$ can
be time and position dependent in the flow-space.

For the considered problem, the convective flow due to the temperature
difference plays a key role, therefore in the Navier-Stokes equation, we have
to take into account the temperature dependence of the density and describe
also the time evolution of the temperature inside the fluid. The evolution of
temperature and the temperature dependence of the density, are considered by
the following equations:

$\begin{split}\rho=\frac{\rho_{0}}{1+(T-T_{0})\cdot\alpha}\\\ \frac{\partial
T}{\partial t}=D\cdot\Delta T-(\mathbf{u}\cdot\nabla)\cdot T\end{split}$ (2)

In the above equations, $\rho_{0}$ is the density at $T_{0}$, $T$ is the
temperature of the fluid at a given spatial position and in a given time-
moment, $D$ is the diffusion constant and $T_{0}$ is the ambient temperature.
The numerical solution of the coupled systems of partial differential
equations (1)-(2) was done by using the "FEniCS" software package Langtangen
and Logg (2017). FEniCS is an open-source platform developed for solving
Partial Differential Equation (PDE) systems. We chose this platform because it
has high-level programming interfaces (C ++, Python), the shape of the
equations in the program code is similar to their symbolic form and the
program is optimized for a wide range of hardware from laptops to high-
performance clusters.

###  II.1 The simulation code

FEniCS uses finite element methods to solve PDEs. As an example in Appendix A
we illustrate how to solve the simple 2D Poison equation with FEniCS. For our
specific problem we first deal with the term describing the evolution of
temperature:

$\begin{split}\frac{\partial T}{\partial t}=D\cdot\Delta
T-(\mathbf{u}\cdot\nabla)\cdot T\end{split}$ (3)

This equation contains a time derivative, so in addition to the coordinates we
also have to discretize time. This is done by the Euler method as follows:

$\begin{split}\frac{T(t+dt)-T(t)}{dt}=D\cdot\Delta
T(t)-(\mathbf{u}\cdot\nabla)\cdot T(t),\end{split}$ (4)

where the indexes for the $T$ temperature label the discretized time-steps. We
then bring each term to the left hand side of the equation, we multiply the
equation by a $\tau$ test function, and integrate the equation over the entire
simulated domain:

$\int_{\Omega}[T(t+dt)-T(t)-D\cdot\Delta T(t)\cdot
dt+(\mathbf{u}\cdot\nabla)\cdot T(t)\cdot dt]\cdot\tau\ d\Omega=0$ (5)

The equations from above contains a second-order derivative for the
coordinates, which is eliminated by partial integration:

$\int_{\Omega}(\nabla^{2}\cdot T(t))\cdot\tau\
d\Omega=\int_{\partial\Omega}\left(\frac{\partial
T(t)}{\partial\mathbf{n}}\right)\cdot\tau\ ds-\int_{\Omega}\nabla
T(t)\cdot\nabla\tau\ d\Omega$ (6)

Here we denoted by $\mathbf{n}$ the unit normal vector to the $\partial\Omega$
surface. The derivative in respect to $\mathbf{n}$ is defined as:

$\frac{\partial T}{\partial\mathbf{n}}=(\nabla T)\cdot\mathbf{n}$ (7)

Rewriting equation (5) using the above result and the fact that under the
Dirichlet and free boundary conditions the surface integral disappears, we
obtain the final form:

$\int_{\Omega}([T(t+dt)-T(t)+(\mathbf{u}\cdot\nabla)\cdot T(t)\cdot
dt]\cdot\tau\ +D\cdot\nabla T(t)\cdot\nabla\tau\cdot dt\ )\,d\Omega=0$ (8)

The incompressible Navier-Stokes equation (1) was solved using the IPCS
(Incremental Pressure Correction Scheme) scheme Logg  _et al._ (2012). The
IPCS method consists of 3 steps but before specifying the steps we introduce
the following functions and notation:

$\begin{split}[\varepsilon(\mathbf{u})]=\frac{1}{2}\cdot([\nabla\otimes\mathbf{u}]+[\nabla\otimes\mathbf{u}]^{T})\\\
[\sigma(\mathbf{u},p)]=2\cdot\mu\cdot\varepsilon(\mathbf{u})-p\cdot I\\\
{\langle\mathbf{f},\mathbf{g}\rangle}_{\Omega}=\int_{\Omega}\mathbf{f}\cdot\mathbf{g}\
d\Omega\\\ {\langle[A],[B]\rangle}_{\Omega}=\int_{\Omega}[A]:[B]\
d\Omega\end{split}$ (9)

The $\otimes$ product is defining a matrix with the following elements:

$\nabla\otimes\mathbf{u}=\left[\frac{\partial u_{j}}{\partial x_{i}}\right]$
(10)

We denoted by $[...]$ a square matrix, by $[...]^{T}$ the transpose of a
matrix and by $:$ the inner product of matrices:

$[A]:[B]\equiv\sum_{i,j}A_{ij}\,B_{ij}$ (11)

Using the $\varepsilon$ and $\sigma$ functions and the specified notation, the
steps of the method will be described in the followings. First we reconsider
the Navier-Stokes equation in a discretised form on the $\Omega$ velocity
space:

$\displaystyle\rho\cdot{\left\langle\frac{\mathbf{u}^{*}-\mathbf{u}(t)}{dt},\mathbf{v}\right\rangle}_{\Omega}+\rho\cdot{\left\langle\mathbf{u}(t)\cdot\nabla\mathbf{u}(t),\mathbf{v}\right\rangle}_{\Omega}+{\left\langle\left[\sigma\left(\frac{\mathbf{u}(t)+\mathbf{u}^{*}}{2},p(t)\right)\right],\left[\varepsilon(\mathbf{v})\right]\right\rangle}_{\Omega}+$
(12) $\displaystyle+{\left\langle
p(t)\cdot\mathbf{n},\mathbf{v}\right\rangle}_{\partial\Omega}-{\left\langle\mathbf{n}\cdot[\nabla\otimes\left(\frac{\mathbf{u}(t)+\mathbf{u}^{*}}{2}\right)]^{T},\mathbf{v}\right\rangle}_{\partial\Omega}=\rho\cdot{\left\langle\mathbf{g},\mathbf{v}\right\rangle}_{\Omega}$

Here $\mathbf{v}$ is a test function, for more information on choosing this
one should consult Langtangen and Logg (2017) The first step of the method is
the calculation of an intermediate velocity $\mathbf{u}^{*}$ from which the
pressure will be determined. Then, the pressure is determined in the i-th step
in equation:

$\begin{split}{\langle\nabla p(t+1),\nabla q\rangle}_{\Omega}={\langle\nabla
p(t),\nabla
q\rangle}_{\Omega}-\frac{{\langle\nabla\mathbf{u}^{*},q\rangle}_{\Omega}}{dt}\cdot\rho\end{split}$
(13)

In the last step, the velocity in the $t+dt$ time step is determined based on
the pressure and the intermediate velocity:

$\begin{split}{\langle\mathbf{u}(t+dt),\mathbf{v}\rangle}_{\Omega}={\langle\mathbf{u}^{*},\mathbf{v}\rangle}_{\Omega}-\frac{dt\cdot{\langle\nabla(p(t+dt)-p(t)),\mathbf{v}\rangle}_{\Omega}}{\rho}\end{split}$
(14)

In the equation from above $q$ is a test function for the pressure. The method
described above for solving the incompressible Navier-Stokes equation is
implemented in 2D. In order to solve the equations numerically, we need
boundary conditions in addition to discretization. We used Dirichlet and free
boundary conditions. The Dirichlet boundary condition means that the value of
the quantity at a given point is fixed. In the case of the free boundary
condition, the derivative as a function of the coordinates of the quantity at
the given point is 0.

We tested visually our 2D simulation environment in two simple problems. First
we intended to reproduce the Karman vortices in the flow of a fluid around an
obstacle (Appendix B). Second, we simulated the expansion and rising of a
heated sphere, verifying the code for non-homogeneous temperature conditions
as well (Appendix C). The test simulations reproduced the expected realistic
behavior for these known problems, giving confidence for the correct
implementation of the relevant equations discussed above.

###  II.2 Simulating the rising hot air column

In the followings we give the details for implementing the simulations aiming
to reproduce the characteristic oscillations observed in a rising gas column.
The boundary conditions introduced for velocity and temperature will be
justified, and we explain how the time series of the characteristic
oscillations were obtained and how the oscillation frequency was calculated.

We consider the inflow geometry presented in Figure 2 leading to the flow
illustrated in Figure 3 (a). On the sidewalls, the value of the velocity is
fixed to 0, on the lower boundary, the $x$-direction component of the velocity
is considered as 0 and the $y$-direction component is given by the following
parabolic-like kernel (see Figure 2)

$\displaystyle
v_{y}(x,0)=c_{1}f\left(x,\frac{d}{2}\right)\left(\frac{d}{2}-x\right)\left(x+\frac{d}{2}\right)+c_{2}f\left(x,\frac{H}{2}\,(1-\frac{1}{30})\right),$

with:

$\displaystyle
f(a,b)=\frac{1}{e^{-c_{3}\cdot(b+a)}+1}-\frac{1}{e^{c_{3}\cdot(b-a)}+1}.$ (16)

Figure 2: An example for the y-component of velocity, $v_{y}(x,0)$, and
temperature profile, $T(x,0)$, of the heated fluid column at the bottom
boundary of the simulated space. The following parameters were used:
$c_{1}=1600\ m^{-1}\cdot s^{-1},\ c_{2}=0.1\ m\cdot s^{-1},\ c_{3}=2000\
m^{-1},\ d=0.08\ m,\ T_{0}=300\,K,\ H=0.3\ m$

In the equation from above, $d$ denotes the nozzle diameter, $H$ denotes the
width of the simulated space, the parameters $c_{1},c_{2}$ determine the
incoming flow rate of the fluid, and $c_{3}$ is a tuning parameter governing
the cut in each profile.

At the upper boundary (height $L$), free boundary conditions are applied for
the $y$ component of the velocity and for the $x$ component Dirichlet
condition is applied, i.e. $v_{x}(x,L)=0$. For pressure, Dirichlet boundary
conditions were used in the upper part of the simulated volume,
$p(x,L)=g\cdot\rho_{0}\cdot l$, and free boundary condition for the other
boundaries. The temperature on the walls is fixed to $T_{0}$. On the upper
boundary we consider $T_{0}$ if the $y$-direction component of the velocity is
negative, otherwise free boundary conditions are used. The temperature at the
lower boundary is determined by using the equation:

$\displaystyle T(x,0)=T_{0}+T_{heating}f\left(x,\frac{d}{2}\right)$ (17)

Here, the $T_{heating}$ temperature governs the extra heat amount introduced
in the system. For simplicity reasons we have used in all the presented
results $T_{heating}=T_{0}$. In the first attempts at the upper part of the
simulated box, free boundary condition were considered for the velocity.
However in such cases unexpected instabilities occurred and after a certain
time the heated fluid column was pushed on one of the sidewalls. We have
examined carefully this phenomenon and concluded that a self-amplifying effect
is responsible for its development. Due to the convective flow, the amount of
fluid leaving the simulation box is larger than the volume of fluid flowing
into the simulation box through the lower boundary. Since the fluid is
incompressible, fluid must flow back into the simulation box through the upper
boundary. Since there is always an asymmetry in the profile of the fluid
inflow this will slightly deflect the outflowing column. In the direction of
the deflection, the inflow area decreases, so the asymmetry in the fluid
inflow increases. An increase in asymmetry over time will result in the fluid
flowing along one of the walls. This is the simple explanation of the observed
instabilty.

Two methods were used to eliminate these instabilities. The first method is to
flow a fluid of ambient temperature $T_{0}$ at a constant rate on both sides
of the heated air column. Since the flow is two-dimensional, the fluid flowing
on a given side can only leave on the same side and this will always provide a
minimum distance from the wall for the rising jet. The second method is to
allow only the $y$-direction component of the velocity at the upper boundary.
Combining this two methods will eliminate the tendency of the jet to approach
one of the sidewalls.

For the upper boundary a proper boundary condition has to be applied also for
the inflowing fluid temperature. At the upper boundary an inflow is also
necessary in order to respect the incompressibility of the fluid. Since the
temperature of the outflowing fluid varies over a wide range we can’t apply
Dirichlet boundary condition to the whole upper boundary because this would
cause unmanageable gradients. Avoiding large gradients due to large
temperature differences was solved by applying boundary condition only to
those points where the $y$-direction component of the velocity became
negative.

The time series for the relevant hydrodynamical parameters were generated and
we followed the temperature distribution in the simulated space. For
realistically chosen parameters, it was shown that the model is capable to
produce an oscillation similar to the one observed in the case of the Helium
column. Interestingly, it was found that such oscillations are possible even
on largely different length-scales. The observed oscillation is shown in
Figure 3 a, b. where we illustrate the temperature space at subsequent time
moments. For the simulations presented in Figure 3, we used the parameters
spcified in the figure caption.

(a)

(b)

Figure 3: Oscillation of a heated air column in snapshots. The images show the
temperature space at the specified $t$ time moments for two different length
scales (Figures a. and b.). The gravity acts in the negative direction of the
$y$-axis and the parameters of the simulation were chosen as follows: (a)$\
\alpha=0.33\cdot 10^{-2}\ K^{-1},\,\rho_{0}=1.2\ Kg\cdot m^{-2},\,T_{0}=300\
K,\,D=10^{-4}\ m^{2}\cdot s^{-1}\cdot K^{-1},\,g_{y}=-9.81\ m\cdot
s^{-2},\,\mu=1.96\cdot 10^{-5}\ Kg\cdot s^{-1},\,c_{2}=0.1\ m\cdot
s^{-1},\,c_{3}=2000\ m^{-1},\,d=0.08\ m,\,c_{1}=1600\ m^{-1}\cdot s^{-1}$,
(b) $\alpha=10^{-3}\ K^{-1},\,\rho_{0}=1\ Kg\cdot m^{-2},\,T_{0}=300\
K,\,D=5\cdot 10^{-2}\ m^{2}\cdot s^{-1}\cdot K^{-1},\,g_{y}=-9.81\ m\cdot
s^{-2},\,\mu=5\cdot 10^{-2}\ Kg\cdot s^{-1},\,c_{2}=0.4\ m\cdot
s^{-1},\,c_{3}=5\ m^{-1},\,d=8\ m,\,c_{1}=0.375\ m^{-1}\cdot s^{-1}$

For a quantitative evaluation of the simulated dynamics, the Otsu method was
applied for the 2D temperature field. To obtain the time series, in uniform
time intervals the Otsu processed pixels were summed up to a certain height,
after this the obtained time series was divided by the average value of the
time series. The oscillation frequency was calculated in a similar manner with
the experiments, based on the above generated time series. In the first step a
Fourier transform was applied to the time series and then the value of the
frequency belonging to the largest peak was determined as the relevant
oscillation frequency.

With the implemented simulation code we examined how the inflow rate (yield)
of the heated fluid column and the nozzle diameter affects its oscillation
frequency. We also investigated the collective behavior for the oscillation of
two columns placed nearby each other.

###  II.3 Numerical results for the oscillation frequency

The effect of flow yield and nozzle diameter was examined on two different
length-scales. To study the flow yield we used the parameter sets (a), (b)
introduced above and the nozzle diameters were $d=0.08\ m$ and $d=8\ m$
respectively. For constant $c_{3}$, $c_{2}$, and $d$ parameters the yield
(flow debit) of the heated fluid depends only on $c_{1}$:

$\displaystyle\Phi=\int_{-\frac{d}{2}}^{\frac{d}{2}}v_{y}(x,0)\ dx=$ (18)
$\displaystyle=\int_{-\frac{d}{2}}^{\frac{d}{2}}[c_{1}f\left(x,\frac{d}{2}\right)\left(\frac{d}{2}-x\right)\left(x+\frac{d}{2}\right)+c_{2}f\left(x,\frac{H}{2}(1-\frac{1}{30})\right)]dx$

The computed oscillation frequency of the heated fluid column as a function of
the $\Phi$ parameter is plotted in Figure 5a and Figure 4a. One will observe
that the oscillation frequency increases as the flow rate $\Phi$ increases and
this increasing trend can be well approximated by a linear fit in good
agreement with the experimental results plotted in Figure 6.

Figure 4: Simulation results for the smaller length-scale. Figure (a) shows
the oscillation frequency of the heated fluid column as a function of the flow
yield $\Phi$ fixed by equation (18) for $d\ =\ 0.08\ m$ inflow diameter.
Figure (b) shows the oscillation frequency of the heated fluid column as a
function of the $d$ nozzle diameter. The other parameters used are the same as
the ones specified in the caption of Figure 3, the value of $c_{1}$ for the
different nozzle diameters are given in Table 1.

Figure 5: Simulation results for the larger length-scale. Figure (a) shows the
oscillation frequency of the heated fluid column as a function of the flow
yield $\Phi$ fixed by equation (18) for $d\ =8\ m$ inflow diameter. Figure (b)
shows the oscillation frequency of the heated fluid column as a function of
the $d$ nozzle diameter. The other parameters are the same as the ones
specified in the caption of Figure 3, the value of $c_{1}$ for the different
nozzle diameters are given in Table 2. Figure 6: Experimentally observed
oscillation frequency of a Helium column as a function of the yield (flow
debit) obtained for a setup with nozzle diameter of 2 cm. With increasing flow
yield, the frequency of the oscillation increases in an almost linear manner.
The plot is done by using the results presented in Gergely  _et al._ (2021).
Error bars show the uncertainty in the measurement of the flow yield.

The effect of nozzle diameter on the oscillation frequency was investigated at
a constant inflow yield. Since the yield $\Phi$ depends on $d$ according to
equation (18), for different nozzle diameters we must rescale the parameters
$c_{1}$ so that the flow rate remains constant. For the smaller length-scale
simulations, we used $\Phi_{1}=0.076\ m^{2}/s$ flow yield and for the larger
scale simulations, we used $\Phi_{2}=29\ m^{2}/s$ flow yield. To keep constant
the flow yield for different nozzle diameters, we varied the value of the
$c_{1}$ parameter. The $c_{1}$ values for the different nozzle diameters are
shown in Table 1 and Table 2.

d [m] | 0.04 | 0.06 | 0.08 | 0.1 | 0.12 | 0.14
---|---|---|---|---|---|---
$c_{1}\ [m^{-1}\cdot s^{-1}]$ | 6781 | 1952 | 800 | 398 | 223 | 136
Table 1: Value of the $c_{1}$ parameter for different nozzle diameters $d$, in order to keep the flow rate $\Phi_{1}=0.076\ m^{2}/s$. d [m] | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---
$c_{1}\ [m^{-1}\cdot s^{-1}]$ | 0.45 | 0.3 | 0.21 | 0.153 | 0.11 | 0.088
Table 2: Value of the $c_{1}$ parameter for different nozzle diameters $d$, in
order to keep the flow rate $\Phi=29\ m^{2}/s$.

For both length scales a decreasing trend of the oscillation frequency as a
function of the nozzle diameter was observed. Results in such sense are
plotted in Figure 5b and Figure 4b, the trend is in good agreement with the
experimental results from Figure 7.

Figure 7: Experimentally observed oscillation frequency of a Helium column as
a function of the nozzle diameter for a yield of $\Phi=$46 $\pm$ 2.3 $l/min$.
The plot is done by using the results presented in Gergely  _et al._ (2021).
Error bars show the variations of the experimentally used nozzle diameter in
different directions.

###  II.4 Numerical results for the collective behavior

We turn now our interest on reproducing the experimentally observed collective
behaviour in form of anti-phase synchronization.

The dimensions of the simulation boxes used to study the collective behavior
are the following: 46 m wide (H=46 m) and 30 m high (L=30 m) for the large
length-scale and 0.3 m wide (H=0.3 m) and 0.15 m high (L=0.15 m) at the
smaller length-scale. At the lower boundary, the $x$ component of the inflow
fluid velocity is 0, and the $y$ component is given by the following kernel
function:

$\displaystyle
v_{y}(x,0)=c_{1}f\left(x-d_{0},\frac{d}{2}\right)\left(d_{0}-x+\frac{d}{2}\right)\left(x-d_{0}+\frac{d}{2}\right)+$
$\displaystyle+c_{1}f\left(x+d_{0},\frac{d}{2}\right)\left(\frac{d}{2}-d_{0}-x\right)\left(x+d_{0}+\frac{d}{2}\right)+c_{2}f\left(x,\frac{H}{2}-c_{4}\right)$

Here the value of $c_{4}$ is $0.5\ m$ for the large scale system and $0.005\
m$ for the small scale system.

Figure 8: An example for the y-component velocity, $v_{y}(x,0)$, and
temperature, $T(x,0)$, profiles of the heated fluid columns at the bottom
boundary of the simulated space. The following parameters were
used:$c_{1}=0.45\ m^{-1}\cdot s^{-1},\ c_{2}=0.4\ m\cdot s^{-1},\ c_{3}=5\
m^{-1},\ d=7\ m,\ d_{0}=6\ m,\ T_{0}=300\ K,\ H=46\ m$.

This leads to an inflow profile with two peaks where the centers are separated
at a distance of $2d_{0}$, as it is illustrated in Figure 8. The temperature
profile is adjusted accordingly:

$\displaystyle
T(x,0)=T_{0}+T_{heating}f\left(x-d_{0},\frac{d}{2}\right)+T_{heating}f\left(x+d_{0},\frac{d}{2}\right)$

We used the same simulation parameters as before and fixed $d=7\ m$,
$c_{1}=0.45\ m^{-1}\cdot s^{-1}$ values for the large length-scale and
$d=0.04\ m$, $c_{1}=6400\ m^{-1}\cdot s^{-1}$ values for the small length-
scale system. Again, for the presented results we considered
$T_{heating}=T_{0}$. Experimental results from Gergely  _et al._ (2021) shows
that at small separation distance collective behavior in form of counter-phase
synchronization appears. A snapshot for a simulated stable collective behavior
is visible in Figure 9, reproducing successfully this counter phase
synchronization on the smaller length-scale. Similar behavior is observable
for the larger length-scales as well.

Figure 9: Counter-phase synchronization of two nearby heated columns. Computer
simulation results with the following parameters: $\ \alpha=0.33\cdot 10^{-2}\
K^{-1},\,\rho_{0}=1.2\ Kg\cdot m^{-2},\,T_{0}=300\ K,\,D=10^{-4}\ m^{2}\cdot
s^{-1}\cdot K^{-1},\,g_{y}=-9.81\ m\cdot s^{-2},\,\mu=1.96\cdot 10^{-5}\
Kg\cdot s^{-1},\,c_{2}=0.1\ m\cdot s^{-1},\,c_{3}=2000\ m^{-1},\,d=0.04\
m,\,c_{1}=6400\ m^{-1}\cdot s^{-1},\ H=0.3\ m$

For the pictures processed with the Otsu method the collective oscillation of
nearby heated fluid columns are shown in Figures 10 a and b, for the small and
large length-scales, respectively.

Figure 10: Simulated time series for the oscillations of nearby heated fluid
columns. Motion of the interface is detected by the Otsu method at the same
height from the nozzle. For not too high separation distances ((a) $2\cdot
d_{0}=0.03\ m$ , (b) $2\cdot d_{0}=5\ m$) a clear counter-phase
synchronization is observable. The parameters of the simulations are
(a): $\ \alpha=0.33\cdot 10^{-2}\ K^{-1}$, $\rho_{0}=1.2\ Kg\cdot m^{-2}$,
$T_{0}=300\ K$, $D=10^{-4}\ m^{2}\cdot s^{-1}\cdot K^{-1}$, $g_{y}=-9.81\
m\cdot s^{-2}$, $\mu=1.96\cdot 10^{-5}\ Kg\cdot s^{-1}$, $c_{2}=0.1\ m\cdot
s^{-1}$, $c_{3}=2000\ m^{-1}$, $d=0.04\ m$, $c_{1}=6400\ m^{-1}\cdot s^{-1}$,
$2\cdot d_{0}=0.03\ m$.
(b): $\alpha=10^{-3}\ K^{-1}$, $\rho_{0}=1\ Kg\cdot m^{-2}$, $T_{0}=300\ K$,
$D=5\cdot 10^{-2}\ m^{2}\cdot s^{-1}\cdot K^{-1}$, $g_{y}=-9.81\ m\cdot
s^{-2}$, $\mu=5\cdot 10^{-2}\ Kg\cdot s^{-1}$, $c_{2}=0.4\ m\cdot s^{-1}$,
$c_{3}=5\ m^{-1}$, $d=8\ m$, $c_{1}=0.375\ m^{-1}\cdot s^{-1}$, $2\cdot
d_{0}=0.03\ m$ and we have fixed $d=7$ and $c_{1}=0.45\ m^{-1}\cdot s^{-1}$.

For the indicated smaller separation distances an almost perfect counter-phase
synchronization develops. For larger separation distances the phases of the
oscillations will begin to shift relative to each other and no clear phase-
difference blocking is observable.

For the simulations performed on the smaller length-scale, corresponding to
the experimental conditions in Gergely  _et al._ (2021) we computed the
synchronization order parameter which is meant to characterize the collective
oscillation. We used the same $z$ synchronization order parameter as the one
used in Gergely  _et al._ (2020) and Gergely  _et al._ (2021). The
computationally derived synchronization parameter is plotted in Figure 11 a.
It’s values in the neighborhood of -1 indicates that we have counter-phase
synchronization for the studied distances. In Figure 11b we also show the
oscillation frequency of the two synchronized heated fluid columns as a
function of their separation distance. This frequency decreases as we increase
the separation distance between the columns, similarly with what has been
reported in our experiments for the Helium columns Gergely  _et al._ (2021).

Figure 11: Simulation results for the collective behavior of two heated
columns. Figure (a) the synchronization order parameter of two interacting
heated fluid columns and Figure (b) shows the collective oscillation
frequency, both as a function of the separation distance between the columns.
The following simulation parameters were considered: $\ \alpha=0.33\cdot
10^{-2}\ K^{-1},\,\rho_{0}=1.2\ Kg\cdot m^{-2},\,T_{0}=300\ K,\,D=10^{-4}\
m^{2}\cdot s^{-1}\cdot K^{-1},\,g_{y}=-9.81\ m\cdot s^{-2},\,\mu=1.96\cdot
10^{-5}\ Kg\cdot s^{-1},\,c_{2}=0.1\ m\cdot s^{-1},\,c_{3}=2000\
m^{-1},\,d=0.04\ m,\,c_{1}=6400\ m^{-1}\cdot s^{-1}$

##  III Discussion and Conclusions

In our previous study Gergely  _et al._ (2021), we used both experimental and
theoretical approaches to investigate whether the hydrodynamic instabilities
that occur in rising gas columns is responsible also for the oscillations
observed for the diffusion flames Gergely  _et al._ (2020). It was shown that
indeed this is the case: Helium columns ascending in air from a circular
nozzle produce similar oscillations with the ones observed in diffusion
flames. Also, the similar collective behavior of these oscillations (counter-
phase synchronization) for Helium columns and flickering candle flames suggest
that the hydrodynamic processes by their own are enough to explain these
phenomenon.

In our previous work Gergely  _et al._ (2021), for modeling the observed
oscillations , a simplified but analytically treatable hydrodynamical approach
was used. The model predicted the right trends for the oscillation frequencies
as a function of the relevant parameters, but was unsuitable to approach the
collective behavior. In the present follow-up, we offered an improved modeling
by considering a 2D numerical hydrodynamics computer simulation where for
computational simplicity instead of ascending Helium columns, heated fluid
columns were considered. This approach proved to be successful for reproducing
all the experimentally observed features. For a constant nozzle diameter the
numerical hydrodynamics approach lead to an oscillation frequency that
increased roughly linearly with the flow yield, in agreement with the
experimental results. For constant flow yield the numerical results suggested
a decreasing trend of the oscillation frequency as a function of the nozzle
diameter, confirming the experimental results. The exact shape of the
simulated trend was slightly different however from the one observed in the
experiments. The main reason for this discrepancy is probably the reduction of
the real 3D problem to a 2D topology. Finally, the presented computer
simulations were successful also in reproducing the counter-phase
synchronization of two heated fluid columns placed nearby each other. The
computed trends for the synchronization order-parameter and the collective
frequency were also in agreement with the experimental results obtained for
rising Helium columns in air.

It worth mentioning here that the computer simulations were performed both for
laboratory and for a much larger length-scale than the experiments. The
qualitative agreement between the results (trends and collective behavior) on
these different length-scales suggests that the investigated phenomenon is
more general than it was thought to be, and might have further, yet unexplored
connections.

Acknowledgment We acknowledge Dr. M. Ercsey-Ravasz, Dr. B. Molnár, and Prof.
dr. G.C. Silaghi, all from the Babes-Bolyai University in helping us with the
computational resources for running our simulation codes.

## Appendix

### A. Solving the 2D Poison equation in FEniCS

Here we illsutrate through the 2D Poison equation how to solve a PDE in
FEniCS. The Poisson equation can be given in the following form:

$-\Delta\varphi=f$ (21)

If we have a simple rectangular space then the above equation can be easily
given in the finite element form

$-\frac{\varphi_{i-1,j}-2\varphi_{i,j}+\varphi_{i+1,j}}{h^{2}}-\frac{\varphi_{i,j-1}-2\varphi_{i,j}+\varphi_{i,j+1}}{h^{2}}=f_{i,j},$
(22)

however, with this simple and intuitive approach, we soon run into problems
because even for a circle it is impossible to map the boundary with an
acceptably small number of squares. The FEniCS program Langtangen and Logg
(2017) uses a triangular grid instead of a square grid to cover the simulated
space, in which case we can always select the grid so that the grid points are
on the boundary surfaces. Discretization alone does not solve the equation,
the next question is how to determine the solution at each lattice point. As a
first step, we write the $\varphi(x_{1},x_{2})$ function in the following
form:

$\varphi(x_{1},x_{2})=\sum_{i=0}^{N}c_{i}\phi_{i}(x_{1},x_{2})$ (23)

In the above equation, $\phi_{i}(x_{1},x_{2})$ is a given k-th order
polynomial, $c_{i}$ are the coefficients that determine $\varphi(x_{1},x_{2})$
and N+1 is the number of the grid points. The $c_{i}$ coefficients are
determined by multiplying the Poison equation by N+1 different
$v(x_{1},x_{2})$ so-called test functions and integrating the product over the
whole domain to obtain N + 1 linearly independent equations from which the
$c_{i}$ coefficients can be calculated. All this can be formally given in the
following form:

$\int_{\Omega}-v\,\Delta\varphi\ d\Omega=\int_{\Omega}fv\ d\Omega$ (24)

We denoted by $d\Omega=dx_{1}dx_{2}$.The above form of PDE is called the weak
formulation of the equation and this is what is calculated by the FEniCS
program. The second-order coordinate derivative in the above equation means
that the polynomials used need to be twice differentiable. Because the use of
polynomials with large degrees requires more memory and computation, we always
strive to keep the degree of polynomials to a minimum. In the above equation,
the reduction of the order of derivatives can be done by Gauss-Green
integration as follows:

$\int_{\Omega}-v\,\Delta\varphi\ d\Omega=\int_{\Omega}\nabla v\nabla\varphi\
d\Omega-\int_{\partial\Omega}\frac{\partial\varphi}{\partial{\bf n}}v\ d{s}$
(25)

Since we use Dirichlet boundary conditions for the Poisson problem, the value
of $v$ at the boundary is 0, so the equation (24) can be written in the
following form:

$\int_{\Omega}\nabla v\nabla\varphi\ d\Omega=\int_{\Omega}f\,v\ d\Omega$ (26)

We have seen above how to rewrite the 2D Poisson problem in a form that can be
solved with the FEniCS program, now we show the implementation of the solution
in Python.

⬇

from fenics import * #

import numpy as np #Numpy is required for error calculation

import matplotlib.pyplot as plt #We plot the result with the matplotlib

Nx=10#The number of grid points in the x directions

Ny=10#The number of grid points in the y directions

mesh=UnitSquareMesh(Nx,Ny)

V=FunctionSpace(mesh,’P’,1)#space containing first degree polynomials

fi_D=Expression(’1+x[0]*x[0]+2*x[1]*x[1]’,degree=2)

#boundary conditions equation (on boundarys fi(x,y)=x^2+2y^2+1)

def boundery(x,on_boundary):

return on_boundary

bc=DirichletBC(V,fi_D,boundery)#boundary conditions

fi=TrialFunction(V)

v=TestFunction(V)

f=Constant(-6)

a=dot(grad(fi),grad(v))*dx#right side of equation

L=f*v*dx#left side of equation

fi=Function(V)

solve(a==L,u,bc)#solve the equation

c = plot(interpolate(fi, V), mode=’color’)

plt.colorbar(c)

plot(fi)

plt.savefig(’result1.png’)

plt.show()

vertex_v_ud=fi_D.compute_vertex_values(mesh)

vertex_v_u=fi.compute_vertex_values(mesh)

err_max=np.max(np.abs(vertex_v_ud-vertex_v_u))

print("maximum␣error:␣",err_max)

The above program solves equation (21) on the unit square of
$\\{(0,0),(1,1)\\}$. The largest difference between the theoretically expected
and the numerically obtained value was of the order of the precision of the
numerical representation of the numbers, giving us confidence for the use of
the numerical solution.

### B. Test for the 2D hydrodynamics simulations. Karman vortices

The first phenomenon we aimed to reproduce using our fluid dynamics simulation
is the formation of Karman vortices in the flow of fluids around an obstacle.
With this test we aimed to check visually whether the Navier-Stokes equation
has been correctly implanted, since the temperature of the fluid at all points
is considered fixed: $T_{0}$. Therefore for this test the density in the
simulated volume is constant.

In these simulations, the length of the simulated volume was considered as
$2.2$m, the width of the simulated volume is $0.41$m and in the middle of the
simulated coordinate space a circular obstacle with a radius of $0.05$m is
placed. The coordinates of the centre of this obstacle was taken at (0.2m,
0.2m). The density of the fluid was taken as unity (1kg/m2), the viscosity is
$0.001$kg/s, and no gravitational field is considered. For input (the left
region of the space in Figure 12) we considered that the velocity in the $y$
direction is $0$, the velocity in the $x$ direction has a parabolic profile
with a maximum value of $1.5$m/s. On the horizontal walls and on the boundary
of the obstacle we consider no-slip conditions, thus the velocity is fixed to
$0$. For the output (right-sight region) we have also imposed for the $y$
direction velocity to be zero, and the pressure at the output is fixed also to
$0$. Elsewhere there are all free boundary conditions for the pressure. The
temperature is fixed at $T_{0}=300$K at input and on the horizontal walls.

The velocity vector spaces obtained from the simulations are shown in Figure
12 for $4$ time moments as it is indicated on the left side of the images. One
can observe that the simulation reproduces successfully the expected Karman
vortices.

Figure 12: Velocity fields obtained from the simulation of the Karman vortices
at four consecutive time moments.

### C. Test for the 2D hydrodynamics simulations. Heat induced mushroom cloud

In this second test we aimed to implement the density and temperature
evolution of a heated gas sphere in gravitational field. It is expected that
the shape of the heated gas will follow the known dynamics of a mushroom cloud
in a nuclear explosion.

In the performed simulations, the density at $T_{0}$ was chosen as unity
($1$kg/m2), the ambient temperature $T_{0}$ is $300$K, the $\alpha$ parameter
in equation (2) is $\alpha=0.001\ K^{-1}$, the thermal diffusion constant $D$
is $0.3$m2/s, the initially heated sphere temperature is $600$K, the fluid
viscosity is taken as $\mu=0.05$kg/s, the gravitational acceleration is
$g_{y}=-9.81$m/s2 and the size of the simulated volume is $30$m2 in both the
$x$ and $y$ directions.

At the bottom and walls of the simulated space, the velocity is fixed to $0$,
at the upper boundary, we fix the $x$ component of the velocity to $v_{x}=0$.
For the pressure, the value of $g\cdot\rho_{0}\cdot l$ is fixed for the upper
boundary, and free boundary conditions are applied to all the other
boundaries. For temperature, free boundary condition is applied at the upper
boundary if the $y$ direction component of the velocity is positive otherwise
the temperature is fixed to $T_{0}=300$K units on the upper boundary and on
the side-walls. The temperature is fixed to a higher value of $450$K units at
the bottom-wall of the simulated area. This is necessary in order to make the
resulting flow visible in the temperature space. Initially, the velocity in
the whole simulated volume is $0$ and the volume contains a sphere (disk) with
a radius of $5$m in which the temperature is $600$K. The centre of the sphere
is at the coordinates (0m, 7m), the temperature around the sphere is fixed to
$300$K, and the temperature between the centre and the surface of the disk is
given by an interpolation with a sigmoid function.

Figure 13: Snapshots of the evolution of a heated gas sphere in time. The
parameters and details of the simulation can be found in the text.

The time-evolution of the temperature map derived from the simulation is shown
in Figure 13. The effect of thermal diffusion can be observed in the first two
frames. As a result of this diffusion the initially sharp boundary line
between the high and low temperature regions gets blurred. Subsequent frames
show the displacement due to convective flow and as a result of this the
characteristic mushroom cloud shape is formed.

## References

  * Gergely  _et al._ (2021) A. Gergely, C. Paizs, R. Tötös, and Z. Néda, “Oscillations and collective behavior in convective flows,” Physics of Fluids 33, 124104 (2021).
  * Yuan, Durox, and Villermaux (1994) T. Yuan, D. Durox, and E. Villermaux, “An analogue study for flame flickering,” Experiments in Fluids 17, 337–349 (1994).
  * Chamberlin and Rose (1948) D. Chamberlin and A. Rose, “The flicker of luminous flames,” Proceedings of the Symposium on Combustion 1-2, 27–32 (1948).
  * Durox  _et al._ (1995) D. Durox, T. Yuan, F. Baillot, and J. Most, “Premixed and diffusion flames in a centrifuge,” Combustion and Flame 102, 501–511 (1995).
  * Durox, Yuan, and Villermaux (1997) D. Durox, T. Yuan, and E. Villermaux, “The effect of buoyancy on flickering in diffusion flames,” Combustion Science and Technology 124, 277–294 (1997).
  * Huang  _et al._ (1999) Y. Huang, Y. Yan, G. Lu, and A. Reed, “On-line flicker measurement of gaseous flames by image processing and spectral analysis,” Measurement Science and Technology 10, 726–733 (1999).
  * Kitahata  _et al._ (2009) H. Kitahata, J. Taguchi, M. Nagayama, T. Sakurai, Y. Ikura, A. Osa, Y. Sumino, M. Tanaka, E. Yokoyama, and H. Miike, “Oscillation and synchronization in the combustion of candles,” The Journal of Physical Chemistry A 113, 8164–8168 (2009).
  * Ghosh  _et al._ (2010) S. Ghosh, S. Mondal, T. Mondal, A. Mukhopadhyay, and S. Sen, “Dynamic characterization of candle flame,” International Journal of Spray and Combustion Dynamics 2, 267–284 (2010).
  * Okamoto  _et al._ (2016) K. Okamoto, A. Kijima, Y. Umeno, and H. Shima, “Synchronization in flickering of three-coupled candle flames,” Scientific Reports 6 (2016), 10.1038/srep36145.
  * Chen  _et al._ (2019) T. Chen, X. Guo, J. Jia, and J. Xiao, “Frequency and phase characteristics of candle flame oscillation,” Scientific Reports 9 (2019), 10.1038/s41598-018-36754-w.
  * Gergely  _et al._ (2020) A. Gergely, B. Sándor, C. Paizs, R. Tötös, and Z. Néda, “Flickering candle flames and their collective behavior,” Scientific Reports 10 (2020), 10.1038/s41598-020-78229-x.
  * Settles (2001) G. S. Settles,  _Schlieren and Shadowgraph Techniques : Visualizing Phenomena in Transparent Media_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2001).
  * Leptuch and Agrawal (2006) P. A. Leptuch and A. K. Agrawal, “High-speed rainbow schlieren visualization of an oscillating helium jet undergoing gravitational change,” Journal of Visualization 9, 101–109 (2006).
  * Otsu (1979) N. Otsu, “A threshold selection method from gray-level histograms,” IEEE Transactions on Systems, Man, and Cybernetics 9, 62–66 (1979).
  * FreddyCree (2020) FreddyCree,  _Wikipedia, Otsu’s method_ (2017 (accessed July 11, 2020)).
  * Gergely (2020a) A. Gergely, (2020 (accessed November 10, 2020)a), https://youtube.com/playlist?list=PLamJmxTyiR_3sy-fYUEDEXW4NtbsnsXnW.
  * Joel H. Ferziger and Street (2020) M. P. Joel H. Ferziger and R. L. Street,  _Computational Methods for Fluid Dynamics, 4th edition_ (Springer International Publishing, 2020).
  * Wendt (2009) J. F. Wendt,  _Computational Fluid Dynamics_ (Springer Berlin Heidelberg 2009, 2009).
  * V. Venkatakrishnan and R. Chakravarthy (1998) M. D. S. V. Venkatakrishnan and e. R. Chakravarthy,  _Barriers and Challenges in Computational Fluid Dynamics_ (Springer Science+Business Media Dordrecht, 1998).
  * Correa  _et al._ (2019) P. G. Correa, J. MacIntyre, J. Gomba, M. Cachile, J. Hulin, and A. Auradou, “Three-dimensional flow structures in x-shaped junctions: Effect of the reynolds number and crossing angle,” Physics of Fluids 31, 043606 (2019).
  * Gergely (2020b) A. Gergely,  _Heated air column_ (2020 (accessed November 10, 2020)b), https://youtube.com/playlist?list=PLamJmxTyiR_0fNk5bzTD5GyRF_dOrmB2E.
  * Langtangen and Logg (2017) H. P. Langtangen and A. Logg,  _Solving PDEs in Python_ (Springer, 2017).
  * Logg  _et al._ (2012) A. Logg, K.-A. Mardal, G. N. Wells, _et al._ ,  _Automated Solution of Differential Equations by the Finite Element Method_ (Springer, 2012).


# Multi-task Supervised Learning via Cross-learning ††thanks: $1$: University
of Pennsylvania, United States, e-mails: $\\{$jcervino,
<EMAIL_ADDRESS>$2$: Universidad de la República, Uruguay,
e-mail<EMAIL_ADDRESS>$3$: Massachusetts Institute of Technology,
e-mail<EMAIL_ADDRESS>This work is supported by NSF-Simons MoDLTheorinet
and Uruguay’s ANII FSE 1-2019-1-157459.

Juan Cerviño1, Juan Andrés Bazerque2, Miguel Calvo-Fullana3 and Alejandro
Ribeiro1

###### Abstract

In this paper we consider a problem known as multi-task learning, consisting
of fitting a set of classifier or regression functions intended for solving
different tasks. In our novel formulation, we couple the parameters of these
functions, so that they learn in their task specific domains while staying
close to each other. This facilitates cross-fertilization in which data
collected across different domains help improving the learning performance at
each other task. First, we present a simplified case in which the goal is to
estimate the means of two Gaussian variables, for the purpose of gaining some
insights on the advantage of the proposed cross-learning strategy. Then we
provide a stochastic projected gradient algorithm to perform cross-learning
over a generic loss function. If the number of parameters is large, then the
projection step becomes computationally expensive. To avoid this situation, we
derive a primal-dual algorithm that exploits the structure of the dual
problem, achieving a formulation whose complexity only depends on the number
of tasks. Preliminary numerical experiments for image classification by neural
networks trained on a dataset divided in different domains corroborate that
the cross-learned function outperforms both the task-specific and the
consensus approaches.

###### Index Terms:

Supervised learning, multi-task learning, optimization.

##  I Introduction

Supervised learning is one of the traditional problems studied in statistical
learning. At its core, it consists of learning a function mapping inputs to
outputs based on provided input-output pairs [1]. In practical situations, it
is often the case that in order to learn this function (which we also refer as
a _task_), we have limited input-output pairs available (called _training
data_). It is also common for different tasks to be related in some sense, and
thus, one could attempt to exploit this relationship to improve the individual
tasks’ performance. The study of this problem, in general is known as _multi-
task learning_ [2]. While expected in the limited sample regime, it is also
applicable in the asymptotic case, as exemplified by Stein’s paradox [3],
which shows that when estimating more than three parameters from Gaussian
random variables, a combined estimator exists that has lower Mean Square Error
(MSE) than any separate estimator, even if the random variables are
independent.

Multi-task learning has been widely used in practice, specially finding
success in computer vision [4, 5, 6], and language applications [7, 8].
Applications aside, approaches to multi-task learning can be grouped in two
categories, based on whether they explicitly model or not the relationship
between tasks [9]. Out of those assuming a priori knowledge of the tasks,
prevalent approaches are based on a commonly shared underlying task
representation [10]. This underlying task representation often takes the form
of a sparse, low rank representation across tasks [11, 12], or a manifold
representation [13]. Some other approaches do not assume previous knowledge,
and learn the relationship between tasks directly from data, generally by
performing clustering [14]. In the context of support vector machines, some
works have bounded the pair-wise difference between the elements of the
classifiers’ weights [15, 16]. Some other works focus on finding Pareto
optimal solutions between tasks [17].

In this paper, we take a constrained approach to the multi-task learning
problem. Our formulation is based on the cross-learning framework [18, 19].
Originally used for learning policies in a reinforcement learning scenario,
its principles can be also applied to the multi-task supervised learning
problem. The proposed cross-learning approach consists on bounding the
distance between the parameters of the learned functions (in our case,
classifiers or regression functions). We show, via an illustrative Gaussian
sample mean example, that the use of this centrality measure can guarantee to
outperform both the separate and consensus approaches. Furthermore, since the
cross-learning method is constrained, it requires a projection at each step
which can be computationally expensive (e.g., for neural networks with
millions of parameters). To overcome this hurdle, we propose to solve the
projection step on the dual domain, which reduces the dimension of the problem
to the number of tasks. Finally, we present numerical results in a complex
classification task, which show that our cross-learning strategy outperforms
the consensus and agnostic estimators.

##  II Problem Formulation

We address the problem of learning $N$ functions from data coming from
different tasks. Let function $f:{\mathcal{X}}\times\Theta\to{\mathcal{Y}}$ be
the map between input space ${\mathcal{X}}\subset{\mathbb{R}}^{P}$ and output
space ${\mathcal{Y}}\subset{\mathbb{R}}^{Q}$ parameterized by
$\theta\in\Theta\subset{\mathbb{R}}^{S}$. We seek to minimize the loss
function $\ell$ given $N$ datasets containing pairs
$(x_{j},y_{j})\in({\mathcal{X}},{\mathcal{Y}})$ drawn according to a joint
probability $p_{i}(x,y)$, for $j=1,\ldots,M_{i}$, where $M_{i}$ and
$p_{i}(x,y)$ are the number of samples and generating probability
corresponding to each of the $i=1,\ldots,N$ tasks. Our objective is to obtain
functions $f(x,\theta_{i})$ that stay close together in the optimization space
via the introduction of a measure of centrality $\epsilon$, and central
parameters $\theta_{g}$. Then, we pose the cross-learning problem

$\displaystyle\\{\theta_{i}^{*}\\},\theta_{g}^{*}=\underset{\\{\theta_{i}\\},\theta_{g}\in\Theta}{\operatornamewithlimits{argmin}}\quad$
$\displaystyle\sum_{i=1}^{N}{\mathbb{E}}_{p_{i}(x,y)}[\ell\left(y,f(x,\theta_{i})\right)]$
(PCL) subject to
$\displaystyle\left\|\theta_{i}-\theta_{g}\right\|\leq\epsilon,\quad
i=1,\ldots,N.$

The selection of the centrality parameter $\epsilon$ entails a trade-off
between learning the specific tasks individually or sharing the data
corresponding to all tasks. Indeed, if $\epsilon$ is large enough, the
constraint is always inactive and the functions are learned separately, being
agnostic to the information provided by the data corresponding to the other
tasks. On the other extreme, if $\epsilon=0$, all policies are required to be
equal and thus the cross-learning setup reduces to consensus. In practice, the
best parameter $\epsilon$ to use is not known, and to estimate it one can
resort to cross-validation techniques [1].

Before presenting an algorithm to solve the cross-learning problem, and to
gain insights on the advantages proposed formulation (PCL), we consider the
more tractable problem of estimating the means of two Gaussian variables. This
example will allow us to show that by selecting a nonzero $\epsilon$
judiciously, the cross-learning formulation shown in (PCL) can outperform both
the agnostic and consensus estimators.

###  II-A An illustrative case of sample means

Consider the problem of estimating the means $\mu_{x}$ and $\mu_{y}$ of two
independent real-valued Gaussian random variables, from their samples
$X_{j}\sim\mathcal{N}(\mu_{x},\sigma^{2}),\ j=1,\ldots,M$ and
$Y_{j}\sim\mathcal{N}(\mu_{y},\sigma^{2}),\ j=1,\ldots,M$. If $\mu_{x}$ and
$\mu_{y}$ are estimated separately, being agnostic of the information provided
by the other sample, their maximum likelihood estimators are the sample means
$\bar{X}_{M}$ and $\bar{Y}_{M}$, and the corresponding mean squared errors
equal $\sigma^{2}_{M}:=\sigma^{2}/M$.

Now, if we know $|\mu_{x}-\mu_{y}|=\epsilon_{0}$ a priori, then we can use the
cross-learning algorithm to trade off bias for variance and improve the
estimation accuracy. Indeed, we will prove

######  Claim 1.

The Mean Squared Error (MSE) of the of the estimates obtained with the cross-
learning algorithm for $\epsilon=\epsilon_{0}$ satisfy
$\textrm{mse}_{\text{CL}}(\epsilon_{0})\leq(3/4)\sigma^{2}_{M}$.

In words, the agnostic estimator can be outperformed by a factor of $3/4$, at
least, regardless of the value of $\epsilon_{0}$ and $\sigma_{M}$. This is
true even when $M\to\infty$, where the agnostic estimator becomes consistent.
Although this result assumes exact knowledge of the ground truth
$\epsilon=\epsilon_{0}$, we will characterize the MSE
$\textrm{mse}_{\text{CL}}(\epsilon)$ as a continuous function of $\epsilon$,
and demonstrate that there is a range of values around $\epsilon_{0}$ for
which the cross-learning estimator outperforms the agnostic one. Furthermore,
we also argue that cross-learning outperforms consensus, corresponding to
$\epsilon=0$. Specifically,

######  Claim 2.

There exists $\epsilon>0$ such that
$\textrm{mse}_{\text{CL}}(\epsilon)<\textrm{mse}_{\text{CL}}(0)$, with strict
inequality for all values of $\epsilon_{0}$ and $\sigma_{M}$.

In proving these claims, we will assume without loss of generality that
$\mu_{x}=\mu_{y}+\epsilon_{0}$. Specifying the cross-learning estimator in
(PCL) for the problem at hand, we have

$\displaystyle(\hat{\mu}_{x}^{\text{CL}},\hat{\mu}_{y}^{\text{CL}})=\underset{\bar{\mu}_{x},\bar{\mu}_{y}}{\operatornamewithlimits{argmin}}\quad$
$\displaystyle\sum_{i=1}^{M}\left[\left(X_{i}-\bar{\mu}_{x}\right)^{2}+\left(Y_{i}-\bar{\mu}_{y}\right)^{2}\right]$
(2a) subject to $\displaystyle|\bar{\mu}_{x}-\bar{\mu}_{y}|\leq\epsilon$ (2b)

In this case, the problem admits the following closed-form solution,

$\displaystyle\hat{\mu}_{x}^{\text{CL}}=\begin{cases}\bar{X}_{M},&\textrm{ if
}|\bar{X}_{M}-\bar{Y}_{M}|<\epsilon\\\
(\bar{X}_{M}+\bar{Y}_{M})/2+\epsilon/2,&\textrm{ if
}\bar{X}_{M}>\bar{Y}_{M}+\epsilon\\\
(\bar{X}_{M}+\bar{Y}_{M})/2-\epsilon/2,&\textrm{ if
}\bar{X}_{M}<\bar{Y}_{M}-\epsilon\end{cases}$ (3)

where the cases correspond to the constraint being activated or not. In order
to compute the mean squared error, it is convenient to put
$\hat{\mu}_{x}^{\text{CL}}$ in terms of the variables
$\bar{Z}_{M}:=(\bar{X}_{M}+\bar{Y}_{M})/2$ and
$\bar{W}_{M}:=(\bar{X}_{M}-\bar{Y}_{M})/2$. Under these definitions, (3)
reduces to $\hat{\mu}_{x}^{\text{CL}}=\bar{Z}_{M}+T(\bar{W}_{M})$, where

$\displaystyle T(\bar{W}_{M})$
$\displaystyle:=\bar{W}_{M}\mathds{1}[|\bar{W}_{M}|<\epsilon/2]$
$\displaystyle+(\epsilon/2)\mathds{1}[\bar{W}_{M}>\epsilon/2]-(\epsilon/2)\mathds{1}[\bar{W}_{M}<-\epsilon/2].$
(4)

Since $\bar{X}_{M}$ and $\bar{Y}_{M}$ are independent, with the same variance,
then $\bar{Z}_{M}$ and $\bar{W}_{M}$ also are, thus

$\displaystyle\textrm{mse}_{\text{CL}}(\epsilon)$
$\displaystyle={\mathbb{E}}[(\hat{\mu}_{x}^{\text{CL}}-\mu_{x})^{2}]$ (5)
$\displaystyle={\mathbb{E}}[(\bar{Z}_{M}-\mu_{z})^{2}]+{\mathbb{E}}[(\bar{T}(W_{M})-\mu_{w})^{2}]$
(6)
$\displaystyle=\sigma_{M}^{2}/2+E[(\bar{W}_{M}-\epsilon_{0}/2)^{2}\mathds{1}[|\bar{W}_{M}|<\epsilon/2]$
$\displaystyle+{\mathbb{E}}[(1/4)(\epsilon-\epsilon_{0})^{2}\mathds{1}[\bar{W}_{M}>\epsilon/2]]$
$\displaystyle+{\mathbb{E}}[(1/4)(\epsilon+\epsilon_{0})^{2}\mathds{1}[\bar{W}_{M}<-\epsilon/2]].$
(7)

With $\bar{Z}_{M}$ and $\bar{W}_{M}$ being Gaussian, the expected values can
be found in terms of the error function
$\textrm{erf}(x):=(2/\sqrt{\pi})\int_{0}^{x}e^{-t^{2}}dt$,

$\displaystyle\textrm{mse}_{\text{CL}}(\epsilon)$
$\displaystyle=\frac{\sigma_{M}^{2}}{2}\left(1+\frac{\alpha
e^{-\alpha^{2}}-\beta
e^{-\beta^{2}}}{\sqrt{\pi}}+\frac{e(\beta)-\textrm{erf}(\alpha)}{2}\right)$
$\displaystyle+\frac{\sigma_{M}^{2}}{2}\left(\alpha^{2}(1+\textrm{erf}(\alpha))+\beta^{2}(1-\textrm{erf}(\beta))\right),$
(8)

with $\alpha:=-(\epsilon+\epsilon_{0})/2\sigma_{M}$ and
$\beta:=(\epsilon-\epsilon_{0})/2\sigma_{M}$.

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$$6$$0.6$$0.8$$1$$1.2$$1.4$$1.6$Centrality
($\epsilon$)Mean Squared Error (MSE)Monte CarloClosed Form Figure 1: Mean
squared error with respect to the centrality parameter $\epsilon$ of the
cross-learning problem for one variable. The illustrative case for the
estimation of the means of two Gaussian variables is shown. The two variables
have same variance $\sigma=1$ and the distance between their means is given by
$\epsilon_{0}=2$.

Figure 1 shows the mean squared error in (8) depicted in red as a function of
$\epsilon$ for a particular choice $\epsilon_{0}=2$ and $\sigma=1$. Starting
from the consensus case at $\epsilon=0$, the error reduces to a minimum value
by $\epsilon_{0}$, being $3/4$ lower than the limiting error
$\lim_{\epsilon\to\infty}\textrm{mse}_{\text{CL}}(\epsilon)=\sigma_{M}$. This
limit results from the constraint being activated with probability zero, which
yields $\hat{\mu}_{x}^{\text{CL}}=\bar{X}_{M}$, as $\epsilon$ grows unbounded.
The red curve in is double checked by comparing it to Monte Carlo estimates of
the error, over $100{,}000$ realizations of $\hat{\mu}_{x}^{\text{CL}}$,
depicted by black dots in Figure 1.

To see that Claim 1 is true regardless the values of $\sigma_{M}$ and
$\epsilon_{0}$, put $\epsilon=\epsilon_{0}$ in (8), so that $\beta=0$ and
$\alpha=-\epsilon_{0}/\sigma_{M}$, thus

$\displaystyle\textrm{mse}_{\text{CL}}(\epsilon_{0})=$
$\displaystyle\frac{\sigma_{M}^{2}}{2}\left(1+\frac{\alpha
e^{-\alpha^{2}}}{\sqrt{\pi}}-\frac{\textrm{erf}(\alpha)}{2}+\alpha^{2}(1+\textrm{erf}(\alpha))\right).$

It remains to prove that the term inside brackets in the previous expression
does not exceed $3/2$. Indeed, the function

$f(x)=1+\frac{xe^{-x^{2}}}{\sqrt{\pi}}-\frac{\textrm{erf}(x)}{2}+x^{2}(1+\textrm{erf}(x))$

attains $f(-\infty)=3/2$, $f(0)=1$, and $df/dx=2x(1+\textrm{erf}(x))<0$ in the
interval $(-\infty,0]$, which yields $f(x)\in(1,3/2),\text{for all}\
x\in(-\infty,0)$, holding in particular for
$x=-\epsilon_{0}/\sigma_{M}=\alpha$, as desired.

To see Claim 2 is also true, we compute the derivative of (8) with respect to
$\epsilon$ and show that it is negative at $\epsilon=0$,

$\displaystyle
D(\epsilon):=\frac{d}{d\epsilon}\textrm{mse}_{\text{CL}}(\epsilon)=$
$\displaystyle\sigma_{M}\left(\beta(1-\textrm{erf}(\beta))-\alpha(1+\textrm{erf}(\alpha))\right).$

For $\epsilon=0$ $\alpha=\beta=-\epsilon_{0}/2\sigma_{M}<0$, hence
$D(0)=(\epsilon_{0}/\sigma_{M})\textrm{erf}(-\epsilon_{0}/2\sigma_{M})<0$,
which is strictly negative as desired.

Algorithm 1 Cross-learning algorithm

1:Initialize $\theta_{g}=0,$ $\theta^{0}_{i}=\theta^{0},\ i=1,\ldots,N$

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

3: for $i=1,\ldots,N$ do

4: Obtain $\hat{\nabla}_{\theta_{i}}\ell(y,f(x,\theta_{i}^{k}))$.

5: end for

6:
$(\\{\theta_{i}^{k+1}\\},\theta_{g}^{k+1})=\mathcal{P}_{\mathcal{C}}\left[\\{\theta^{k}_{i}+\eta^{k}\hat{\nabla}_{\theta_{i}}\ell(y,f(x,\theta_{i}^{k})),\theta_{g}^{k}\right]$

7:until convergence

##  III Algorithm Construction

In order to solve the cross-learning problem (PCL) in a general setup, we can
resort to a stochastic projected gradient descent scheme. To take gradient
steps on the objective function, we can take derivatives of the loss function
with respect to the parameters,

$\displaystyle\bar{\theta}_{i}^{k}=\theta_{i}^{k}-\eta^{k}\hat{\nabla}_{\theta_{i}}\ell(y,f(x,\theta_{i}^{k})),$
(9)

where $\hat{\nabla}_{\theta_{i}}\ell(y,f(x,\theta_{i}^{k}))$ is a stochastic
version of the gradient associated to the loss function at a data point, and
$\eta^{k}$ possibly constant, is a learning step. After each gradient step is
taken, the restriction may not be satisfied and thus a projection must be
enforced. We define the cross-learning projection $\mathcal{P}_{\mathcal{C}}$
as,

$\displaystyle\mathcal{P}_{\mathcal{C}}\left[\\{\bar{\theta}_{i}\\},\bar{\theta}_{g}\right]=\underset{\\{\theta_{i}\\},\theta_{g}}{\operatornamewithlimits{argmin}}$
$\displaystyle\sum_{i=1}^{N}\left\|\theta_{i}-\bar{\theta}_{i}\right\|^{2}+\left\|\theta_{g}-\bar{\theta}_{g}\right\|^{2}$
(10a) s.t.
$\displaystyle\left\|\theta_{i}-\theta_{g}\right\|^{2}\leq\epsilon^{2},i=1,\ldots,N$
(10b)

By combining the gradient step given in equation (9) with the projection (10),
we obtain the cross-learning algorithm illustrated in Algorithm 1. As a
projected gradient descent form, it can be shown to converge to the optimal
value of the cross-learning problem (PCL) in the case of a convex problems
[20]. In general, the cross-learning problem (PCL) may not be convex due to
both the loss function $\ell$ and the parameterization function $f(x,\theta)$
used. However, recent results have shown that even in those cases, problems
akin to (PCL) have tractable duality gaps [21], motivating primal-dual
approaches. Nonetheless, for the specific case of the cross-learning Algorithm
1, previous results have shown its converge in high probability to a
neighborhood of a first-order stationary point of problem (PCL) in the context
of reinforcement learning [19].

Algorithm 2 Cross-learning projection

1:Initialize $\mu_{i}=\mu_{i}^{0}$ and prescribe error $\delta>0$

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

3: for $i=1,\dots,N$ do

4:
$\mu_{i}^{k+1}=[\mu_{i}^{k}+\alpha^{k}\partial_{\mu_{i}}\mathcal{L}(\theta_{i},\theta_{g},\mu^{k}_{i})]_{+}$
with eq. (14)

5: end for

6:until
$|\langle\hat{\partial}_{{\mu}}\mathcal{L}(\theta_{i},\theta_{g},\mu^{k}_{i}),{\mu}^{k}\rangle|\leq\delta$

###  III-A Projection in the dual domain

On a neural network, the function parameterization vectors $\theta_{i}$ may
have millions of parameters rendering projection (10) challenging in practice.
However, we can exploit the fact that only one constraint is added per
function, resorting to a dual domain algorithm which has only one variable per
function. With Lagrange multipliers $\mu_{i}\geq 0$, we can write the
Lagrangian $\mathcal{L}(\theta_{i},\theta_{g},\mu_{i})$ of problem (PCL) as,

$\displaystyle\mathcal{L}(\theta_{i},\theta_{g},\mu_{i})$
$\displaystyle=\sum_{i=1}^{N}\left\|\theta_{i}-\bar{\theta}_{i}\right\|^{2}+\left\|\theta_{g}-\bar{\theta}_{g}\right\|^{2}$
$\displaystyle+\sum_{i=1}^{N}\mu_{i}(\left\|\theta_{i}-\theta_{g}\right\|^{2}-\epsilon^{2}).$
(11)

Upon defining $\lambda_{i}=\frac{\mu_{i}}{1+\mu_{i}}$,
$z=[1,\lambda_{1},\dots,\lambda_{N}]$ and $a=|z|_{1}$, the primal minimizers
of the Lagrangian given $\mu$ are

$\displaystyle\theta_{g}$
$\displaystyle=\frac{1}{a}(\bar{\theta_{g}}+\sum_{i=1}^{N}\lambda_{i}\bar{\theta}_{i}),$
(12) $\displaystyle\theta_{i}$
$\displaystyle=(1-\lambda_{i})\bar{\theta}_{i}+\lambda_{i}\theta_{g}.$ (13)

Hence, we can obtain the subgradient of the Lagrangian with respect to the
Lagrangian multipliers $\mu_{i}$, by substituting (12) and (13),

$\displaystyle\partial_{\mu_{i}}\mathcal{L}(\theta_{i},$
$\displaystyle\theta_{g},\mu_{i})=\|\theta_{i}-\theta_{g}\|^{2}-\epsilon^{2}$
(14)
$\displaystyle=(1-\lambda_{i})^{2}\left\|\bar{\theta}_{i}-\frac{1}{a}(\bar{\theta_{g}}+\textstyle\sum_{i=1}^{N}\lambda_{i}\bar{\theta}_{i})\right\|^{2}-\epsilon^{2}.$

By expanding the norm in (14), the gradient of the Lagrangian depends on the
parameters $\bar{\theta}_{i}$ through their inner products
$\bar{\theta}_{i}^{T}\bar{\theta}_{i^{\prime}}$, and these products can be
computed once when initializing the primal-dual algorithm. Other than that,
Algorithm 2 operates on the reduced dimension of the dual variables.

The subgradient (14) is then embedded in the projection described in Algorithm
2. As $\mu_{i}$ must be nonnegative [22], we project them to the nonnegative
orthant at every step. The stopping condition of Algorithm 2 is the duality
gap of the projection (10), which can be made arbitrarily small. Note that the
cross-learning projection (10) is a convex problem, by selecting a non-
summable and square summable step-size [20], once the algorithm is halted, the
optimal parameters $\\{\theta_{i}\\},\theta_{g}$ of can be recovered using the
primal solutions (12)$-$(13).

##  IV Numerical Results

Art

Clipart

Product

World

(a) Alarm

(b) Bike

(c) Glasses

(d) Pen

(e) Speaker

Figure 2: Example images from $5$ of the $65$ categories from the $4$ domains
composing the Office-Home dataset [23]. The $4$ domains are Art, Clipart,
Product and Real World. In total, the dataset contains $15{,}500$ images of
different sizes and shapes.

In this section we test our cross-learning framework on a classification
problem with real data. Our goal is to classify images belonging to $P$
different categories, and the problem is divided in $N$ tasks corresponding to
images belonging to $N$ different domains. Specifically, we use the Office-
Home dataset [23]. It consists of $N=4$ different domains; Art: an artistic
representation of the object, Clipart: a clip art reproduction, Product: an
image of a product for sale, and Real World: pictures of the object captured
with a camera. The overall dataset contains $15{,}500$ RGB images divided in
$P=65$ categories, with five examples given in Figure 2, including Alarm,
Bike, Glasses, Pen, and Speaker. Notice that within each category there are
images belonging to each of the domains. The minimum number of images per
domain and category is $15$ and the image size varies from the smallest image
size of $18\times 18$ to the largest being $6500\times 4900$ pixels. We
preprocessed the images by normalizing them and fitting their size to
$224\times 224$ pixels.

We use our cross-learning strategy for classification with the intuition that
seeing an image of one domain could help identifying an image of the same
category in other domain. For instance, there is a large amount of speakers
pictured as products, which could help identify Clipart versions of the
speakers which are available in much smaller numbers. We also want to check
how the cross-learning classifier compares to using a single classifier that
merges the dataset disregarding the different domains.

We use neural networks as the classifiers $f(x,\theta_{i})$, $i=1,\ldots,N$
for cross-learning, with the architecture being based on AlexNet [24] with a
reduction on the size of the last fully connected layers to $256$ neurons per
layer, corresponding to $\theta\in{\mathbb{R}}^{S}$, with $S=4{,}911{,}745$.
In this case, it is crucial to make use of the cross-learning projection (10)
in the dual domain, thus reducing the dimensionality from $S$ to $N=4$
variables. Furthermore, we split the dataset in two parts, using $4/5$ of the
images for training and $1/5$ for testing. We train each neural network in a
domain according to Algorithm 1. For the gradient step, we use a stepsize
$\eta=0.001$ and we take one image per step. As it is standard for image
classification, we use the cross-entropy loss [1]. We train the $N=4$ neural
networks according to the cross-learning algorithm for different values of the
centrality measure $\epsilon$. We also consider the case of consensus
($\epsilon=0$) which is equivalent to merging the images from all domains and
training a single neural network. Additionally, we train the neural networks
separately one in each domain, which corresponds to $\epsilon=\infty$.

$0$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$32$$34$$36$$38$Centrality
($\epsilon$)Accuracy ($\%$) Figure 3: Accuracy in the test set of a classifier
trained with the cross-learning algorithm for different values of the
parameter $\epsilon$. Mean values and one standard deviation band are plotted,
corresponding to the accuracy over all the $65$ categories and $4$ different
domains, obtained after after $30$ epochs. The best overall accuracy is
achieved at $\epsilon=10^{-6}$. The consensus case corresponds to $\epsilon=0$
and the agnostic case corresponds to $\epsilon=\infty$ (not plotted, achieving
a lower accuracy of $24.44\%$).

$0$$20$$40$$60$$80$$100$$120$$140$$160$$180$$200$$220$$240$$186.61\%$$39.20\%$$26.98\%$$54.82\%$Improvement
($\%$)ArtClipartProductReal World (a) Percentage of improvement in the
accuracy of the classifier between the task-specific agnostic classifier
($\epsilon=\infty$) and the best cross-learning classifier.

$0$$5$$10$$15$$20$$25$$30$$35$$40$$45$$50$$26.38\%$$5.13\%$$6.65\%$$0.74\%$Improvement
($\%$) (b) Percentage of improvement in the accuracy of the classifier between
the task-specific consensus classifier ($\epsilon=0$) and the best cross-
learning classifier.

Figure 4: Percentage of improvement in the accuracy of the best cross-learning
classifier with respect to the agnostic, and consensus classifiers.

In Figure 3, we present the result of these experiments using the
classification accuracy of the trained classifier as figure of merit. We
corroborate the intuitive idea drawn by looking at the images (cf. Figure 2).
Namely, that domains are correlated and adding samples from different domains
improves the overall performance of the agnostic policy, only trained on its
own specific samples. This is exemplified by consensus ($\epsilon=0$)
outperforming the agnostic training ($\epsilon=\infty$), with an accuracy of
$33.97\%$ against $24.44\%$. More importantly, the appropriate choice of the
$\epsilon$ parameter in the cross-learning method outperforms both consensus
and agnostic approaches. In particular, Figure 3 shows that the maximum
performance is achieved at $\epsilon=10^{-6}$, indicating that keeping the
parameters of the neural networks close, yet not merging them, induces an
improvement in the overall accuracy. In this sense, these experimental results
with real data recover our theoretical findings about the behavior of the
cross-learning estimator in the Gaussian example of Section II-A.

The advantage of cross-learning can be better seen in Figure 4, as in all $4$
domains there exists a value of $\epsilon$ that outperforms both the consensus
and agnostic counterparts across all of them. In particular, for the domain
that has the minimum number of samples (Art, with $2{,}427$), the accuracy
almost triples with cross-learning when compared with the agnostic classifier.
This should not come as a surprise, as due to its limited number of samples,
it admits a larger margin of improvement by including samples from additional
domains. Likewise, considerable improvements are still found when compared
with the consensus classifier.

##  V Conclusion

In this paper, we introduced a cross-learning framework for multi-task
supervised learning. The goal is to learn a set of functions, either for
classification or regression, corresponding to different tasks. The proposed
strategy entails the fitting of the parameters of each function to their task-
specific data, while keeping the parameters of all functions close to each
other. We derived intuition from an illustrative case with Gaussian
distributions, corroborating the advantage of our formulation compared both to
the consensus or the separate approach. For generic set tasks, we have derived
a stochastic projected gradient algorithm, together with a dual implementation
of the projection for reducing the problem dimension to the number of tasks.
Numerical experiments on a dataset with different domains validate that the
cross-learning approach is capable of outperforming both the domain specific
and the consensus classifiers.

## References

  * [1] T. Hastie, R. Tibshirani, and J. Friedman, _The elements of statistical learning: data mining, inference, and prediction_. Springer Science & Business Media, 2009.
  * [2] R. Caruana, “Multitask learning,” _Machine learning_ , vol. 28, no. 1, pp. 41–75, 1997.
  * [3] C. Stein, “Inadmissibility of the usual estimator for the mean of amultivariate normal distribution,” in _Proceeding of the fourth Berkeley symposium on mathematical statistics and probability_ , vol. 1. University of California Press, 1956, pp. 197–206.
  * [4] S. Liu, E. Johns, and A. J. Davison, “End-to-end multi-task learning with attention,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 1871–1880.
  * [5] I. Misra, A. Shrivastava, A. Gupta, and M. Hebert, “Cross-stitch networks for multi-task learning,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 3994–4003.
  * [6] A. Kendall, Y. Gal, and R. Cipolla, “Multi-task learning using uncertainty to weigh losses for scene geometry and semantics,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 7482–7491.
  * [7] P. Liu, X. Qiu, and X.-J. Huang, “Adversarial multi-task learning for text classification,” in _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics_ , 2017, pp. 1–10.
  * [8] D. Dong, H. Wu, W. He, D. Yu, and H. Wang, “Multi-task learning for multiple language translation,” in _Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing_ , 2015, pp. 1723–1732.
  * [9] Y. Zhang and Q. Yang, “A survey on multi-task learning,” _arXiv preprint arXiv:1707.08114_ , 2017.
  * [10] S. Ben-David and R. Schuller, “Exploiting task relatedness for multiple task learning,” in _Learning Theory and Kernel Machines_. Springer, 2003, pp. 567–580.
  * [11] A. Argyriou, T. Evgeniou, and M. Pontil, “Convex multi-task feature learning,” _Machine learning_ , vol. 73, no. 3, pp. 243–272, 2008.
  * [12] T. Evgeniou and M. Pontil, “Regularized multi–task learning,” in _Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining_ , 2004, pp. 109–117.
  * [13] A. Agarwal, S. Gerber, and H. Daume, “Learning multiple tasks using manifold regularization,” in _Advances in neural information processing systems_ , 2010, pp. 46–54.
  * [14] L. Jacob, J.-p. Vert, and F. R. Bach, “Clustered multi-task learning: A convex formulation,” in _Advances in neural information processing systems_ , 2009, pp. 745–752.
  * [15] T. Kato, H. Kashima, M. Sugiyama, and K. Asai, “Multi-task learning via conic programming,” in _Advances in Neural Information Processing Systems_ , 2008, pp. 737–744.
  * [16] ——, “Conic programming for multitask learning,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 22, no. 7, pp. 957–968, 2009.
  * [17] O. Sener and V. Koltun, “Multi-task learning as multi-objective optimization,” in _Advances in Neural Information Processing Systems_ , vol. 31. Curran Associates, Inc., 2018\.
  * [18] J. Cerviño, J. A. Bazerque, M. Calvo-Fullana, and A. Ribeiro, “Meta-learning through coupled optimization in reproducing kernel hilbert spaces,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 4840–4846.
  * [19] ——, “Multi-task reinforcement learning in reproducing kernel hilbert spaces via cross-learning,” _arXiv preprint arXiv:2008.11895_ , 2020.
  * [20] D. P. Bertsekas, _Nonlinear Programming_. Athena Scientific, 1999.
  * [21] L. F. Chamon, S. Paternain, M. Calvo-Fullana, and A. Ribeiro, “The empirical duality gap of constrained statistical learning,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 8374–8378.
  * [22] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge University Press, 2009.
  * [23] H. Venkateswara, J. Eusebio, S. Chakraborty, and S. Panchanathan, “Deep hashing network for unsupervised domain adaptation,” in _(IEEE) Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017.
  * [24] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.