#  Anisotropic van der Waals Dispersion Forces in Polymers:
Structural Symmetry Breaking Leads to Enhanced Conformational Search

Mario Galante  Alexandre Tkatchenko  Department of Physics and Materials
Science, University of Luxembourg, L-1511 Luxembourg City, Luxembourg

###### Abstract

The modeling of conformations and dynamics of (bio)polymers is of primary
importance for understanding physicochemical properties of soft matter.
Although short-range interactions such as covalent and hydrogen bonding
control the local arrangement of polymers, non-covalent interactions play a
dominant role in determining the global conformations. Here we focus on how
the inclusion of many-body effects in van der Waals (vdW) dispersion influence
the molecular dynamics of model polymers and a poly-alanine chain studied with
semiempirical quantum mechanics. We show that delocalized force contributions
are key to explore the conformational landscape, as they induce an anisotropic
polarization response which efficiently guides the conformation towards
globally optimized structures. This is in contrast with the commonly used
pair-wise interatomic Lennard-Jones-like potentials, which result in
structures that optimize local regularity rather than global symmetries. The
many-body vdW formalism results in a significant reduction in the roughness of
the energy landscape and drives the conformational search towards geometries
that are consistent with global spatial symmetries.

††preprint: APS/123-QED

The prediction of structure formation in polymers under realistic
thermodynamical conditions is a challenging problem that is relevant for
fields ranging from medicine and biology to physics and engineering. For
example, the folding of proteins into their secondary structure underpins the
regulation of many functions of cellular activity Bu and Callaway (2011). In
addition, the ubiquitous presence of polymeric materials in modern
applications is due to their variety of macroscopic properties which are in
turn dictated by polymer morphology and its dynamics Adhikari and Majumdar
(2004); Dersch  _et al._ (2005); Ulery  _et al._ (2011). However,
conformational stabilities are easily undermined by thermal fluctuations and
the solvent environment, therefore the knowledge of significant portions of
the potential energy surface (PES) of the system become essential Napolitano
(2015); Fraser  _et al._ (2009); Uversky (2009). Models based on
phenomenological definitions of inter-fragment interactions Souza and et al
(2021) have consistently proved to be invaluable tools for accurate
predictions of dynamical properties in biological systems, as they allow to
probe time and length scales inaccessible to atomistic simulations.
Nevertheless, modern accurate electronic structure calculations have become
feasible for systems of biological relevance. The analysis of electron
correlations at such scales can therefore shed new light on mechanisms driving
dynamical processes, and therefore be inspiration for the further development
of accurate force fields Unke and et al (2021); Unke  _et al._ (2021).

The determination of an appropriate first-principles description for the
molecular interactions that are relevant for large and relatively flexible
systems has thus become a timely and challenging issue. The general physical
model Tran and Ramprasad (2020); Abeln  _et al._ (2019) includes bonded
interactions together with non-covalent potentials to account for
electrostatics and van der Waals (vdW) dispersion interactions. From a
physical perspective, these are the result of the interaction between
instantaneous charge fluctuations Hermann  _et al._ (2017); Stöhr and
Tkatchenko (2019) and play a dominant role in determining the global
arrangement because of their strikingly slow spatial decay. Traditionally,
they are included by means of the pair-wise (PW) approximation, where the
energy decays strictly with the $6$th power of the interatomic distance and
its strength is dictated by the local chemical environments Tkatchenko and
Scheffler (2009); Sato and Nakai (2009); Grimme  _et al._ (2010). However,
recent studies Dobson (2014); Ambrosetti  _et al._ (2016); Stöhr and
Tkatchenko (2019); Hauseux  _et al._ (2020) have shown that the poor
versatility that derives from the lack of many-body effects makes the PW
approximation not suitable to adequately capture the complexity of vdW
interactions in a variety of systems. However, the implications of noncovalent
many-body interactions on molecular dynamics have not been studied.

In this Letter we focus on the qualitative features of the conformational
search for PW and many-body dispersion (MBD) potentials by means of molecular
dynamics (MD) simulations for a model carbon-based polymer. We show that the
many-body character of vdW forces enables one to capture anisotropic
polarizability effects, which consistently drive the molecular structure to
optimize the global symmetry of the conformation according to geometrical
constraints. On the contrary, the shorter range of the forces within the PW
approximations, together with their independence on the global atomic
configuration, yields conformations that are less compact and an increased
roughness of the PES. We then present the results of MD simulations for an
Ala_15 polypeptide where the electronic interactions are treated with semi-
empirical quantum mechanics, which confirms the trends observed for the model
polymer system.

We begin by considering a toy model carbon chain represented by the model
Hamiltonian

$\displaystyle\mathcal{H}$
$\displaystyle=\sum_{i=1}^{N-1}k\Big{(}\big{|}\big{|}\bm{R}_{i}-\bm{R}_{i+1}\big{|}\big{|}-R_{0}\Big{)}^{2}$
$\displaystyle+\sum_{i=1}^{N}\sum_{j\neq i,i\pm 1}a\
\text{exp}\Big{\\{}-b\big{|}\big{|}\bm{R}_{i}-\bm{R}_{j}\big{|}\big{|}\Big{\\}}+E_{\rm{vdW}},$

with $\bm{R}_{i}$ representing the position of the $i$th carbon atom,
$R_{0}=$$1.52\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, $k=39.57$ a.u. ,
$a=1878.38$ a.u. and $b=5.49$ a.u. . The first term represents an harmonic
covalent bond between the nearest-neighbours, while the second summation
introduces a Born-Mayer repulsive potential between non-nearest-neighbours
that prevents artificial agglomeration of the polymer. We note that the
parameters were chosen in order to ensure a solid backbone, however they do
not aim to faithfully represent the chemical environment. In a similar spirit,
we choose here to entirely neglect bond angle and dihedral constraints in
order to further reduce the number of degrees of freedom and favours a
qualitative comparison of the PES resulting from different formulations for
the vdW energy, $E_{\rm{vdW}}$. These assumptions are therefore instrumental
to detect the signatures characteristic of MBD in a system of reduced
complexity, which then enables targeted analysis in realistic simulation
conditions. Dispersion interactions can be generally quantified by means of
the adiabatic connection fluctuation-dissipation (ACFD) formula Hermann  _et
al._ (2017); Tkatchenko  _et al._ (2013),

$\displaystyle E_{\text{vdW}}^{\textrm{MBD}}$
$\displaystyle=\frac{\hbar}{2\pi}\int_{0}^{\infty}du\int_{0}^{1}d\lambda\int\int
d\bm{r}d\bm{r}^{\prime}$
$\displaystyle\times\text{tr}\Big{[}\big{(}\bm{\alpha}_{\lambda}(\bm{r},\bm{r}^{\prime},iu)-\bm{\alpha}_{0}(\bm{r},\bm{r}^{\prime},iu)\big{)}\bm{T}(\bm{r},\bm{r}^{\prime})\Big{]},$

where $\bm{r}$, $\bm{r}^{\prime}$ are spatial variables, $u$ is the frequency
of the electric field, tr denotes the trace of the polarizability tensors,
$\bm{\alpha}$, and the integral over the coupling strength, $\lambda$,
expresses the adiabatic switching-on of the electron-electron interaction.
Here we approximate the ACFD equation within the many-body dispersion approach
Tkatchenko  _et al._ (2012, 2013) and we define the coupling tensor, $\bm{T}$,
as the dipole-dipole potential corresponding to two Gaussian charge
distributions. The polarizability tensor of the system corresponding to
uncorrelated atoms, $\bm{\alpha}_{0}$, can be written as the sum of the
polarizabilities localized at each atomic site,

$\bm{\alpha}_{\text{0}}(\bm{r},\bm{r}^{\prime})\approx\sum_{i}\delta^{3}(\bm{r}-\bm{R}_{i})\delta^{3}(\bm{r}-\bm{r}^{\prime})\mathbb{1}_{3}\alpha_{i}^{\text{TS}},$
(3)

where $\alpha_{i}^{\rm{TS}}$ accounts for all short range correlations and is
here taken as the isotropic atomic polarizability for a generic carbon
calculated within the Tkatchenko-Scheffler (TS) scheme Tkatchenko and
Scheffler (2009) , $\alpha^{\text{TS}}_{\text{C}}=10.44$ a.u. The
polarizability of the interacting system, $\bm{\alpha}_{\lambda}$, can then be
calculated through the Dyson-like screening equation,

$\bm{\alpha}_{\lambda}=\sum_{n=0}^{\infty}\langle\bm{\alpha}_{\text{0}}(-\lambda\bm{T}\bm{\alpha}_{\text{0}})^{n}\rangle,$
(4)

which allows to construct a non-diagonal tensor starting from the local atomic
polarizabilites. We remark that the solution of Eqs. (Anisotropic van der
Waals Dispersion Forces in Polymers: Structural Symmetry Breaking Leads to
Enhanced Conformational Search) and (4) within these approximations
corresponds to solving the Schrödinger equation for a system of quantum
harmonic oscillators centered at the atomic sites and coupled via the dipole-
dipole potential Hermann  _et al._ (2017). Therefore, the MBD formalism
accounts for the energy of global charge oscillations expressed through linear
combinations of the atomic dipolar displacement vectors. In other terms, the
MBD energy is the sum of the interaction energies of the $3N$ normal modes,
with $N$ being the total number of oscillators. In contrast, the PW
approximation corresponds to truncating the series of Eq. (4) to the second
order and considers only unscreened dipole-dipole interactions,

$E^{\text{PW}}_{\text{vdW}}=-\frac{1}{2}\sum_{i\neq
j}f_{\text{damp}}(R_{ij})\frac{C_{6,ij}}{R_{ij}^{6}},$ (5)

where $f_{\text{damp}}$ is a Fermi-type damping function that eliminates short
range contributions and $C_{6,ij}=\frac{3}{\pi}\hbar\int
du\alpha^{\text{TS}}_{i}(iu)\alpha^{\text{TS}}_{j}(iu)$ .

Figure 1:  (Color Online) Dependence of the averaged pair-wise (a) and many-
body (b) force components for the interatomic distance for a conformation
obtained through molecular dynamics for a $100$ atoms chain (shown in the
inset).

We begin by analyzing the differences between the long range behavior of the
components of the PW and MBD forces for a representative conformation of a 100
atoms chain extracted from a MD trajectory [Fig. 1]. For simplicity, we
constrain the atoms at each edge of the chain to oscillate around their
equilibrium position through a quadratic potential. The PW force between a
given pair of atoms can be calculated analytically from Eq. 5 and therefore
follows by construction the $R^{-7}$ power law, becoming negligible within
$15\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. On the contrary, the MBD force is
the result of the competition between all normal oscillation modes of the
coupled dipoles. Here we numerically calculate the contribution acting on the
$i$-th site coming from the $j$-th atom by using the formula,

$F_{ij,\alpha}^{\textrm{MBD}}=-\frac{\partial}{\partial
R_{i}^{\alpha}}\sum_{m=1}^{3N}\left[\varepsilon_{m}\chi_{j}^{m}\right],\quad\chi_{j}^{m}=\sum_{\alpha=x,y,z}\big{(}\chi_{j,\alpha}^{m}\big{)}^{2}$
(6)

where $\varepsilon_{m}$ is the eigenvalue of the $m$-th MBD mode, and
$\chi_{j,\alpha}^{m}$ is the $\alpha=x,y,z$ component of the displacement
vector of the $j$-th dipole within the same mode. From panel (b) it is evident
that the average contribution to MBD forces is entropically scattered
throughout the whole distance range and no clear decay rate can be singled
out. This can be attributed to the elevated degree of delocalization of the
MBD eigenmodes, as shown in Ref. Ambrosetti  _et al._ (2016). Moreover, the
persistence of MBD interactions for long distances has already been reported
in model and realistic systems Hauseux  _et al._ (2022); Stöhr and Tkatchenko
(2019), and their structure dependence will be analyzed in future
publications. Furthermore, there is no evident relation between the MBD force
components and the (cylindrical-like) spatial geometry of the conformation.
This is in contrast with the PW case, where the longitudinal ($x$) component
is more pronounced than the others. We also remark that the typical
contribution at a given distance is three orders of magnitude higher for MBD
than for PW forces, although in both cases contributions of opposite sign
balance to yield an average atomic force of
$\approx$0.01\text{\,}\mathrm{eV}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}$$.
This can be associated to the interplay between the contributions due to the
$3N$ normal modes, consistently with the analysis of the atomic force response
presented in Ref. Haseux2021.

Figure 2: (Color Online) Evolution of the average atomic polarizability,
$\alpha$, (top) and of the vdW energy (bottom) along an MBD optimization path.
Panel (a): the orange solid line represents the input atomic polarizability,
$\alpha_{C}^{\text{TS}}$, while the dotted, dashed and dot-dashed lines
correspond to the $x$, $y$ and $z$ components of the average atomic
polarizability, respectively. Panel (b): the dashed blue (solid orange) curve
represents the trend of the MBD (PW) energy with the progression of the
optimization. The black squares denote the vdW energy and optimization step of
the remaining three conformations displayed.

We further analyze this by evaluating the polarizability tensor by means of
the Dyson-like equation of Eq. (4) for $\lambda=1$ (its trace corresponding to
experimentally measured polarizability). In Fig. 2 (a) we display the average
atomic polarizability along the Cartesian directions,

$\alpha_{\text{tot}}^{a}=N^{-1}\sum_{ij}\sum_{b}\alpha_{\lambda=1}^{ab}(\bm{R}_{i},\bm{R}_{j}),$
(7)

calculated for a series of configurations for a $40$ atoms chain obtained
through the Newton-Raphson optimization with respect to the MBD energy of an
input two dimensional conformation, while maintaining fixed the coordinates of
the edge atoms. It is clear that the screened atomic total polarizabilty
remains mostly constant, similarly to the PW atomic polarizability that is
fixed by definition to $\alpha^{\textrm{TS}}_{\textrm{C}}$. However, the
components of the total screened polarizability gradually change to reproduce
the the cylindrical-like symmetry imposed by the constrain on the edge atoms,
with the longitudinal ($x$) component growing dominant with respect to the
others. We remark that this is analogous to the asymmetry of the force
components observed in Fig. 1, and it is evident that the anisotropic
character of MBD polarizabilities is responsible for driving structures
towards three-dimensional conformations. In fact, the change in the symmetry
of the total polarizability goes on equal footing with the lowering of the MBD
energy, as shown in Fig. 2 (b). For comparison, we display the PW energy
calculated for the same structures. We find that the energy gain along the
sequence is significantly less pronounced for PW than for MBD, with the trend
of the former displaying two local minima. This behaviour hints to the fact
that MBD forces may favour a more efficient conformational optimization than
their PW counterpart. We note that the increase of the efficiency of
conformational seaches, measured in terms of the decrease of energy landscape
roughness, has already been linked with the increase of interaction range in
the context of pair-wise potentials Hoare and McInnes (1976, 1983); Doye and
Wales (1996); Rechtsman  _et al._ (2005, 2006).

Figure 3: (Color Online) Results of MD simulations (100 per vdW model) at
$5\text{\,}\mathrm{K}$ for a $40$ atoms chain starting from the 2D
conformation defined by $n=5$ and $f=2$. Orange (blue) lines represent the
evolution of the total energy with time for the PW (MBD) approximation of the
vdW energy. The structures on the right represent the steady-state PW
conformations with the highest (top) and the lowest (bottom) energy, while the
structure on the left is a typical MBD steady-state structure. The inset shows
the Gaussian distributions centered on steady-state energies for a each vdW
model, rescaled by $E_{\text{min}}=\text{min}\\{E_{i}\\}$ and normalized so
that they peak at $1$.

To confirm our hypothesis, we performed $100$ MD simulations for each vdW
model, all run at $5\text{\,}\mathrm{K}$ and starting from the same 2D
conformation of the $40$ atom polymer displayed in Fig. 2. In Fig. 3 we
display the time evolution up to $100\text{\,}\mathrm{ps}$ of the total energy
for all MD trajectories, with the inset showing the Gaussian-broadened
distribution of the steady-state energies. We exclude from the latter the 5
highest energies as they are not statistically representative. It is clear
that PW trajectories yield a large variety of steady-state energies, with the
lowest distinctly separated from the others, while MBD energies are more
condensed towards the lowest minimum, $E_{\text{min}}$. To quantify this
difference we calculated the roughness of the potential energy landscape,
defined as the root-mean-squared energy distance from $E_{\text{min}}$,
$\langle(E-\text{min}\\{E\\})^{2}\rangle^{1/2}$ Milanesi  _et al._ (2012),
which gives $\sim 137\ k_{\text{B}}T$ and $\sim 63\ k_{\text{B}}T$ for PW and
MBD, respectively. Moreover, 19 of the MBD trajectories undergo a transition
after the initial transient ($t\gtrsim$$20\text{\,}\mathrm{ps}$), while this
occurs only in 10 cases for PW. Although the height of the energy barriers is
roughly similar ($\approx$$10\text{\,}\mathrm{meV}$), MBD transitions have
average longer duration ($\approx$$34\text{\,}\mathrm{ps}$ against
$23\text{\,}\mathrm{ps}$) and result in a larger energy decrease
($\approx$$0.21\text{\,}\mathrm{eV}$ against $0.15\text{\,}\mathrm{eV}$).
Therefore, we can conclude that the introduction of anisotropic
polarizabilities in the model polymer results in a two-fold decrease of PES
roughness and enchances the likelyhood of structural transitions towards lower
energy local minima.

Figure 4: (Color Online) Top: Structure formation for poly-alanine (n=15).
Time dependence of the anisotropy of the gyration ratio normalized to the
initial geometries for a model Carbon chain (b) and for poly-alanine (c). We
note that the absence of angular contributions in the covalent bond of the
model implies a stronger stability to thermal fluctuations. On the contrary,
the poly-alanine systems are simulated whithin a full treatment of the
electronic structure, which therefore strongly compete with vdW interactions
creating the oscillations displayed in panel (c).

As further validation, we present the results of
$300\text{\,}\mathrm{p}\mathrm{s}$ simulations at $300\text{\,}\mathrm{K}$ of
a linear poly-alanine chain with 15 residues in gas phase. At each time-step
the electronic structure is treated within the DFTB+vdW method, which includes
covalent bonding, electrostatics/polarization, Pauli repulsion and vdW
interactions Hourahine  _et al._ (2020). Despite the different simulation
conditions, we find that the dynamical structure formation for each vdW
approximation displays analogous features across the two systems. For the
model polymer, we find that the conformations corresponding to the lowest
steady-state energies for PW has a regular helical shape. However, the latter
are asymmetric with respect to all three cartesian coordinates, which is not
consistent with the cylindrical symmetry imposed by the constraint on the
positions of the edge atoms. On the contrary, MBD generally yields knot-like
structures where the local regularity is sacrificed in favour of a
cylindrical-like symmetry. This originates in the versatility of MBD in
addressing the anisotropy of the polarization response, which guides the
conformational search towards geometries that are consistent with the
geometrical constraints imposed by the environment. Similar features are
encountered in the case of the folding of the poly-alanine chain: Figure 4 (a)
displays the conformation with the lowest steady-state energy (out of 5 runs
per model) for PW and MBD. Here no constraint on the chain edges was applied,
and in both cases a local helical structure is obtained. Nevertheless, in the
PW case the chain maintains a globally linear structure while the MBD one
folds on itself, optimizing the geometrical symmetry and resulting in a total
energy gain of $0.285\text{\,}\mathrm{eV}$. We remark that more compact,
rather than seemingly linear, structures are consistent with previous studies
on conformational searches lead by global non-specific interactions Osmanovic
and Rabin (2016); Chakrabartty and Baldwin (1995); Huyghues-Despointes  _et
al._ (1995); Sung (2017). The obtained poly-alanine structures in the gas
phase cannot be compared to protein structures formed under physiological
conditions, which would require more sophisticated simulations in presence of
solvent and ions. Nevertheless, a force field trained on DFTB+MBD trajectories
was recently applied to poly-alanine and crambin in solution with results
consistent with experimental data and showing significant discrepancies with
classical force fields Unke and et al (2022).

We conclude by comparing the conformational stabilities at the steady-state
for the lowest steady-state energy trajectories for each model and for the two
systems. Panels (b) and (c) of Fig. 4 show the time-dependence of the
anisotropy of the gyration ratio, normalized at the value for the initial
structure, which measures the time evolution of the structural compactness. On
the one hand, the two vdW approximations yield distinctively different
conformational anisotropy for the model polymer. On the other hand, the
conformational anisotropy of the poly-alanine has similar absolute minima,
however the PW model yields oscillations with larger amplitude and lower
frequency. This indicates that the MBD conformations here observed are
systematically more stable than the PW ones, in consistency with the higher
entropic stability of MBD-optimized structures already demonstrated in the
context of molecular crystals Reilly and Tkatchenko (2014). The more
pronounced instability of the poly-alanine systems can attributed to the
competition between dispersion interactions and the angular component of the
covalent bonding, that is not included in the model polymer. Nevertheless, we
note that the steady-state time-average of the anisotropy of the
polarizability tensor in poly-alanine was calculated to be $\approx 1.5$ times
higher for PW rather than MBD, in full consistency with the discussion around
Fig. 2. In conclusion, we have shown that the inclusion of many-body effects
in van der Waals dispersion interactions allows to account for anisotropic
polarization responses, which depend on the global conformation. Such effect
yields an increased accuracy with respect to the PW approximation, however our
results do not invalidate the numerous findings obtained neglecting many-body
effects. Nevertheless, our analysis shows that anisotropic dispersion
interactions can be key to achieve conformations consistent with the
geometrical constraints imposed by the environment, and future efforts will be
dedicated to solvated proteins in connection to physiological conditions.

###### Acknowledgements.

The authors acknowledge support from the ERC Consolidator Grant “BeStMo”. MG
is grateful to Marco Pezzutto and Matteo Gori for many fruitful discussions.
The calculations on the poly-alanine systems were carried out with the HPC
facilities at the University of Luxembourg (see hpc.uni.lu).

## References

  * Bu and Callaway (2011) Z. Bu and D. J. Callaway, in  _Protein Structure and Diseases_ , Advances in Protein Chemistry and Structural Biology, Vol. 83, edited by R. Donev (Academic Press, 2011) pp. 163–221.
  * Adhikari and Majumdar (2004) B. Adhikari and S. Majumdar, Prog. Polym. Sci. 29, 699 (2004).
  * Dersch  _et al._ (2005) R. Dersch, M. Steinhart, U. Boudriot, A. Greiner, and J. Wendorff, Polym. Adv. Technol. 6, 276 (2005).
  * Ulery  _et al._ (2011) B. Ulery, L. Nair, and C. Laurencin, Polym. Adv. Technol. 49, 832 (2011).
  * Napolitano (2015) S. Napolitano,  _Non-Equilibrium Phenomena in Confined Soft Matter_ (Springer, 2015).
  * Fraser  _et al._ (2009) J. S. Fraser, M. W. Clarkson, S. C. Degnan, R. Erion, D. Kern, and T. Alber, Nature 462, 669 (2009).
  * Uversky (2009) V. N. Uversky, in  _Protein Folding and Misfolding: Neurodegenerative Diseases_ , edited by J. Ovádi and F. Orosz (Springer, 2009) pp. 21–75.
  * Souza and et al (2021) P. C. T. Souza and et al, Nature Methods 382 (2021).
  * Unke and et al (2021) O. T. Unke and et al, Nat. Comm. 12, 7273 (2021).
  * Unke  _et al._ (2021) O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko, and K.-R. Müller, Chemical Reviews 121, 10142 (2021), https://doi.org/10.1021/acs.chemrev.0c01111 .
  * Tran and Ramprasad (2020) H. D. Tran and R. Ramprasad, J. Phys. Chem. Lett. 11, 5823 (2020).
  * Abeln  _et al._ (2019) S. Abeln, K. A. Feenstra, and J. Heringa, in  _Encyclopedia of Bioinformatics and Computational Biology_ , Vol. 2, edited by S. Ranganathan (Elsevier, 2019) p. 497.
  * Hermann  _et al._ (2017) J. Hermann, R. A. DiStasio Jr., and A. Tkatchenko, Chem. Rev. 117, 4714 (2017).
  * Stöhr and Tkatchenko (2019) M. Stöhr and A. Tkatchenko, Sci. Adv. 5, eaax0024 (2019).
  * Tkatchenko and Scheffler (2009) A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 73005 (2009).
  * Sato and Nakai (2009) T. Sato and H. Nakai, J. Chem. Phys. 131, 224104 (2009).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 138, 154104 (2010).
  * Dobson (2014) J. F. Dobson, Int. J. Quantum Chem. 114, 1157 (2014).
  * Ambrosetti  _et al._ (2016) A. Ambrosetti, N. Ferri, R. A. DiStasio Jr., and A. Tkatchenko, Science 351, 1171 (2016).
  * Hauseux  _et al._ (2020) P. Hauseux, T. Nguyen, A. Ambrosetti, K. S. Ruiz, S. P. A. Bordas, and A. Tkatchenko, Nat. Commun. 11, 1651 (2020).
  * Tkatchenko  _et al._ (2013) A. Tkatchenko, A. Ambrosetti, and R. A. DiStasio Jr., J. Chem. Phys. 138, 074106 (2013).
  * Tkatchenko  _et al._ (2012) A. Tkatchenko, R. A. DiStasio Jr., R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012).
  * Hauseux  _et al._ (2022) P. Hauseux, A. Ambrosetti, S. P. A. Bordas, and A. Tkatchenko, Phys. Rev. Lett. 128, 106101 (2022).
  * Hoare and McInnes (1976) M. R. Hoare and J. McInnes, J. Chem. Soc. 61, 12 (1976).
  * Hoare and McInnes (1983) M. R. Hoare and J. McInnes, Adv. Phys. 32, 791 (1983).
  * Doye and Wales (1996) J. P. K. Doye and D. J. Wales, J. Chem. Phys. 105, 84228 (1996).
  * Rechtsman  _et al._ (2005) M. C. Rechtsman, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 95, 228301 (2005).
  * Rechtsman  _et al._ (2006) M. C. Rechtsman, F. H. Stillinger, and S. Torquato, Phys. Rev. E 73, 011406 (2006).
  * Milanesi  _et al._ (2012) L. Milanesi, J. P. Waltho, C. A. Hunter, D. J. Shaw, G. S. Beddard, G. D. Reid, S. Dev, and M. Volk, PNAS 109, 19563 (2012).
  * Hourahine  _et al._ (2020) B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Y. Deshaye, T. Dumitrică, A. Dominguez, S. Ehlert, M. Elstner, T. van der Heide, J. Hermann, S. Irle, J. J. Kranz, C. Köhler, T. Kowalczyk, T. Kubař, I. S. Lee, V. Lutsker, R. J. Maurer, S. K. Min, I. Mitchell, C. Negre, T. A. Niehaus, A. M. N. Niklasson, A. J. Page, A. Pecchia, G. Penazzi, M. P. Persson, J. Řezáč, C. G. Sánchez, M. Sternberg, M. Stöhr, F. Stuckenberg, A. Tkatchenko, V. W.-z. Yu, and T. Frauenheim, The Journal of Chemical Physics 152, 124101 (2020).
  * Osmanovic and Rabin (2016) D. Osmanovic and Y. Rabin, J. Chem. Phys. 144, 205104 (2016).
  * Chakrabartty and Baldwin (1995) A. Chakrabartty and R. L. Baldwin, Adv. Prot. Chem. 46, 141 (1995).
  * Huyghues-Despointes  _et al._ (1995) B. M. P. Huyghues-Despointes, T. M. Klinger, and R. L. Baldwin, Biochemistry 34, 13267 (1995).
  * Sung (2017) S.-S. Sung, Prot. Sci. 26, 2003 (2017).
  * Unke and et al (2022) O. T. Unke and et al, (2022), preprint at https://arxiv.org/pdf/2205.08306.pdf.
  * Reilly and Tkatchenko (2014) A. M. Reilly and A. Tkatchenko, Phys. Rev. Lett. 113, 0655701 (2014).


# Sampling-Frequency-Independent Audio Source Separation Using Convolution
Layer Based on Impulse Invariant Method ††thanks: This work was supported by
JSPS KAKENHI Grant Number JP20K19818.

Koichi Saito†, Tomohiko Nakamura†, Kohei Yatabe‡, Yuma Koizumi⋆, Hiroshi
Saruwatari†
†Graduate School of Information Science and Technology, The University of
Tokyo, Tokyo, Japan
‡Department of Intermedia Art and Science, Waseda University, Tokyo, Japan
⋆NTT Corporation, Tokyo, Japan

###### Abstract

Audio source separation is often used as preprocessing of various
applications, and one of its ultimate goals is to construct a single versatile
model capable of dealing with the varieties of audio signals. Since sampling
frequency, one of the audio signal varieties, is usually application specific,
the preceding audio source separation model should be able to deal with audio
signals of all sampling frequencies specified in the target applications.
However, conventional models based on deep neural networks (DNNs) are trained
only at the sampling frequency specified by the training data, and there are
no guarantees that they work with unseen sampling frequencies. In this paper,
we propose a convolution layer capable of handling arbitrary sampling
frequencies by a single DNN. Through music source separation experiments, we
show that the introduction of the proposed layer enables a conventional audio
source separation model to consistently work with even unseen sampling
frequencies.

###### Index Terms:

Audio source separation, analog-to-digital filter conversion, deep neural
networks

##  I Introduction

Figure 1: Architectures of (a) Conv-TasNet [3] and (b) proposed model, and (c)
illustration of proposed SFI convolution layer.

Audio source separation is a technique for extracting individual sources from
a mixture signal. It is one of the fundamental techniques for various audio
applications including music remixing, automatic music transcription, and
automatic speech recognition. The recent development of source separation has
been built upon machine leaning techniques using the deep neural network (DNN)
[1, 2, 3, 4, 5, 6, 7, 8, 9]. Since source separation is often utilized as a
preprocess of another task, one of the ultimate developmental goals is to
construct a single universal DNN that can be utilized as the preprocessor for
any application. To realize such an almighty source separator, every variety
of applications and conditions must be handled by a single DNN.

One important but often unnoticeable variety of audio signals is sampling
frequency. It is usually application specific, and hence a preprocessor must
be designed for the sampling frequency specified by the following application.
For example, for music remixing and editing, $44.1$ and $48$ kHz are usually
used as sampling frequencies to cover the entire human audible range [5, 8].
This is because these applications are aimed at human listeners. In contrast,
the applications aimed at the recognition of the contents contained in audio
signals do not require such full-band information. For example, beat tracking
may use $16$ kHz [10], automatic music transcription may use $11.025$ and
$22.05$ kHz [11, 12], and automatic speech recognition may use $8$ and $16$
kHz [13, 14, 15]. A versatile preprocessor must be able to handle signals
sampled at all of these sampling frequencies.

However, ordinary DNNs cannot handle audio signals sampled at various sampling
frequencies. Conventional DNN-based models work well for the sampling
frequency specified by the training data [1, 2, 3, 4, 5, 6, 7, 8, 9]. The
parameters of a DNN are trained to adapt for the training dataset, and thus
there is no guarantee of applicability for signals that are sampled at the
other (unseen) sampling frequencies. This is because the layers utilized in a
DNN are not designed for multiple sampling frequencies. In fact, the sampling
frequency has not been considered as a parameter of a DNN, but it is
implicitly given by the training dataset. In order to realize a DNN that
consistently works for any sampling frequency, a DNN must be designed as a
sampling-frequency-independent (SFI) network.

In this paper, we propose an SFI convolution layer for the handling of
arbitrary sampling frequencies by a single DNN. The key idea behind the
proposed layer is to consider the connection between a digital filter and a
convolution layer. From a signal processing viewpoint, we can interpret a
convolution layer as a collection of time-reversed digital finite impulse
response (FIR) filters. Therefore, a filter design technique can be utilized
to design a convolution layer. In this paper, we consider the impulse
invariant method (see Chap. 7 in [16]), in which a digital filter is designed
by sampling an analog filter. On the basis of this analog-to-digital filter
conversion, we introduce latent analog filters into a convolution layer. Since
an analog filter is independent of sampling frequency, we can construct an SFI
convolution layer via the analog representation of a filter, where the impulse
invariant method determines its sampling frequency afterward. The proposed SFI
layer can be trained by parametrizing the analog filter as a differentiable
function. By incorporating the proposed layer into one of the state-of-the-art
source separation models, we also propose an SFI audio source separation
model.

##  II Conventional Models

###  II-A Conv-TasNet [3]

Conv-TasNet is a recent time-domain DNN for audio source separation that works
well for speech [3] and music source separation [17, 5]. Since the
architectures of Conv-TasNet are slightly different in those three papers, we
adopt the Conv-TasNet architecture for music source separation defined in
[17], as illustrated in Fig. 1(a). Conv-TasNet consists of a pair of an
encoder and a decoder and $C$ source-specific masking modules, where $C$
denotes the number of sources. The encoder and decoder imitate a traditional
time-frequency transform (e.g., the short-time Fourier transform) and its
inverse transform. The encoder transforms a monaural time-domain signal into
an $N$-channel latent representation by a one-dimensional (1D) convolution
layer (with kernel size $L$ and stride $W$) and the rectified linear unit
(ReLU). Each masking module estimates a mask for the target source from the
latent representation. It comprises $R$ convolution blocks, each of which
consists of $X$ 1D dilated convolution layers with an exponentially increasing
dilation factor. The details of the convolution block are shown in [3]. The
decoder converts the masked latent representations into the separated time-
domain signals by a 1D transposed convolution layer with kernel size $L$ and
stride $W$.

###  II-B Multi-phase Gammatone Filter [7]

In [7], the multi-phase gammatone filter (MP-GTF) was introduced to design the
weights of the convolution layer of the encoder of Conv-TasNet, which improved
the speech separation performance. The impulse response of the MP-GTF is given
by

$g^{\mathrm{(MP\text{-}GTF)}}(t)=a\,t^{p-1}e^{-2\pi bt}\cos{(2\pi ft+\phi)},$
(1)

where $a$ denotes the amplitude, $p$ the filter order, $b$ the bandwidth, $f$
the center frequency, and $\phi$ the phase shift. The parameter $b$ is given
by $b=\mathrm{ERB}(f)/1.57$, where $\mathrm{ERB}(f)=24.7+f/9.265$. By sampling
$L$ points from $g^{\mathrm{(MP\text{-}GTF)}}(t)$ for various $f$ and $\phi$,
we obtain $N$ discrete-time impulse responses of length $L$ and concatenate
them along the channel axis to form the weights as a $1\times N\times L$
tensor. The convolution layer of the encoder is followed by the ReLU
nonlinearity, which blocks the negative values of the MP-GTF output and
discards the information of the input signal. To avoid a lack of information,
$g^{\mathrm{(MP\text{-}GTF)}}(t)$ is used together with its phase-reversed
version, i.e., with the phase shift $\phi+\pi$ [7].

###  II-C Multiple-sampling-frequency Training [13, 14, 15, 17]

There exist a few methods of training DNNs using audio signals sampled at
multiple sampling frequencies [13, 14, 15, 17]. In [13, 14], an automatic
speech recognition (ASR) model was trained using audio signals sampled at $8$
and $16$ kHz, where the part of input features corresponding to the missing
frequency band was padded by zeros. In [15], to compensate the missing
frequency band, an ASR model was jointly trained with a bandwidth expansion
model. The music source separation model presented in [17] was constructed by
stacking three Conv-TasNets that account for sampling frequencies of $8$,
$16$, and $32$ kHz. The Conv-TasNets for $16$ and $32$ kHz estimate the source
signals of the target sampling frequencies, referring to the masked latent
representations obtained with the lower sampling frequencies.

While these training methods are valid for the trained sampling frequencies,
they are not guaranteed to work with unseen sampling frequencies. In contrast,
we explicitly define an SFI structure to handle any sampling frequency without
retraining as shown later in Section III.

##  III Proposed Model

###  III-A Sampling-frequency-independent (SFI) Convolution Layer

To realize an SFI network, we introduce latent analog filters and analog-to-
digital filter conversion into a convolution layer. By interpreting the
weights of a convolution layer as a collection of time-reversed digital FIR
filters, we consider them from a signal processing viewpoint. The digital
filters are inherently sampling frequency dependent, whereas analog filters
are SFI owing to their definition in the continuous time domain. Focusing on
this fact, we introduce latent analog filters behind a convolution layer so
that its weights can be adjusted using the sampling frequency of an input
signal.

As shown in Fig. 1(c), the proposed layer consists of the usual 1D convolution
layer and impulse responses of $M^{\mathrm{(in)}}M^{\mathrm{(out)}}$ analog
filters defined in the continuous time domain, where $M^{\mathrm{(in)}}$ and
$M^{\mathrm{(out)}}$ are the input and output channel sizes, respectively. The
generating process of the proposed layer consists of three steps. Given the
sampling frequency of an input signal, the proposed layer (i) generates a
discrete-time impulse response of length $L$ from each analog filter, (ii)
stacks the time-reversed versions of these discrete-time impulse responses to
form the weights as an $M^{\mathrm{(in)}}\times M^{\mathrm{(out)}}\times L$
tensor, and (iii) works as the usual convolution layer using them. Since steps
(i) and (ii) depend only on the sampling frequency and the continuous-time
impulse responses, we only need to perform them once (before the features are
input) whenever the sampling frequency changes.

For step (i), we employ the impulse invariant method to generate digital FIR
filters from their analog counterparts. Note that while this method is
originally for the design of infinite impulse response filters, we can use it
for the digital FIR filter design. Let us denote the sampling period as $T$, a
discrete time index as $l=1,\cdots,L$, and the continuous time as
$t\in\mathbb{R}$. The impulse invariant method generates a discrete-time
impulse response $h[l]$ from an analog filter $g(t)$ so that the sampled
instants coincide:

$h[l]=Tg(lT).$ (2)

Changing $T$ yields an impulse response for different sampling frequencies
$1/T$. By stacking the generated impulse responses, the weights for the
convolution layer is obtained in step (ii). Similarly, an SFI version of a
transposed convolution layer (SFI transposed convolution layer) is given by
changing the convolution layer in the SFI convolution layer to the transposed
convolution layer.

For the analog filter $g(t)$, we can use the MP-GTF in Eq. (1). The
continuous-time impulse responses can be different for each channel, and
hence, hereafter, a channel subscript $m$ is added to
$g^{\mathrm{(MP\text{-}GTF)}}(t)$ and its parameters:
$g^{\mathrm{(MP\text{-}GTF)}}_{m}(t),a_{m},p_{m},b_{m},f_{m}$, and $\phi_{m}$.
Whereas all parameters of $g^{\mathrm{(MP\text{-}GTF)}}_{m}(t)$ were fixed in
[7], we propose to train $f_{m}$ and $\phi_{m}$ jointly with the other DNN
components by the commonly-used backpropagation algorithm.

The gradient of $h[l]$ can be computed in the same manner as the usual
convolution layer. Since the gradient of $h[l]$ equals that of $g(lT)$ owing
to Eq. (2) and $g^{\mathrm{(MP\text{-}GTF)}}_{m}(t)$ is differentiable with
$f_{m}$ and $\phi_{m}$, the gradients of the trainable parameters of
$g^{\mathrm{(MP\text{-}GTF)}}_{m}(t)$ can be computed by the chain rule. These
computations can be easily implemented by defining the forward computation
process of the proposed layer owing to the automatic differentiation mechanism
installed in modern deep learning frameworks (e.g., PyTorch and TensorFlow).

###  III-B Aliasing Reduction Technique

Since the impulse invariant method simply performs sampling to an analog
filter, aliasing occurs in the obtained digital filters. As reported in [18,
19, 6], aliasing causes degradation of the DNN performance, and thus we
introduce an aliasing reduction technique. Since the energy of aliased
components is distributed above the Nyquist frequency, we propose to set the
weights of the $m$th channel to zero whenever the center frequency $f_{m}$ of
the corresponding analog filter is above the Nyquist frequency. This aliasing
reduction technique is important for the proposed layer when it is utilized
with low sampling frequencies, as shown later in Section IV.

###  III-C Application of Proposed Layers to Conv-TasNet

As shown in Fig. 1(b), we build an SFI audio source separation model by
incorporating the proposed layers into Conv-TasNet [3]. The convolution layer
of the encoder and the transposed convolution layer of the decoder are
respectively replaced with the SFI convolution and transposed convolution
layers. The masking modules are the same as in [3].

For our model, we should modify the kernel size $L$ and stride $W$ in
accordance with the sampling frequency during inference. As described in
Section II-A, the encoder and decoder can be interpreted as the time-frequency
transform and its inverse transform. With this interpretation, $L$ and $W$
correspond to the frame length and the frame shift, respectively. Hence, as
the sampling frequency doubles, $L$ and $W$ should double to make the
representation consistent for the masking module. For this reason, we
determine $L$ and $W$ for each target sampling frequency to keep the frame
length and shift unchanged in the continuous time domain. By replacing all
convolution layers in the masking modules with the SFI convolution layers,
this issue might be resolved. However, we left it as a future work because
additional care regarding the combination with other layers (e.g., group
normalization [20]) must be considered.

##  IV Experimental Evaluation

TABLE I: Features of proposed methods and Conv-TasNet
---
Method | $g_{m}(t)$ | Samp. freq. | Aliasing
adapt. | reduction
Conv-Tasnet [3] | - | No | -
T-MP-GTF | $g_{m}^{\mathrm{(MP\text{-}GTF)}}(t)$ | No | No
Proposed | $g_{m}^{\mathrm{(MP\text{-}GTF)}}(t)$ | Yes | No
Proposed+ | $g_{m}^{\mathrm{(MP\text{-}GTF)}}(t)$ | Yes | Yes
TABLE II: Hyperparameters of masking modules used in experiments
---
Symbol | Description | Value
$N$ | # of channels of latent representation | $440$
$B$ | # of channels in bottleneck and | $160$
residual paths’ $1\times 1$ convolution blocks
$Sc$ | # of channels in skip-connection paths’ | $160$
$1\times 1$ convolution blocks
$H$ | # of channels in convolution blocks | $160$
$P$ | Kernel size in convolution blocks | $3$

Figure 2: SDRs of Conv-TasNet and proposed models for test data at various
sampling frequencies. These SDRs and error bars respectively denote averages
and standard errors over results obtained with four random seeds. Red line
shows trained sampling frequency.

###  IV-A Experimental Settings

To evaluate the efficacy of the proposed method, we conducted music source
separation on the MUSDB18-HQ dataset [21], which consists of $86$ training,
$14$ validation, and $50$ test tracks. Each track contains separate recordings
of four musical instruments (vocals, bass, drums, and other), i.e., $C=4$. The
training and validation tracks were down-sampled to $16$ kHz, and we created
the test data by down- and up-sampling the test tracks to several target
sampling frequencies, $F_{s}^{\mathrm{(target)}}=8,12,\ldots,48$ kHz. As an
evaluation metric, we used the median signal-to-distortion ratios (SDRs)
computed with the BSSEval v4 toolkit [22].

We used the same data augmentation techniques as those of [17]: random
cropping of the $8$ s training audio segments, random amplification within
$[0.75,1.25]$, random selection of the left or the right channel, and random
intertrack shuffling of the instruments in half of the minibatch. We also
applied standardization (zero mean and unit variance) to the tracks.

We compared the proposed model (Proposed) and that using the aliasing
reduction technique (Proposed+) with Conv-TasNet and its variant (T-MP-GTF)
whose encoder and decoder instead use the trainable extension of MP-GTF as
their weights of the convolution and transposed convolution layers,
respectively. T-MP-GTF was included to separately evaluate the trainable
extension of the MP-GTF and the sampling frequency adaptation. We applied all
models to the audio signals of the unseen sampling frequencies without
resampling them at the trained sampling frequency in order to examine the
effects of the sampling frequency mismatch and the proposed sampling frequency
adaptation. Table I summarizes the features of these models. For Proposed and
Proposed+, we determined $L$ and $W$ to be $5.0$ and $2.5$ ms at the sampling
frequency of $16$ kHz, respectively, as described in Section III-C, whereas we
set $L=80$ and $W=40$ for the other models. For all models, we set $X=6$ and
$R=2$. The hyperparameters of the masking modules are shown in Table II, where
the symbols correspond to those used in the literature of Conv-TasNet (see
Table 1 in [3]).

For $g_{m}^{\mathrm{(MP\text{-}GTF)}}(t)$, we trained $f_{m}$ and $\phi_{m}$
for $m=1,\ldots,220$ jointly with the entire network, and constrained these
parameters for the other $m$s so that $f_{m+220}=f_{m}$ and
$\phi_{m+220}=\phi_{m}+\pi$, as described in Section III-A. We initialized
$f_{m}$ and $\phi_{m}$ as in [7]; let us denote $48$ frequencies distributed
uniformly in the equivalent rectangular bandwidth (ERB) scale [23] from $50$
to $8000$ Hz by $f^{\mathrm{(center)}}_{i}$, where $i=1,\cdots,48$ is the
center frequency index and $f_{i}<f_{i+1}$ for all $i$. We initialized $f_{m}$
as $f_{m}=f^{\mathrm{(center)}}_{\lfloor m/K\rfloor+1}$ for $m=1,\ldots,140$
($m=141,\ldots,220$), where $K=5$ ($K=4$, respectively). The phase shifts
$\phi_{m}$ of the filters with the same $f^{\mathrm{(center)}}_{i}$ were
initialized to be uniformly distributed in $[0,\pi)$. The other parameters
were set as $a_{m}=1$, $p_{m}=2$. As in [7], these filters were normalized so
that they have the same $l^{2}$ norm.

For training, we used the RAdam optimizer [24] with a weight decay rate of
$5.0\times 10^{-4}$ and the Lookahead mechanism [25] with a synchronization
period of $6$ and a slow weights step of $0.5$. Gradient clipping with the
maximum $L_{2}$-norm of $5$ was applied. The learning rate scheduler presented
in [26] was employed with an initial learning rate of $1.0\times 10^{-3}$ and
a restart period of $200\,000$ iterations. We trained each model with a batch
size of $12$ for $250$ epochs, using the negative scale-invariant source-to-
noise ratio as the loss function, and selected the model with the lowest
validation loss. We applied the trained models to the left and right channels
of the test tracks separately, and scaled the source estimates using
instrument-wise factors to minimize the mean squared error between the input
mixture and the sum of all instrument estimates, resulting from the scale
invariance of the loss function [17].

###  IV-B Results

Fig. 2 shows the separation performances for the test data with sampling
frequencies from $8$ to $48$ kHz. At the trained sampling frequency, $16$ kHz,
the proposed models including T-MP-GTF achieved higher SDRs than Conv-TasNet
for all instruments, showing the effectiveness of the proposed trainable
extension of the MP-GTF. Interestingly, the proposed models gave much lower
standard errors than Conv-TasNet, which was not reported in [7]. This
observation reveals that the use of the trainable MP-GTF makes Conv-TasNet
robust to the initialization of the DNN parameters.

As the sampling frequency moved away from $16$ kHz, the SDRs of Conv-TasNet
and T-MP-GTF greatly decreased and they failed to separate the sources. By
contrast, the proposed models with the sampling frequency adaptation, Proposed
and Proposed+, provided similar SDRs for the $12$\- to $48$-kHz-sampled data
and outperformed the other models by a large margin, particularly at $20$ kHz
and higher, even though they were trained only with the $16$-kHz-sampled data.
This result clearly shows that the proposed sampling frequency adaptation
plays a crucial role in achieving the consistent performance.

Fig. 3 shows the magnitudes of the frequency responses of the filters of the
trained SFI convolution layer at sampling frequencies of $8,16,$ and $32$ kHz.
We can confirm that the filters of the $16$ and $32$ kHz sampling frequencies
exhibited consistent frequency responses. Importantly, the filters of the $32$
kHz sampling frequency blocked the frequency components higher than around $8$
kHz. Nevertheless, the proposed model achieved a consistent performance at the
sampling frequencies higher than $8$ kHz, which should be because the dominant
frequency components of the music signals were distributed below $8$ kHz.

For the sampling frequency of $8$ kHz, aliasing occurred from $4$ kHz (see
Fig. 3(a)). This resulted in the performance degradation of Proposed (the
proposed method without the aliasing reduction technique). By contrast,
Proposed+ showed a consistent performance at all sampling frequencies,
demonstrating the effectiveness of the proposed aliasing reduction technique
when the sampling frequency is reduced. For drums and other, Proposed+ gave
slightly lower SDRs than Proposed at the sampling frequency of $12$ kHz, which
might be because the filters with center frequencies near the Nyquist
frequency are helpful for the separation. A further investigation of this
observation remains as a future work.

Figure 3: Magnitudes of frequency responses of first $220$ filters of trained
SFI convolution layer at sampling frequencies of $8,16,$ and $32$ kHz.

##  V Conclusion

We proposed an SFI convolution layer that can be adjusted to an arbitrary
sampling frequency. We focused on the fact that the weights of a convolution
layer can be seen as a collection of digital FIR filters. We explicitly
defined the weights generation process of the convolution layer from the
latent analog filters based on the impulse invariant method. Since the analog
filters do not depend on sampling frequency, the proposed layers can generate
consistent weights for arbitrary sampling frequencies. Furthermore, we built
an SFI audio source separation model by incorporating the proposed layers into
the encoder and decoder of Conv-TasNet. We showed, through music source
separation experiments, that even when trained with only the audio signals
sampled at a specific sampling frequency, the proposed model consistently
worked well with not only the trained sampling frequency but also unseen ones.
Since the proposed layer is a general component for audio processing, it would
also be useful for various audio applications such as speech separation [27,
1, 2, 3, 7].

## References

  * [1] Z. Wang, J. Le Roux, and J. R. Hershey,  “Alternative objective functions for deep clustering,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, 2018, pp. 686–690.
  * [2] D. Yu, M. Kolbæk, Z. Tan, and J. Jensen,  “Permutation invariant training of deep models for speaker-independent multi-talker speech separation,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, 2017, pp. 241–245.
  * [3] Y. Luo and N. Mesgarani,  “Conv-TasNet: Surpassing ideal time–frequency ,magnitude masking for speech separation,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 27, no. 8, pp. 1256–1266, 2019.
  * [4] D. Stoller, S. Ewert, and S. Dixon,  “Wave-U-Net: A Multi-Scale neural network for end-to-end audio source separation,”  in Proceedings of International Society for Music Information Retrieval Conference, 2018, pp. 334–340.
  * [5] A. Défossez, N. Usunier, L. Bottou, and F. Bach,  “Music source separation in the waveform domain,”  arXiv preprint arXiv:1911.13254, 2019.
  * [6] T. Nakamura and H. Saruwatari,  “Time-domain audio source separation based on Wave-U-Net combined with discrete wavelet transform,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, 2020, pp. 386–390.
  * [7] D. Ditter and T. Gerkmann,  “A Multi-Phase Gammatone Filterbank for speech separation via TasNet,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, 2020, pp. 36–40.
  * [8] H. Liu, L. Xie, J. Wu, and G. Yang,  “Channel-wise subband input for better voice and accompaniment separation on high resolution music,”  in Proceedings of INTERSPEECH, 2020.
  * [9] D. Takeuchi, K. Yatabe, Y. Koizumi, Y. Oikawa, and N. Harada,  “Data-driven design of perfect reconstruction filterbank for DNN-based sound source enhancement,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, 2019, pp. 596–600.
  * [10] F. Krebs, S. Böck, M. Dorfer, and G. Widmer,  “Downbeat tracking using beat synchronous features with recurrent neural networks,”  in Proceedings of International Society for Music Information Retrieval Conference, 2016, pp. 129–135.
  * [11] S. Sigtia, E. Benetos, and S. Dixon,  “An end-to-end neural network for polyphonic piano music transcription,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 24, no. 5, pp. 927–939, 2016.
  * [12] F. Pedersoli, G. Tzanetakis, and K. M. Yi,  “Improving music transcription by pre-stacking a U-Net,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, 2020, pp. 506–510.
  * [13] D. Yu, M. Seltzer, J. Li, J. Huang, and F. Seide,  “Feature learning in deep neural networks - studies on speech recognition,”  in Proceedings of International Conference on Learning Representations, 2013.
  * [14] A. Narayanan, A. Misra, K. C. Sim, G. Pundak, A. Tripathi, M. Elfeky, P. Haghani, T. Strohman, and M. Bacchiani,  “Toward domain-invariant speech recognition via large scale training,”  in IEEE Spoken Language Technology Workshop, 2018, pp. 441–447.
  * [15] J. Gao, J. Du, and E. Chen,  “Mixed-bandwidth cross-channel speech recognition via joint optimization of DNN-based bandwidth expansion and acoustic modeling,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 27, no. 3, pp. 559–571, 2019.
  * [16] A. V. Oppenheim, J. R. Buck, and R. W. Schafer,  Discrete-time signal processing,  Prentice Hall, Upper Saddle River, NJ, 2001.
  * [17] D. Samuel, A. Ganeshan, and J. Naradowsky,  “Meta-learning extractors for music source separation,”  in Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, 2020, pp. 816–820.
  * [18] M. D. Zeiler and R. Fergus,  “Visualizing and understanding convolutional networks,”  in Proceedings of European Conference on Computer Vision, 2014, pp. 818–833.
  * [19] Y. Gong and C. Poellabauer,  “Impact of aliasing on deep CNN-based end-to-end acoustic models,”  in Proceedings of INTERSPEECH, 2018, pp. 2698–2702.
  * [20] Y. Wu and K. He,  “Group normalization,”  in Proceedings of European Conference on Computer Vision, 2018.
  * [21] Z. Rafii, A. Liutkus, F.-R. Stöter, S. I. Mimilakis, and R. Bittner,  “MUSDB18-HQ - an uncompressed version of MUSDB18,” 2019.
  * [22] F.-R. Stöter, A. Liutkus, and N. Ito,  “The 2018 signal separation evaluation campaign,”  in Proceedings of International Conference on Latent Variable Analysis and Signal Separation, 2018, pp. 293–305.
  * [23] V. Hohmann,  “Frequency analysis and synthesis using a gammatone filterbank,”  Acta Acustica united with Acustica, vol. 88, no. 03, pp. 433–442, 2002.
  * [24] L. Liu, H. Jiang, P. He, W. Chen, X. Liu, J. Gao, and J. Han,  “On the variance of the adaptive learning rate and beyond,”  in Proceedings of International Conference on Learning Representations, 2020.
  * [25] M. Zhang, J. Lucas, J. Ba, and G. Hinton,  “Lookahead Optimizer: k steps forward, 1 step back,”  in Proceedings of Advances in Neural Information Processing Systems, 2019, pp. 9597–9608.
  * [26] I. Loshchilov and F. Hutter,  “SGDR: Stochastic gradient descent with warm restarts,”  in Proceedings of International Conference on Learning Representations, 2017.
  * [27] H. Sawada, N. Ono, H. Kameoka, D. Kitamura, H. Saruwatari,,  “A review of blind source separation methods: two converging routes to ILRMA originating from ICA and NMF,”  in APSIPA Transactions on Signal and Information Processing, vol.8, no. e12, 14 pages, 2019.


# Over-the-Air Computation with Spatial-and-Temporal Correlated Signals

Wanchun Liu, Xin Zang, Branka Vucetic and Yonghui Li

###### Abstract

††The authors are with School of Electrical and Information Engineering, the
University of Sydney, Australia. Emails: {wanchun.liu, xin.zang,
branka.vucetic<EMAIL_ADDRESS>(_Yonghui Li is the corresponding
author._)

Over-the-air computation (AirComp) leveraging the superposition property of
wireless multiple-access channel (MAC), is a promising technique for effective
data collection and computation of large-scale wireless sensor measurements in
Internet of Things applications. Most existing work on AirComp only considered
computation of spatial-and-temporal independent sensor signals, though in
practice different sensor measurement signals are usually correlated. In this
letter, we propose an AirComp system with spatial-and-temporal correlated
sensor signals, and formulate the optimal AirComp policy design problem for
achieving the minimum computation mean-squared error (MSE). We develop the
optimal AirComp policy with the minimum computation MSE in each time step by
utilizing the current and the previously received signals. We also propose and
optimize a low-complexity AirComp policy in closed form with the performance
approaching to the optimal policy.

###### Index Terms:

Over-the-air computing, wireless sensor networks, multiple-access channel,
spacial-and-temporal correlation.

##  I Introduction

In the past few years, over-the-air computation (AirComp) has emerged as a
promising solution for large-scale wireless sensor data fusion in Internet of
Things applications [1, 2, 3, 4, 5, 6]. In general, an AirComp system consists
of multiple sensors and a receiver to perform data fusion by computing a pre-
determined function of the sensors’ measurement signals. Each sensor of the
AirComp system sends its _pre-processed_ measurement signal simultaneously to
the receiver over a multiple-access channel (MAC). By leveraging the
_superposition_ property of MAC and applying a _post-processing_ function to
the received sum (superimposed) signal, the receiver obtains the desired
function output of the sensors’ signals (see e.g. [2] for details).

From a computation theoretic perspective, pre-processing and post-processing
function design problems of AirComp systems for different fusion functions
have been comprehensively investigated in [7, 8, 9]. From a wireless
communications perspective, due to the imperfection of practical MAC with non-
zero receiver noise, unequal channel coefficients and limited transmit power
budgets, the output signal of a MAC is not a perfect sum of the pre-processed
sensor signals as expected. Thus, most recent work focuses on optimal AirComp
transmitter and receiver policy design for achieving the minimum computation
mean-squared error (MSE) of the sum of pre-processed signals [1, 2, 3, 4]. In
[1], the computation MSE of the sum signal under imperfect channel state
information was investigated. In [2], the optimal single-antenna AirComp
design and the scaling law analysis of computation MSE and transmit power
consumption were derived in terms of the number of sensors. In [3] and [4],
the transmitting and receiving beamforming of multi-antenna AirComp systems
were considered to minimize the computation MSE. In particular, AirComp has
been applied to emerging wireless applications such as wireless distributed
machine learning [5], and over-the-air consensus [6].

Most existing work on AirComp system assumed spatial-and-temporal independent
sensor signals for ease of analysis and optimization [1, 2, 3, 4]. However, in
practice, different sensor measurements are correlated and the current
measurements are normally related to the previous ones. Spatial-and-temporal
sensor signals are commonly considered in the literature of remote state
estimation (see e.g. [10, 11]), but have not been considered in AirComp
systems to the best of our knowledge. In particular, unlike the conventional
time-independent AirComp system, where only the current received signal is
used to estimate the sum signal of the sensors, the time correlation
introduces a new challenge in AirComp system design: how to effectively
utilize both the current and the sequence of previously received signals for
optimal computation of the sum of current sensor signals?

In this paper, we focus on AirComp system design with spatial-and-temporal
correlated sensor signals. The main contributions are summarized as follows.

  * •

We propose an AirComp system with spatial-and-temporal correlated sensor
signals for the first time in the literature, and formulate the optimal
AirComp policy design problem for achieving the minimum computation MSE.

  * •

We derive the optimal AirComp policy to achieve the minimum computation MSE in
each time step by utilizing the current and the previously received signals.

  * •

We propose a low-complexity AirComp policy, where the receiver stores the most
recent $(l+1)$ received signals and applies a linear filter to it as the
computation output. The optimal AirComp filter is derived in closed form. Once
an AirComp filter is given, we show that the optimal transmit power of each
sensor can be obtained by solving a classic quadratically constrained
quadratic program (QCQP) problem. An alternating minimization approach is
adopted to jointly design the AirComp filter and the transmit power.

_Notations_ : $\\{x\\}^{j}_{i}$ denotes the sequence
$\\{x_{i},x_{i+1},\dots,x_{j}\\}$. $\rho(\mathbf{A})$ is the spectral radius
of a square matrix $\mathbf{A}$, i.e., the largest absolute value of its
eigenvalues. $\mathsf{diag}\\{\mathbf{v}\\}$ denotes a diagonal matrix
generated by the vector $\mathbf{v}$. $\mathbf{1}$ is a all-one column vector.
$\mathbf{0}_{i,j}$ is a $i\times j$ all-zero matrix. $[\mathbf{A}]_{i,j}$ is
the element in the $i$th row and $i$th column of $\mathbf{A}$.

##  II AirComp System with Spatial-and-Temporal Correlated Signals

We consider a single-antenna AirComp system with $K$ sensors and a receiver.

###  II-A Spatial-and-Temporal Correlated Signals

Given a sampling period $T$, sensor $k$’s signal at time slot $t$ is
$x_{k,t}\in\mathbb{R},\forall k\in\mathcal{K},t\in\mathbb{N}_{0}$, where
$\mathcal{K}\triangleq\\{1,\cdots,K\\}$. We assume that the sensors’ signals
are spatial-and-temporal correlated, and are modeled as a Gauss-Markov process
[12]:

$\mathbf{x}_{t+1}=\mathbf{A}\mathbf{x}_{t}+\mathbf{w}_{t},t\in\mathbb{N}_{0}$
(1)

where
$\mathbf{x}_{t}\triangleq[x_{1,t},x_{2,t},\dots,x_{K,t}]^{\top}\in\mathbb{R}^{K}$,
$\mathbf{A}\in\mathbb{R}^{K\times K}$ is the process parameter, and
$\mathbf{w}_{t}\triangleq[w_{1,t},w_{2,t},\dots,w_{K,t}]^{\top}\in\mathbb{R}^{K}$
is a white Gaussian noise with zero mean and a covariance matrix
$\mathbf{V}_{w}$. Thus, each sensor’s current signal depends on its and other
sensors’ previous signals. Further, it is assumed that $\rho(\mathbf{A})<1$,
i.e., $\mathbf{x}_{t}$ has a steady state distribution with zero mean and a
(bounded) converged steady-state covariance matrix
$\mathbf{V}_{x}\triangleq\mathsf{E}[\mathbf{x}_{t}\mathbf{x}^{\top}_{t}]$.
From [12], $\mathbf{V}_{x}$ of the Gauss-Markov process (1) is the unique
positive semi-definite solution of the discrete Lyapunov matrix equation below

$\mathbf{A}\mathbf{V}_{x}\mathbf{A}^{\top}-\mathbf{V}_{x}+\mathbf{V}_{w}=\mathbf{0}.$
(2)

In the scenario that $\rho(\mathbf{A})\geq 1$, $\mathbf{x}_{t}$ has an
unbounded covariance matrix when $t\rightarrow\infty$, which is beyond the
scope of the paper.

###  II-B AirComp System and Computation MSE

Each sensor linearly scales its signal by its Tx-scaling factor $b_{k}$, and
then sends $b_{k}x_{k,t}$ to the receiver simultaneously over a MAC at time
slot $t$. Let $h_{k}$ denote the channel coefficient between sensor $k$ and
the receiver. We focus on the scenario that the sensor signals change much
faster than the channel status. Thus, $h_{k}$ is assumed to be time-invariant
[13]. The received signal at time slot $t$ is [2]

$y_{t}=\sum_{k=1}^{K}h_{k}b_{k}x_{k,t}+z_{t},t\in\mathbb{N}_{0},$ (3)

where $z_{t}$ is the receiver-side additive white Gaussian noise (AWGN) with
zero mean and variance $\sigma^{2}_{z}$. For ease of notation but without loss
of generality, we use $b_{k}$ to denote $h_{k}b_{k}$ in the rest of the
analysis. Let $\mathbf{b}\triangleq[b_{1},\dots,b_{K}]^{\top}$ denote the Tx-
scaling vector. Then, (3) can be rewritten as

$y_{t}=\mathbf{b}^{\top}\mathbf{x}+z_{t}.$ (4)

Note that similar to the previous work [2, 14], we focus on a fundamental real
channel scenario, i.e., $h_{k},b_{k},y_{t},z_{t}\in\mathbb{R},\forall
k\in\mathcal{K}$, though the analytical framework and results in the paper can
be readily extended to the complex channel scenario by representing $y_{t}$
and $z_{t}$ as two-dimensional vectors, and $\mathbf{b}$ as a $K\times 2$
matrix with real and imaginary parts.

Different from the scenario with time-independent signals (see e.g., [3]), the
previously received signals $\\{y\\}^{t-1}_{0}$ are correlated with and hence
contain information of the current signal $\mathbf{x}_{t}$. Thus, it is
assumed that the previous received signals are stored at the receiver.

The AirComp system with time-correlated signals in (1) is designed to
calculate the sum of the sensors’ signals $\sum_{k=1}^{K}x_{k,t}$ in each time
slot based on the current and previously received signals $\\{y\\}^{t}_{0}$.
Let $\chi_{t}$ denote the calculated sum signal. Then, the computation
distortion of the AirComp system at time slot $t$ is measured by the
computation MSE defined as

$\mathsf{CMSE}_{t}\triangleq\mathsf{E}\left[\lvert\chi_{t}-\sum_{k=1}^{K}x_{k,t}\rvert^{2}\Big{|}\\{y\\}^{t}_{0}\right],\forall
t\in\mathbb{N}.$ (5)

In the following, we will investigate the optimal AirComp policy $\pi(\cdot)$
to compute $\chi_{t}=\pi(\\{y\\}^{t}_{0})$ with the smallest computation MSE.

##  III Optimal AirComp Policy Design

We introduce a virtual vector
$\mathbf{\tilde{x}}_{t}\triangleq[\tilde{x}_{1,t},\tilde{x}_{1,t},\dots,\tilde{x}_{K,t}]^{\top}\in\mathbb{R}^{K}$.
Although an AirComp policy is to find a scalar as the computation output,
without loss of generality, the policy can be treated as to find a vector
$\mathbf{\tilde{x}}_{t}$ in each time slot and calculate the sum of it, i.e.,
$\chi_{t}\triangleq\sum_{k=1}^{K}\tilde{x}_{k,t}$. To find the optimal
$\chi_{t}$ achieving the minimum computation MSE, one needs to find the
optimal virtual vector $\mathbf{\tilde{x}}_{t}$. From (5), the optimal AirComp
policy design problem is formulated as

$\min_{\mathbf{\tilde{x}}_{t}}\quad\mathsf{CMSE}_{t}\triangleq\mathsf{E}\left[\lvert\sum_{k=1}^{K}\tilde{x}_{k,t}-\sum_{k=1}^{K}x_{k,t}\rvert^{2}\Big{|}\\{y\\}^{t}_{0}\right],\forall
t\in\mathbb{N}.$ (6)

The target function in problem (6) can be rewritten as

$\displaystyle\mathsf{CMSE}_{t}=\mathsf{E}\left[\mathbf{1}^{\top}(\mathbf{\tilde{x}}_{t}-\mathbf{x}_{t})(\mathbf{\tilde{x}}_{t}-\mathbf{x}_{t})^{\top}\mathbf{1}\Big{|}\\{y\\}^{t}_{0}\right]$
(7)
$\displaystyle=\mathsf{E}\left[\mathbf{1}^{\top}\mathbf{\tilde{x}}_{t}\mathbf{\tilde{x}}_{t}^{\top}\mathbf{1}-2\times\mathbf{1}^{\top}\mathbf{\tilde{x}}_{t}\mathbf{{x}}_{t}^{\top}\mathbf{1}+\mathbf{1}^{\top}\mathbf{{x}}_{t}\mathbf{{x}}_{t}^{\top}\mathbf{1}\Big{|}\\{y\\}^{t}_{0}\right]$
$\displaystyle=\\!\\!\int_{\mathbb{R}^{K}}\\!\\!\\!\\!\left(\mathbf{1}^{\top}\mathbf{\tilde{x}}_{t}\mathbf{\tilde{x}}_{t}^{\top}\mathbf{1}-2\times\\!\mathbf{1}^{\top}\mathbf{\tilde{x}}_{t}\mathbf{{x}}_{t}^{\top}\mathbf{1}+\mathbf{1}^{\top}\mathbf{{x}}_{t}\mathbf{{x}}_{t}^{\top}\mathbf{1}\right)\\!\\!f\\!\\!\left(\mathbf{x}_{t}|\\{y\\}^{t}_{0}\right)\\!\\!\mathrm{d}\mathbf{x}_{t}$
$\displaystyle=\mathbf{1}^{\top}\mathbf{\tilde{x}}_{t}\mathbf{\tilde{x}}_{t}^{\top}\mathbf{1}\\!+\\!\\!\int_{\mathbb{R}^{K}}\\!\\!\\!\\!\left(-2\times\\!\mathbf{1}^{\top}\mathbf{\tilde{x}}_{t}\mathbf{{x}}_{t}^{\top}\mathbf{1}\\!+\\!\mathbf{1}^{\top}\mathbf{{x}}_{t}\mathbf{{x}}_{t}^{\top}\mathbf{1}\right)\\!\\!f\\!\\!\left(\mathbf{x}_{t}|\\{y\\}^{t}_{0}\right)\\!\\!\mathrm{d}\mathbf{x}_{t},$

where $f\left(\mathbf{x}_{t}|\\{y\\}^{t}_{0}\right)$ is the conditional
probability density function of $\mathbf{x}_{t}$.

It can be shown that $\mathsf{CMSE}_{t}$ is a convex function of
$\mathbf{\tilde{x}}_{t}$ as the Hessian matrix of $\mathsf{CMSE}_{t}$ is
$2\times\mathbf{1}\mathbf{1}^{\top}$, which is positive semi-definite. Since
any local minimum of a convex function is also a global minimum [15], we only
need to find a local minimum of $\mathsf{CMSE}_{t}$. Then, by letting the
derivative equal to zero, we have

$\frac{\partial\mathsf{CMSE}_{t}}{\partial\mathbf{\tilde{x}}_{t}}=2\times\mathbf{1}\mathbf{1}^{\top}\mathbf{\tilde{x}}_{t}-2\times\mathbf{1}\mathbf{1}^{\top}\int_{\mathbb{R}^{K}}\mathbf{x}_{t}f\left(\mathbf{x}_{t}|\\{y\\}^{t}_{0}\right)\mathrm{d}\mathbf{x}_{t}=\mathbf{0}.$
(8)

Although (8) has many solutions, it is clear that

$\mathbf{\tilde{x}}_{t}=\mathbf{\hat{x}}_{t}\triangleq\int_{\mathbb{R}^{K}}\mathbf{x}_{t}f\left(\mathbf{x}_{t}|\\{y\\}^{t}_{0}\right)\mathrm{d}\mathbf{x}_{t}$
(9)

is a solution achieving the minimum $\mathsf{CMSE}_{t}$. Note that
$\mathbf{\hat{x}}_{t}$ in (9) is the conditional expectation of the current
signal $\mathbf{x}_{t}$. Based on the linear estimation theory [12],
$\mathbf{\hat{x}}_{t}$ is obtained by the classic Kalman filter (KF) in an
iterative prediction-correction way as below:

$\displaystyle\textbf{Prediction:
}\hat{\mathbf{x}}_{t|t-1}=\mathbf{A}\hat{\mathbf{x}}_{t-1},$ (10a)
$\displaystyle\textbf{Prediction MSE Matrix:
}\\!\mathbf{M}_{t|t-1}\\!\\!=\\!\mathbf{A}\mathbf{M}_{t-1}\mathbf{A}^{\top}\\!\\!\\!+\\!\mathbf{V}_{w},$
(10b) $\displaystyle\textbf{Kalman Gain Vector:
}\mathbf{K}_{t}\\!=\\!{\mathbf{M}_{t|t-1}\mathbf{b}}\left({\sigma_{z}^{2}+\\!\mathbf{b}^{\top}\mathbf{M}_{t|t-1}\mathbf{b}}\right)^{\\!-1}\\!,$
(10c) $\displaystyle\textbf{Correction:
}\hat{\mathbf{x}}_{t}=\hat{\mathbf{x}}_{t|t-1}+\mathbf{K}_{t}(y_{t}-\mathbf{b}^{\top}\hat{\mathbf{x}}_{t|t-1}),$
(10d) $\displaystyle\textbf{Estimation MSE matrix:
}\mathbf{M}_{t}\\!=\\!(\mathbf{I}\\!-\\!\mathbf{K}_{t}\mathbf{b}^{\top}\\!)\mathbf{M}_{t|t-1}.$
(10e)

In particular, $\mathbf{M}_{t}$ is the estimation error covariance matrix of
$\mathbf{\hat{x}}_{t}$, i.e.,

$\mathbf{M}_{t}\triangleq\mathsf{E}\left[(\mathbf{\hat{x}}_{t}-\mathbf{x}_{t})(\mathbf{\hat{x}}_{t}-\mathbf{x}_{t})^{\top}\Big{|}\\{y\\}^{t}_{0}\right].$
(11)

Taking (9) into (7), the minimum $\mathsf{CMSE}_{t}$ is equal to
$\mathbf{1}^{\top}\mathbf{M}_{t}\mathbf{1}$. Therefore, we have the following
result.

######  Theorem 1.

Given the Tx-scaling factor vector $\mathbf{b}$, the optimal AirComp policy
for achieving the minimum computation MSE is to first apply the KF defined in
(10) to estimate the sensor signal vector $\mathbf{x}_{t}$ and then calculate
the sum of the estimate $\mathbf{\hat{x}}_{t}$. The minimum computation MSE is
$\mathbf{1}^{\top}\mathbf{M}_{t}\mathbf{1}$, where $\mathbf{M}_{t}$ is give in
(10e).

######  Remark 1.

It is non-trivial that the KF-based state estimation method leads to the
optimal computation of the sum of the original signals. The KF is optimally
designed to achieve the minimum state estimation MSE, i.e.,
$\mathsf{E}\left[(\mathbf{\hat{x}}_{t}-\mathbf{x}_{t})^{\top}(\mathbf{\hat{x}}_{t}-\mathbf{x}_{t})\right]=\mathsf{Tr}\left(\mathbf{M}_{t}\right)$
[12], which is not equal to the design target of the AirComp problem (6),
i.e.,
$\mathsf{E}\left[\mathbf{1}^{\top}(\mathbf{\hat{x}}_{t}-\mathbf{x}_{t})(\mathbf{\hat{x}}_{t}-\mathbf{x}_{t})^{\top}\mathbf{1}\right]=\mathbf{1}^{\top}\mathbf{M}_{t}\mathbf{1}$.
Thus, without the rigorous proof of Theorem 1, one cannot tell whether the KF-
based estimation-then-sum method is the optimal or not.

Also, it is interesting to see that in the AirComp system with time correlated
signals, one needs to recover (estimate) each sensor’s signal first based on
the previous and current observations (received signals) and then calculate
the sum signal; while in the conventional time independent scenario [2, 14],
the optimal policy is to estimate the sum signal directly rather than the
original signals.

Taking (10b) into (10c) and then into (10d), it is clear that the computation
of $\mathbf{\hat{x}}_{t}$ requires $K\times K$ matrix multiplications. So the
computation complexity of the optimal AirComp policy in each time step is
$\mathcal{O}(K^{3})$ based on the big O notation, when the total number of
sensors $K$ is very large. Further, since $\mathbf{M}_{t}$ is a complex and
implicit function of the Tx-scaling factor $\mathbf{b}$, it is difficult to
design the optimal $\mathbf{b}$ for achieving the minimum computation MSE,
$\mathbf{1}^{\top}\mathbf{M}_{t}\mathbf{1}$.

In the following, we will propose and optimize a low complexity AirComp
policy, based on which an optimal Tx-scaling factor can be obtained.

##  IV Low-Complexity AirComp Policy and Transmission Power Control

###  IV-A Low-Complexity AirComp Policy

We propose a low-complexity AirComp policy: the receiver stores the most
recent $(l+1)$ received signals

$\mathbf{y}_{t}\triangleq[y_{t-l},y_{t-l+1},\dots,y_{t}]^{\top},\ t\geq
l,l\geq 1,$ (12)

and applies a linear filter
$\mathbf{g}\triangleq[g_{0},g_{1},\dots,g_{l}]^{\top}\in\mathbb{R}^{l+1}$ to
it as the computation output $\chi_{t}$, i.e.,

$\chi_{t}=\mathbf{g}^{\top}\mathbf{y}_{t}.$ (13)

Intuitively, a longer filter utilizing more previous received signals achieves
a better computation MSE. Compared with the optimal AirComp policy in Theorem
1 with third-degree polynomial complexity, the policy (13) has a linear
computation complexity of $\mathcal{O}(l)$ and is invariant with the total
number of sensors.

The computation MSE of the low-complexity policy is

$\widetilde{\mathsf{CMSE}}_{t}=\mathsf{E}\left[|\mathbf{g}^{\top}\mathbf{y}_{t}-\mathbf{1}^{\top}\mathbf{x}_{t}|^{2}\right].$
(14)

Using the Gauss-Markov property of $\\{\mathbf{x}_{t}\\}$ in (1), the $i$th
previously received signal can be written as

$y_{t-i}=\mathbf{b}^{\top}\mathbf{A}^{i}\mathbf{x}_{t-l}+c_{i,t},i\in\\{0,1,\dots,l\\}$
(15)

where

$c_{i,t}=z_{t-l+i}+\sum_{m=1}^{i}\mathbf{b}^{\top}\mathbf{A}^{m-1}\mathbf{w}_{t-l+m-1}.$
(16)

Taking (15) into (12), we have

$\mathbf{y}_{t}=\mathbf{Mx}_{t-l}+\mathbf{c}_{t},$ (17)

where $\mathbf{c}_{t}\triangleq[c_{0,t},c_{2,t},\dots,c_{l,t}]^{\top}$, and

$\mathbf{M}\triangleq\left[(\mathbf{A}^{0})^{\top}\mathbf{b},(\mathbf{A}^{1})^{\top}\mathbf{b},\dots,(\mathbf{A}^{l})^{\top}\mathbf{b}\right]^{\top}.$
(18)

Further, $\mathbf{x}_{t}$ can be represented as a function of
$\mathbf{x}_{t-l}$ and the noise terms as

$\mathbf{x}_{t}=\mathbf{A}^{l}\mathbf{x}_{t-l}+\sum_{i=1}^{l}\mathbf{A}^{l-i}\mathbf{w}_{t-l+i-1}.$
(19)

Taking (17) and (19) into (14), it can be obtained that

$\displaystyle\widetilde{\mathsf{CMSE}}_{t}$ (20)
$\displaystyle=\mathsf{E}\\!\\!\left[\\!\Big{|}\mathbf{g}^{\top}\\!(\mathbf{Mx}_{t-l}+\mathbf{c}_{t})\\!-\\!\mathbf{1}^{\top}\\!(\mathbf{A}^{l}\mathbf{x}_{t-l}\\!+\\!\sum_{i=1}^{l}\\!\\!\mathbf{A}^{l-i}\mathbf{w}_{t-l+i-1})\Big{|}^{2}\\!\right].$

To find the optimal low-complexity AirComp filter $\mathbf{g}$, we have the
following optimization problem:

$\min_{\mathbf{g}}\qquad\eqref{eq:mse_2}.$ (21)

By introducing and deriving the covariance matrices
$\mathbf{V}_{c}\triangleq\mathsf{E}\left[\mathbf{c}_{t}\mathbf{c}_{t}^{\top}\right]$
and
$\mathbf{C}_{i}\triangleq\mathsf{E}\left[\mathbf{c}_{t}\mathbf{w}_{t-i}^{\top}\right],i\in\\{1,2,\dots,l\\}$,
with the properties that $\\{z_{t}\\}$ and $\\{\mathbf{w}_{t}\\}$ are
independent white Gaussian processes, the target function can be simplified as

$\displaystyle\widetilde{\mathsf{CMSE}}_{t}=\mathbf{g}^{\top}\mathbf{M}\mathbf{V}_{x}\mathbf{M}^{\top}\mathbf{g}\\!-\\!\mathbf{g}^{\top}\mathbf{M}\mathbf{V}_{x}(\mathbf{A}^{l})^{\top}\mathbf{1}\\!-\\!\mathbf{1}^{\top}\mathbf{A}^{l}\mathbf{V}_{x}\mathbf{M}^{\top}\mathbf{g}$
(22)
$\displaystyle+\mathbf{1}^{\top}\mathbf{A}^{l}\mathbf{V}_{x}(\mathbf{A}^{l})^{\top}\mathbf{1}+\mathbf{g}^{\top}\mathbf{V}_{c}\mathbf{g}+\sum_{i=1}^{l}\mathbf{1}^{\top}\mathbf{A}^{l-i}\mathbf{V}_{w}(\mathbf{A}^{l-i})^{\top}\mathbf{1}$
$\displaystyle-\mathbf{g}^{\top}\sum_{i=1}^{l}\mathbf{C}_{l+1-i}(\mathbf{A}^{l-i})^{\top}\mathbf{1}-\mathbf{1}^{\top}\sum_{i=1}^{l}\mathbf{A}^{l-i}\mathbf{C}_{l+1-i}^{\top}\mathbf{g},$

which is time independent. Thus, we drop the time index in
$\widetilde{\mathsf{CMSE}}_{t}$. It is clear that $\widetilde{\mathsf{CMSE}}$
has a quadratic form and is strictly convex, and the optimal $\mathbf{g}$ can
be obtained directly by solving the linear equation:

$\begin{split}&\frac{\partial\widetilde{\mathsf{CMSE}}}{\partial\mathbf{g}}\\!=\\!\mathbf{0}\\!=\\!2\mathbf{g}^{\top}\mathbf{M}\mathbf{V}_{x}\mathbf{M}^{\top}\\!-\\!\mathbf{1}^{\top}\mathbf{A}^{l}\mathbf{V}_{x}^{\top}\mathbf{M}^{\top}\\!-\\!\mathbf{1}^{\top}\mathbf{A}^{l}\mathbf{V}_{x}^{\top}\mathbf{M}^{\top}\\\
&+2\mathbf{g}^{\top}\mathbf{V}_{c}-\mathbf{1}^{\top}\sum_{i=1}^{l}\mathbf{A}^{l-i}\mathbf{C}_{l+1-i}^{\top}-\mathbf{1}^{\top}\sum_{i=1}^{l}\mathbf{A}^{l-i}\mathbf{C}_{l+1-i}^{\top}.\end{split}$
(23)

Then, we have the following result.

######  Theorem 2.

Given the Tx-scaling factor $\mathbf{b}$, the optimal length-$(l+1)$ linear
filter $\mathbf{g}$ of the low-complexity AirComp policy is

$\mathbf{g}^{\top}\\!=\\!\\!\left(\\!\mathbf{1}^{\top}\mathbf{A}^{l}\mathbf{V}_{x}^{\top}\mathbf{M}^{\top}\\!\\!+\\!\mathbf{1}^{\top}\\!\sum_{i=1}^{l}\mathbf{A}^{l-i}\mathbf{C}_{l+1-i}^{\top}\\!\right)\\!\\!(\mathbf{M}\mathbf{V}_{x}\mathbf{M}^{\top}\\!+\\!\mathbf{V}_{c})^{-1},$
(24)

where

$\mathbf{V}_{c}=\mathsf{diag}\\{\underbrace{\sigma^{2}_{z},\dots,\sigma^{2}_{z}}_{l+1}\\}+\mathbf{\tilde{V}}_{c},$
(25) $[\mathbf{\tilde{V}}_{c}]_{i,j}=\left\\{\begin{aligned} &0,&&i=0\text{ or
}j=0,\\\ &\mathbf{b}^{\top}\mathbf{U}_{i,j}\mathbf{b},&&1\leq i\leq j\\\
&v_{j,i},&&i>j\geq 1,\end{aligned}\right.$ (26)
$\mathbf{U}_{i,j}=\sum_{u=1}^{i}\mathbf{A}^{i-u}\mathbf{V}_{w}(\mathbf{A}^{j-u})^{\top},1\leq
i\leq j,$ (27)

and

$\mathbf{C}_{i}=\begin{bmatrix}\mathbf{0}_{(l+1-i)\times K}\\\
\mathbf{b}^{\top}\mathbf{A}^{0}\mathbf{V}_{w}\\\
\mathbf{b}^{\top}\mathbf{A}^{1}\mathbf{V}_{w}\\\ \vdots\\\
\mathbf{b}^{\top}\mathbf{A}^{i-1}\mathbf{V}_{w}\end{bmatrix}.$ (28)

###  IV-B Tx-Scaling Factor Design

Given a low-complexity AirComp policy (13), we will optimize the Tx-scaling
factor $\mathbf{b}$ under a power constraint to minimize the computation MSE
in (22). Although the computation MSE in the current form is an implicit
function of $\mathbf{b}$, from the definitions of $\mathbf{M}$,
$\mathbf{V}_{c}$ and $\mathbf{C}_{i}$ in (18), (LABEL:eq:V_c) and (28),
respectively, it is not difficult to see that (22) is a quadratic function of
$\mathbf{b}$. In the following, we will simplify each term of (22) with
respect to $\mathbf{b}$.

For the first term, we have

$\displaystyle\mathbf{g}^{\top}\mathbf{M}\mathbf{V}_{x}\mathbf{M}^{\top}\mathbf{g}$
$\displaystyle=\sum_{i=0}^{l}\sum_{j=0}^{l}g(i)g(j)\mathbf{b}^{\top}\mathbf{A}^{i}\mathbf{V}_{x}(\mathbf{A}^{j})^{\top}\mathbf{b}$
(29)
$\displaystyle=\mathbf{b}^{\top}\left(\sum_{i=0}^{l}\sum_{j=0}^{l}\mathbf{G}_{i}\mathbf{A}^{i}\mathbf{V}_{x}(\mathbf{A}^{j})^{\top}\mathbf{G}_{j}\right)\mathbf{b},$

where $\mathbf{G}_{i}$ denotes
$\mathsf{diag}\\{\underbrace{g_{i},\dots,g_{i}}_{K}\\},i\in\\{0,1,\dots,l\\}$.
Similarly, it is easy to have

$\displaystyle\mathbf{g}^{\top}\mathbf{M}\mathbf{V}_{x}(\mathbf{A}^{l})^{\top}\mathbf{1}=\mathbf{b}^{\top}\sum_{i=0}^{l}\mathbf{G}_{i}\mathbf{A}^{i}\mathbf{V}_{x}(\mathbf{A}^{l})^{\top}\mathbf{1},$
(30)

and

$\begin{split}\mathbf{g}^{\top}\mathbf{V}_{c}\mathbf{g}=\sigma_{z}^{2}\mathbf{g}^{\top}\mathbf{g}+\mathbf{b}^{\top}\left(\sum_{i=1}^{l}\sum_{j=1}^{l}\mathbf{G}_{i}\mathbf{U}^{\prime}_{i,j}\mathbf{G}_{j}\right)\mathbf{b},\end{split}$
(31)

where $\mathbf{U}^{\prime}_{i,j}=\mathbf{U}_{i,j}$ and $\mathbf{U}_{j,i}$ for
cases $i\leq j$ and $i>j$, respectively. Further, from the definition of (28),
it can be obtained that

$\displaystyle\mathbf{g}^{\top}\sum_{i=1}^{l}\mathbf{C}_{l+1-i}(\mathbf{A}^{l-i})^{\top}\mathbf{1}\\!=\\!\mathbf{b}^{\top}\\!\sum_{i=1}^{l}\sum_{j=i}^{l}\mathbf{G}_{j}\mathbf{A}^{j-i}\mathbf{V}_{w}(\mathbf{A}^{l-i})^{\top}\mathbf{1}.$
(32)

Taking (29)-(32) into (22), we have

$\displaystyle\widetilde{\mathsf{CMSE}}=\mathbf{b}^{\top}\left(\sum_{i=0}^{l}\sum_{j=0}^{l}\mathbf{G}_{i}\mathbf{A}^{i}\mathbf{V}_{x}(\mathbf{A}^{j})^{\top}\mathbf{G}_{j}\right)\mathbf{b}$
(33)
$\displaystyle-2\mathbf{b}^{\top}\\!\\!\left(\\!\sum_{i=0}^{l}\\!\mathbf{G}_{i}\mathbf{A}^{i}\mathbf{V}_{x}(\mathbf{A}^{l})^{\top}\mathbf{1}\\!+\\!\sum_{i=1}^{l}\sum_{j=i}^{l}\\!\mathbf{G}_{j}\mathbf{A}^{j-i}\mathbf{V}_{w}(\mathbf{A}^{l-i})^{\top}\mathbf{1}\\!\right)$
$\displaystyle+\mathbf{1}^{\top}\mathbf{A}^{l}\mathbf{V}_{x}(\mathbf{A}^{l})^{\top}\mathbf{1}+\sigma_{z}^{2}\mathbf{g}^{\top}\mathbf{g}+\sum_{i=1}^{l}\mathbf{1}^{\top}\mathbf{A}^{l-i}\mathbf{V}_{w}(\mathbf{A}^{l-i})^{\top}\mathbf{1}.$

We assume that each sensor has an average transmit power limit $P_{k},\forall
k\in\mathcal{K}$. Then, the average transmit power constraint is
$\mathsf{E}\left[|b_{k}x_{k,t}|^{2}\right]\leq h^{2}_{k}P_{k},$ with a
quadratic form as

$\mathbf{b}^{\top}\mathbf{P}_{k}\mathbf{b}\leq
P_{k}/\sigma^{2}_{x,k},k\in\mathcal{K},$ (34)

where $\sigma^{2}_{x,k}\triangleq\mathsf{E}[x^{2}_{k,t}]$, which is the $k$th
diagonal element of the covariance matrix $\mathbf{V}_{x}$, and
$\mathbf{P}_{k}$ is an all-zero matrix except for the $k$th diagonal element,
which equals to one.

Considering the power constraint (34), the optimal Tx-scaling factor design
problem is formulated as

$\min_{\mathbf{b}}\qquad\eqref{eq:mse_b},\qquad\text{subject to
}\eqref{eq:power_con}.$ (35)

(35) is a classic QCQP problem, which can be effectively solved using semi-
definite programming (SDP) [16].

###  IV-C Joint Design of AirComp Policy and Tx-Scaling Factor

The computation MSE in (33) is a non-convex function of the AirComp filter
$\mathbf{g}$ and the Tx-scaling factor $\mathbf{b}$. To the best of the our
knowledge, there is no general approach to find the optimal $\mathbf{g}$ and
$\mathbf{b}$ jointly. Therefore, we adopt the _alternating minimization
approach_ by alternately solving problems (21) and (35) for $\mathbf{g}$ and
$\mathbf{b}$, respectively, while fixing the other. Note that alternating
minimization is a widely applicable and empirically successful approach for
optimization problems with different subsets of variables [17]. The
effectiveness of the approach in the AirComp joint design problem will be
verified via numerical results in the following section.

##  V Numerical Results

In this section, we numerically evaluate the computation MSE of the optimal
and low-complexity policies with different number of sensors $K$. Unless
otherwise stated, we set $\mathbf{A}=\alpha\mathbf{I}$ [11], where
$0<\alpha<1$, $\mathbf{V}_{x}=\mathbf{I}$, $\sigma^{2}_{z}=1$ and
$P_{k}=10,\forall k\in\mathcal{K}$ [14]. A large $\alpha$ means that the
process ${\mathbf{x}_{t}}$ is tightly time-correlated. From (2), we have
$\mathbf{V}_{w}=(1-\alpha^{2})\mathbf{I}$. We assume normalized Rayleigh
fading channels for evaluating the AirComp system performance that averages
over $10^{5}$ random channel realizations.

In Fig. 1, we first plot the computation MSE of low-complexity policy with
AirComp filter $\mathbf{g}$ and Tx-scaling factor $\mathbf{b}$ achieved by the
alternating minimization approach in Section IV-C with $50$ iteration rounds.
We see that the computation MSE decreases with the length of the filter
$(l+1)$ as expected, and a larger sensor number $K$ leads to a higher
commutation MSE. It is interesting to see that $l=1$ has achieved the optimal
computation MSE when $\alpha=0.9$, while it is larger than $5$ when
$\alpha=0.99$. So we need a longer filter to achieve the minimum computation
MSE when the process ${\mathbf{x}_{t}}$ is tightly time-correlated. Then, we
plot the computation MSE achieved by the optimal AirComp policy in Theorem 1
and the optimal Tx-scaling factor of the low-complexity policy mentioned
earlier. We see that the performance of the low-complexity policy approaches
to the optimal one when $l$ is large.

Figure 1: The computation MSE of the optimal and low-complexity AirComp
policies.

In Fig. 2, we evaluate the convergence speed of the alternating minimization
approach of the low-complexity AirComp policy with $l=4$ under an equal
channel power gain $\mathbf{h}=\mathbf{s}_{1}=[1,1,\cdots,1]^{\top}$ of
different sensors and an unequal channel power gain
$\mathbf{h}=\mathbf{s}_{2}$, where $\mathbf{s}_{2}$ is generated by $K$ evenly
spaced values between $0.1$ and $1.9$. We see that the proposed algorithm
converges within $30$ iterations in different scenarios. It can be observed
that a larger number of sensors and a higher divergence of channel conditions
of different channels require more iteration rounds to converge.

Figure 2: The computation MSE of the low-complexity AirComp policy with
different iteration rounds.

##  VI Conclusions

In this letter, we have proposed an AirComp system with spatial-and-temporal
correlated sensor signals, and derived the optimal AirComp policy to achieve
the minimum computation MSE in each time step by utilizing the current and the
previously received signals. We have also proposed and optimized a low-
complexity AirComp policy with a linear AirComp filter, based on which the
optimal transmit power of each sensor has been optimized. An alternating
minimization approach has been adopted to jointly design the AirComp filter
and the Tx-scaling factor.

## References

  * [1] M. Goldenbaum and S. Stanczak, “On the channel estimation effort for analog computation over wireless multiple-access channels,” _IEEE Wireless Commun. Lett._ , vol. 3, no. 3, pp. 261–264, June 2014.
  * [2] W. Liu, X. Zang, Y. Li, and B. Vucetic, “Over-the-air computation systems: Optimization, analysis and scaling laws,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 8, pp. 5488–5502, Aug. 2020.
  * [3] G. Zhu and K. Huang, “MIMO over-the-air computation for high-mobility multimodal sensing,” _IEEE Internet Things J._ , vol. 6, no. 4, pp. 6089–6103, Aug 2019.
  * [4] D. Wen, G. Zhu, and K. Huang, “Reduced-dimension design of mimo over-the-air computing for data aggregation in clustered iot networks,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5255–5268, 2019.
  * [5] M. Mohammadi Amiri and D. Gündüz, “Machine learning at the wireless edge: Distributed stochastic gradient descent over-the-air,” _IEEE Trans. Signal Process._ , vol. 68, pp. 2155–2169, 2020.
  * [6] F. Molinari, S. Stanczak, and J. Raisch, “Exploiting the superposition property of wireless communication for average consensus problems in multi-agent systems,” in _Proc. ECC_ , 2018, pp. 1766–1772.
  * [7] M. Goldenbaum and S. Stanczak, “Robust analog function computation via wireless multiple-access channels,” _IEEE Trans. Commun._ , vol. 61, no. 9, pp. 3863–3877, Sep. 2013.
  * [8] M. Goldenbaum, H. Boche, and S. Stańczak, “Nomographic gossiping for f-consensus,” in _Proc. IEEE WiOpt_ , May 2012, pp. 130–137.
  * [9] O. Abari, H. Rahul, and D. Katabi, “Over-the-air function computation in sensor networks,” _arXiv preprint_ , 2016. [Online]. Available: https://arxiv.org/pdf/1612.02307.pdf
  * [10] K. Huang, W. Liu, M. Shirvanimoghaddam, Y. Li, and B. Vucetic, “Real-time remote estimation with hybrid ARQ in wireless networked control,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 5, pp. 3490–3504, May 2020\.
  * [11] A. S. Leong, S. Dey, G. N. Nair, and P. Sharma, “Power allocation for outage minimization in state estimation over fading channels,” _IEEE Trans. Signal Process._ , vol. 59, no. 7, pp. 3382–3397, Jul. 2011.
  * [12] T. Kailath, A. H. Sayed, and B. Hassibi, _Linear estimation_. Prentice Hall, 2000.
  * [13] W. Liu, P. Popovski, Y. Li, and B. Vucetic, “Wireless networked control systems with coding-free data transmission for Industrial IoT,” _IEEE Internet Things J._ , vol. 7, no. 3, pp. 1788–1801, Mar. 2020.
  * [14] X. Zang, W. Liu, Y. Li, and B. Vucetic, “Over-the-air computation systems: Optimal design with sum-power constraint,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 9, pp. 1524–1528, Sep. 2020.
  * [15] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [16] Z. Luo, W. Ma, A. M. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” _IEEE Signal Process. Mag._ , vol. 27, no. 3, pp. 20–34, May 2010.
  * [17] X. Yu, J. Shen, J. Zhang, and K. B. Letaief, “Alternating minimization algorithms for hybrid precoding in millimeter wave mimo systems,” _IEEE JSTSP_ , vol. 10, no. 3, pp. 485–500, Apr. 2016.