Algorithms for audio inpainting based on probabilistic nonnegative matrix factorization

The work was supported by the Czech Science Foundation (GAČR) Project No. 20-29009S, the French ANITI Project No. ANR-19-PI3A-0004 and the European Research Council (ERC) Project FACTORY No. 6681839.
Part of the work was realized during the stay of O. Mokrý at IRIT, co-financed from the Erasmus+ mobility program.
The authors would like to thank Pavel Rajmic and Pierre-Hugo Vial for their inputs concerning both the development of the methods and the preparation of the manuscript.

but]Ondřej Mokrýmycorrespondingauthor
[mycorrespondingauthor]Corresponding author

inria]Paul Magron

isae]Thomas Oberlin

irit]Cédric Févotte

[but]Brno University of Technology, Faculty of Electrical Engineering and Communications, Department of Telecommunications, Technická 12, 616 00 Brno, Czech Republic
[inria]Université de Lorraine, CNRS, Inria, LORIA, F-54000 Nancy, France
[isae]ISAE-SUPAERO, Université de Toulouse, France
[irit]IRIT, Université de Toulouse, CNRS, Toulouse, France

Audio inpainting, i.e., the task of restoring missing or occluded audio signal samples, usually relies on sparse representations or autoregressive modeling.
In this paper, we propose to structure the spectrogram with nonnegative matrix factorization (NMF) in a probabilistic framework. First, we treat the missing samples as latent variables, and derive two expectation–maximization algorithms for estimating the parameters of the model, depending on whether we formulate the problem in the time- or time-frequency domain. Then, we treat the missing samples as parameters, and we address this novel problem by deriving an alternating minimization scheme.
We assess the potential of these algorithms for the task of restoring short- to middle-length gaps in music signals. Experiments reveal great convergence properties of the proposed methods, as well as competitive performance when compared to state-of-the-art audio inpainting techniques.

alternating minimization,
audio inpainting,
nonnegative matrix factorization

§ INTRODUCTION

Audio inpainting <cit.> is an inverse problem aiming at restoring audio signals degraded by sample loss.
Such a problem typically occurs as a result of packet loss during transmission (packet loss concealment <cit.>) or in digitization of physically degraded media.
Inpainting can also be used to restore signal samples subject to a degradation so heavy that the information about the samples can be considered lost.
Formally, let $\sig\in\CC^\Tdim$ denote the original time-domain signal prior to the degradation (we consider complex-valued signals for the sake of generality).
The goal of inpainting is to estimate $\sigh\in\CC^\Tdim$, given an incomplete observation of $\sig$.
This problem is ill-posed because of the missing samples, and even the observed samples are prone to some measurement error or noise.
Restoring the signal thus requires some
assumptions on the structure of
the original signal,
in order to guide the estimation towards the most desirable solution.

One of the earliest, yet most successful approaches to audio inpainting is to assume the underlying autoregressive (AR) nature of clean audio signals: Janssen's iterative method <cit.> is still among state-of-the-art methods to date.
Great performance is also achieved by Etter's extrapolation-based technique <cit.>, however, it is limited by the need for a sufficiently long clean context to reconstruct each individual gap.
Recently, the class of sparsity-based methods methods has emerged <cit.>.
Generally speaking, these methods solve a regularized inverse problem, where the solution is assumed to have a sparse time-frequency (TF) spectrum, while fitting the observed temporal samples.

A disadvantage of most methods is the local nature of the regularizers.
For example, sparsity-based methods are effective for inpainting of short- to medium-length gaps, typically up to 50 ms, or for restoring randomly subsampled signals <cit.>.
Even for gaps of length of tens of milliseconds, there is a need for strong regularization that not only exploits the local TF sparsity but also the global structure of audio signals.
A step towards using global properties of audio signals is the so-called social sparsity <cit.>, where the significant TF coefficients are expected to occur in particular transient or temporal patterns, or recent approaches based on generative deep neural networks on medium gaps <cit.> and long gaps in the range of seconds <cit.>.

In the present paper, we focus on audio inpainting methods for medium gaps (in the range of tens of milliseconds) that leverage
the low-rank structure of audio signals in the TF domain. Among low-rank models, nonnegative matrix factorization (NMF)
has been intensively used for analysis and decomposition, both in machine learning and signal processing <cit.>.
Concerning audio, NMF can be used to decompose the signal's power or magnitude spectrogram
as the product of nonnegative matrices $\W\H$, where $$ is a dictionary of spectral patterns, and $$̋ contains the temporal activations of these patterns.
Such a decomposition provides a semantically reasonable generative model for audio signals, which is one of the reasons why NMF is among the classical approaches to source separation <cit.>.
It also benefits from being an unsupervised and interpretable method while being computationally cheap.
NMF has also been used as a prior for the restoration of degraded signals <cit.>, thus being successfully applied to audio declipping <cit.>.

The main idea of NMF-based signal reconstruction can be summarized as building an estimation problem where the NMF is used as a prior, which requires to estimate parameters given the incomplete data, usually via the expectation–maximization (EM) algorithm <cit.>.
Formally, the structure
of the spectrogram is introduced via the probabilistic Gaussian composite model <cit.>, which results in applying NMF with the Itakura–Saito divergence.
However, the estimation of the parameters – the matrices $\W$ and $\H$ -- is a non-trivial problem which can be approached in different ways.
One possibility is to derive a generalized EM algorithm to estimate the factorization of the power spectrogram of the clean signal in the maximum likelihood (ML) sense.
This has been previously presented by Bilen, Ozerov and Pérez \cite{BilenOzerovPerez2015:declipping.via.NMF,Bilen2018:NTF_audio_inverse_problems} with application to audio declipping, where it reaches state-of-the-art performance.
However, the method is known to be computationally very demanding \cite[Sec.\,V.D]{ZaviskaRajmicOzerovRencker2021:Declipping.Survey}, and it represents only one of several possible approaches to treating the missing samples (specifically, they are treated as latent variables).
Furthermore, it has not been compared to the state-of-the-art methods in the audio inpainting setting.

Drawing on the work of Bilen et al., we provide a novel generalization of their EM-based approach, where the missing samples are formally treated as latent variables.
We formulate the estimation problem in both time and TF domains and therefore derive two algorithms, among which one is novel.
As a novel approach, we also propose to treat the missing samples as parameters. This leads to a new estimation problem, for which we derive an alternating minimization (AM) scheme.
Even though we built upon the application in audio inpainting, the core of the work is the conceptual development of novel approaches to the NMF-based modeling in signal restoration.
We conduct experiments on the task of restoring
\edt{short to middle gaps} %short- to middle-length gaps
in music signals \edt{(in particular gaps up to 80\,ms)}. These reveal great convergence properties of the proposed methods, as well as competitive performance when compared to state-of-the-art audio inpainting techniques.

The rest of the paper is organized as follows.
We formulate audio inpainting with NMF as an estimation problem in Section \ref{sec:formulation}.
In Section~\ref{sec:mle.1} we
review and extend the approach by Bilen et al.\ \cite{BilenOzerovPerez2015:declipping.via.NMF}.
In Section~\ref{sec:mle.2} we consider the missing samples as parameters, for which we derive a novel AM estimation approach.
Section \ref{sec:experiments} is devoted to experiments and evaluation of the proposed methods.
Finally, Section \ref{sec:conclusion} concludes the paper.

\section{Problem formulation}
\label{sec:formulation}

For the whole derivation of the methods, it is convenient to divide the signal into windowed temporal frames $_n∈^, n=1,…,N$, potentially arranged in a matrix
$= [_1,…,_N]$.
In the case of inpainting, we observe the samples $_n = _n_n$ in each frame $n$ and we aim at obtaining an estimate $_n$ of the whole frame given these observed samples.
The matrix $_n$ is constructed from the identity matrix by omitting the rows corresponding to the missing samples, thus the multiplication with $_n$ shortens the vector, selecting only the observed samples.
Note that in the present work, the \emph{indices} of the missing samples are assumed to be known, thus the matrices $_n, n=1,…,N$ are known.
The whole signal estimate $∈^$ is then obtained by folding all $_n$ together by the overlap-add procedure,

To propose a probabilistic formulation of the inpainting problem, it is necessary to postulate a statistical model for the data at hand.
Since we aim at promoting the low-rank structure of the TF coefficients $= [_1,…,_N] = [s_fn] ∈^$ of the audio signal, we first define the following synthesis model:
\begin{equation}
\fsig_n = \syn\coef_n,\quad n = 1,\dots,N,
\label{eq:reconstruction.operator}
\end{equation}
or, in matrix form, $= $, where $∈^×$ is a linear reconstruction operator.
The matrix $$ typically represents the inverse discrete Fourier transform (DFT).
If the temporal frames are weighted by a windowing function, $$ represents the short-time Fourier transform (STFT) coefficients of the whole signal $$.
We will discuss the effect of particular choices of $$ later in \ref{ssec:equivalence}.

The low-rank structure of the TF coefficients can be formalized within the following assumptions.

\begin{assumption}[Gaussian coefficients]
\label{ass:gaussian}
The time-frequency coefficients
are treated as conditionally mutually independent,\footnote{Note that temporal Markov NMF models such as in \cite{Smaragdis2014:Static.dnymic.source.separation.NMF} could readily be considered but we use independence as a working assumption for ease of presentation.} and each coefficient follows a complex circular zero-mean Gaussian distribution:
\begin{equation}
s_{fn}\sim\mathcal{N}(0, v_{fn}).
\end{equation}
\end{assumption}

Due to the assumed independence of the individual TF coefficients, we can rewrite the distribution of the spectrum of the $n$-th frame as:
\begin{equation}
\coef_n \sim \mathcal{N}(\matr{0}, \D_n),\quad\D_n = \diag \left([v_{fn}]_{f=1,\dots,\Fdim} \right),
\label{eq:distribution.s}
\end{equation}
where the symbol $𝒩$ from now on denotes the multivariate complex Gaussian distribution.

\begin{remark}
\label{rem:dft.independence}
Note that Assumption \ref{ass:gaussian} is stated for the case of a generic linear transform $\syn$ interconnecting the TF coefficients and temporal samples.
In the common case where the TF coefficients are computed from the temporal samples by the DFT (i.e., $\syn$ represents the inverse DFT operator), and the temporal signal is real-valued, the independence assumption needs to be relaxed:
The TF coefficients of real-valued signals are independent only in half of the frequency spectrum, the other half being determined via the Hermitian property of the transform.
However, for the sake of generality, we will assume no particular structure of $\syn$ throughout the article.
\end{remark}

\begin{assumption}[NMF structure of the variances]
\label{ass:nmf}
The variance matrix $\matr{V} = [v_{fn}]$ has the low-rank NMF structure:
\begin{equation}
v_{fn} = \sum_{k=1}^{K} w_{fk}h_{kn},
\end{equation}
where $K$ is small and all parameters are nonnegative.
This model amounts to $\matr{V} = \W\H$ with $\W$ and $\H$ being $\Fdim \times K$ and $K \times N$ nonnegative matrices, respectively \cite[Sec.\,2.2]{BilenOzerovPerez2015:declipping.via.NMF}.
\end{assumption}

To estimate the parameters $,$̋ of the variance matrix, given the observed samples $\FOBS = \{ \fobs_1,\dots, \fobs_N \}$, we employ maximum likelihood (ML) estimation.
The audio inpainting per se, i.e., the computation of the missing samples, is then performed explicitly given the estimated parameters and the reliable samples in a way similar to Wiener filtering.

However, there is not a unique way to formulate the ML problem.
The following sections present two approaches that differ in whether the missing samples are treated as latent variables (Section <ref>) or explicitly treated as parameters
(Section <ref>).
Note that subsection <ref>
reformulates the method proposed by Bilen et al. in <cit.>, whereas the rest of Section <ref> and the whole Section <ref> are new contributions.

§ ML ESTIMATION BY TREATING THE MISSING SAMPLES AS LATENT VARIABLES

We first present the ML formulation, where the goal is to estimate the parameters $\W,\H$ of the distribution of the restored signal in the TF domain, by minimizing the negative log-likelihood of the observed samples $= {_1,…,_N}$, as presented by Bilen et al.\ in \cite{Bilen2018:NTF_audio_inverse_problems}.
The problem can be formalized as
\begin{equation}
\hat\W, \hat\H = \argmin_{\W,\H} -\log p(\,\FOBS\mid\W,\H\,).
\label{eq:mle.1}
\end{equation}

This expression can be broken down for a single frame $n$, where the probability $p$
is given by the distribution of $_n$ from Eq.\,\eqref{eq:distribution.s} and by  the linear observation model as
\begin{equation}
\fobs_n = \Mr_n\fsig_n = \Mr_n\syn\coef_n\sim\mathcal{N}(\matr{0},\Mr_n\syn \D_n \adjoint{\syn}\transp{\Mr}_n),\quad n = 1,\dots,N,
\label{eq:distribution.fobs}
\end{equation}
where the dependence on $,$̋ is contained within the definition of the matrices $\D_n$.

We resort to the EM algorithm <cit.>, in line with Bilen et al. <cit.>.
However, we broaden their work by considering two different settings of the EM algorithm, depending on the domain of the complete data to be estimated.

First, we briefly describe the EM algorithm for the problem (<ref>).
The setting is that the incomplete data corresponds to the observed reliable signal, i.e., $\FOBS = \{\fobs_1, \ldots, \fobs_N\}$ in the framed time domain.
The complete data corresponds to the TF spectrum $\COEF = [\coef_1, \ldots, \coef_N] \in\CC^\TFdim$ of the original signal (which is reflected by the abbreviation ).
Finally, the parameters to be estimated are $\btheta = \{ \W, \H \}$, and $\tilde{\btheta}$ is the current value of the parameters.

Using this setting, the EM algorithm aims at minimizing the functional:

\begin{equation}
Q(\btheta,\tilde{\btheta}) = - \int \log p(\COEF \mid \btheta) p(\COEF\mid\FOBS, \tilde{\btheta}) \,\mathrm{d}\COEF
\label{eq:em.functional}
\end{equation}

by iterating two steps:

* E-step: compute $Q(\btheta,\tilde{\btheta})$,

* M-step: update $\tilde{\btheta}$
by minimizing $Q(\btheta,\tilde{\btheta})$ with respect to $\btheta$.

For the E-step,
we obtain

\begin{equation}
p(\COEF \mid \FOBS, \tilde{\btheta}) = \prod_{n=1}^N p(\coef_n \mid \fobs_n, \tilde{\btheta}) = \prod_{n=1}^N \mathcal{N}(\coef_n\mid\coefh_n, \bSigmah_n),
\end{equation}

\begin{align}
% mean
\coefh_n &= \D_n \adjoint{\syn}\transp{\Mr}_n
\inv{\left(\Mr_n\syn \D_n \adjoint{\syn}\transp{\Mr}_n\right)}
\fobs_n,
\label{eq:e-step.1:mean}
\\
% covariance matrix
\bSigmah_n  &=
\D_n - \D_n \adjoint{\syn}\transp{\Mr}_n
\inv{\left(\Mr_n\syn \D_n \adjoint{\syn}\transp{\Mr}_n\right)}
\Mr_n\syn \D_n,
\label{eq:e-step.1:covariance}
\end{align}

and the matrices $\D_n$ are computed using the current value of the parameters $\tilde{\btheta}$.
The formulas (<ref>) can be derived from the assumed distribution of $\coef_n$ in (<ref>) and from the linear observation model $\fobs_n = (\Mr_n\syn)\coef_n$ (see e.g. <cit.>). Note that these formulas can be seen as Wiener filtering given estimated posterior covariance <cit.> to recover the missing samples.

Now, from (<ref>) and from the independence assumption, it holds that

\begin{equation}
p(\COEF \mid \btheta) = \prod_{n=1}^N \mathcal{N}(\coef_n \mid \matr{0}, \D_n),\quad\D_n = \diag \left([v_{fn}]_{f=1,\dots,\Fdim} \right).
\end{equation}

The M-step, i.e., the minimization of (<ref>) with respect to the parameters $\btheta = \{\W,\H\}$ is equivalent to the minimization of
the Itakura–Saito divergence $D_\mathrm{IS}(\P\mid\W\H)$ <cit.>, where the divergence is defined for matrices $\matr{A} = [a_{ij}],\matr{B} = [b_{ij}]$ as

\begin{equation}
D_\mathrm{IS}(\matr{A}\mid\matr{B}) = \sum_{i,j} d_\mathrm{IS}(a_{ij} \mid b_{ij}) = \sum_{i,j} \left( \frac{a_{ij}}{b_{ij}} - \log\frac{a_{ij}}{b_{ij}} - 1 \right),
\label{eq:is}
\end{equation}

and the matrix $\P = [p_{fn}]$ is the posterior power spectrum given by:

\begin{equation}\label{eq:posterior.power}
p_{fn} = \mathbb{E}\left( \abs{s_{fn}}^2 \mid \fobs_n, \W, \H \right) = \abs{(\coefh_n)_f}^2 + (\bSigmah_n)_{ff}.
\end{equation}

The M-step can thus be performed
by applying (and iterating) the following multiplicative rules <cit.>:

\begin{align}
\W &\leftarrow \W \odot \frac{%
\left((\W\H)^{\odot[-2]} \odot \P\right) \transp{\H}
(\W\H)^{\odot[-1]} \transp{\H}
\label{eq:wh:w}
\\
\H &\leftarrow \H \odot \frac{%
\transp{\W}\left((\W\H)^{\odot[-2]} \odot \P\right)
\transp{\W}(\W\H)^{\odot[-1]}
\label{eq:wh:h}
\end{align}

where $\frac{\matr{A}}{\matr{B}}$ denotes the matrix $\matr{A}\odot\matr{B}^{\odot[-1]}$ and the symbol $\odot$ is used to denote entry-wise multiplication or power.
Note that in the algorithm, the $\H$ update \eqref{eq:wh:h} uses the already updated value of $$ from \eqref{eq:wh:w}.
In practice, the two updates are followed by a normalization step:
The columns of $$ are scaled to have unit norm and the rows of $$̋ are inversely scaled by the same factor, such that the product $\W\H$ does not change.

The whole \EMTF{} algorithm is summarized in Alg.\,\ref{alg:emtf}.
Its output is the estimate of the complete data, i.e., the full TF spectrum $∈^$, together with the estimate of the parameters $$ and $$̋.
Then, the framed time-domain signal $\FSIGH$ is synthesized from $\COEFH$ using the operator $\syn$. Finally, the whole signal estimate $\sigh$ is obtained from $\FSIGH$ using the overlap-add procedure.

As an alternative to , we build an EM algorithm that uses a complete dataset in the time domain (hence the abbreviation ).
As a result, we directly get the posterior distribution of the missing data as a side-product of the algorithm.

This novel algorithm addresses the original problem (<ref>) by minimizing the functional:

\begin{equation}
Q(\btheta,\tilde{\btheta}) = - \int \log p(\FSIG \mid \btheta) p(\FMISS\mid\FOBS, \tilde{\btheta}) \,\mathrm{d}\FMISS,
\end{equation}

where $\FMISS = \{\fmiss_1,\dots,\fmiss_N \}$ represents the missing samples.
These are identified as ${\fmiss_n = \Mm_n\fsig_n}$,
where $\Mm_n$ is the complementary selection matrix to $\Mr_n$, i.e., $\Mm_n$ selects the unrealiable (missing) samples of frame $n$.

It is clear that in some cases, this approach will be equivalent to , such as in the case of one-to-one correspondence between the temporal frame samples and their frequency coefficients.
However, this part discusses a general situation with no other assumptions on the transforms involved, and particular cases will be examined later in Section <ref>.

To derive the steps of the EM algorithm, observe that the relation $\fsig_n = \syn\coef_n$, together with Assumption <ref> (ass:gaussian), directly leads to

\begin{equation}\label{eq:distribution.x}
p(\FSIG\mid\btheta) = \prod_{n=1}^N \mathcal{N}\left(\fsig_n \mid \matr{0}, \syn\D_n\adjoint{\syn} \right).
\end{equation}

Using the selection matrices $\Mr_n$ and $\Mm_n$, we derive the distribution

\begin{equation}
\begin{bmatrix}
\Mr_n\fsig_n \\ \Mm_n\fsig_n
\end{bmatrix}
\begin{bmatrix}
\fobs_n \\ \fmiss_n
\end{bmatrix}
\sim
\mathcal{N}\left(\matr{0},
\begin{bmatrix}
\Mr_n\syn\D_n\adjoint{\syn}\transp{\Mr}_n & \Mr_n\syn\D_n\adjoint{\syn}\transp{\Mm}_n \\[1ex]
\Mm_n\syn\D_n\adjoint{\syn}\transp{\Mr}_n & \Mm_n\syn\D_n\adjoint{\syn}\transp{\Mm}_n
\end{bmatrix}
\begin{bmatrix}
\C_{11} & \C_{12} \\[1ex]
\C_{21} & \C_{22}
\end{bmatrix}
\right),
\end{equation}

from where it is straightforward to derive that $\fmiss_n$, given $\fobs_n$ and the model parameters, follows the distribution
$\mathcal{N}\left(\C_{21}\inv{\C_{11}}\fobs_n,\,\C_{22} - \C_{21}\inv{\C_{11}}\C_{12} \right)$.
Using the results from Eq. (<ref>), we can simply express
$\C_{21}\inv{\C_{11}}\fobs_n = \Mm_n\syn\coefh_n$
$\C_{22} - \C_{21}\inv{\C_{11}}\C_{12} = \Mm_n\syn\bSigmah_n\adjoint{\syn}\transp{\Mm}_n$. This yields:

\begin{equation}\label{eq:em2.from.em1}
p(\FMISS \mid \FOBS, \tilde{\btheta}) = \prod_{n=1}^N p(\fmiss_n \mid \fobs_n, \tilde{\btheta}) = \prod_{n=1}^N
\mathcal{N}(\fmiss_n\mid \Mm_n\syn\coefh_n,\, \Mm_n\syn\bSigmah_n\adjoint{\syn}\transp{\Mm}_n).
\end{equation}

The crucial part of the algorithm is the
ML estimation of the parameters $\W,\H$ in the M-step, given the posterior distribution of the missing temporal samples.
Even though the closed form of $Q(,)$ is available due to the expressions in Eq.\,\eqref{eq:distribution.x} and \eqref{eq:em2.from.em1}, it is expensive to compute and optimize directly.
Thus, we proceed to re-estimate the spectrum corresponding to the signal estimated by Eq.\,\eqref{eq:em2.from.em1} and update $$ and $$̋ as the factorization of this spectrum.

To do this, we introduce an analysis operator, represented by the matrix $\ana\in\CC^{\Fdim\times\Wdim}$, associated to the synthesis operator $\syn$.
So far, the only assumption about the analysis operator is linearity – more details are given below.
Using $\ana$, it is straightforward to derive the posterior distribution of the TF coefficients $\COEF^{\textup{alt}} = \ana\FSIG$ associated to the posterior time-domain samples:

\begin{equation}
p(\COEF^{\textup{alt}} \mid \FOBS, \tilde{\btheta}) = \prod_{n=1}^N p(\coef_n^{\textup{alt}} \mid \fobs_n, \tilde{\btheta}) = \prod_{n=1}^N
\mathcal{N}(\coef_n^{\textup{alt}} \mid \coefh_n^{\textup{alt}},\, \bSigmah_n^{\textup{alt}})
\end{equation}

\begin{equation}\label{eq:e-step.2}
% mean
\coefh_n^{\textup{alt}}	= \ana\syn\coefh_n,
\quad
% covariance matrix
\bSigmah_n^{\textup{alt}} = \ana\syn\bSigmah_n\adjoint{\syn}\adjoint{\ana}.
\end{equation}

Finally, the M-step is equivalent to the M-step defined by the updates in Eq. (<ref>), with the alternative posterior power spectrum (computed from $\coefh_n^{\textup{alt}}$ and $\bSigmah_n^{\textup{alt}}$).
The whole algorithm is summarized in Alg. <ref>.

§.§ On the relations between  and

As previously mentioned, there are some natural choices of the pair $\{ \syn,\ana \}$ which result in the equivalence of the algorithms  and .
Several alternatives are discussed in what follows.

$\syn$ is invertible, $\ana = \inv{\syn}$.

It follows from (<ref>) that $\coefh_n = \coefh_n^\textup{alt}$ and $\bSigmah_n = \bSigmah_n^\textup{alt}$, therefore  and  are identical algorithms.
One noticeable special case is when $\syn$ is unitary, i.e., $\ana = \inv{\syn} = \adjoint{\syn}$.
A popular example is the case of a properly scaled DFT realized formally by multiplication with the unitary matrix $\ana\in\CC^{\Wdim\times\Wdim}$.

$\syn$ is the synthesis operator of a tight frame <cit.>, $\Fdim > \Wdim$, $\ana = \adjoint{\syn}$, $\syn\ana = \Id$, where $\Id$ denotes the identity matrix of appropriate size.

The two algorithms are no longer equivalent, since $\ana\syn$ represents the projection operator on the range space of $\ana$, which is in general different from identity.
For instance, this is the case of a redundant DFT, e.g., with $\Fdim = 2\Wdim$ (twice more frequency bins than the frame length).
In practice, this can be implemented by zero-padding the signal to twice its length and then computing the DFT.
The backward transform is the inverse DFT, followed by cropping the result to the original length.

* $\syn$ is the analysis operator of a tight frame, $\Fdim < \Wdim$, $\ana = \adjoint{\syn}$, i.e., $\ana\syn = \Id$.

In this case, we do not have enough frequency coefficients to reconstruct any signal in the framed time domain.
This means that the time-domain solution (in each frame) is restricted to the range space of $\syn$.
However, the estimation  and  are once again equivalent in this case.

* $\syn$ is arbitrary, $\ana = \pinv{\syn}$.

In this case, the matrix $\ana\syn = \pinv{\syn}\syn$ used in (<ref>) represents the orthogonal projection onto the range space of $\adjoint{\syn}$ (which equals the orthogonal complement of the kernel of $\syn$).
This is in general different from identity, unless the range space of $\adjoint{\syn}$ is the whole coefficient space $\CC^\Fdim$.

* $\ana$ is arbitrary, $\syn = \pinv{\ana}$.

Similarly to the previous option, $\ana\syn = \ana\pinv{\ana}$ represents the orthogonal projection onto the range space of $\ana$ (which equals the orthogonal complement of the kernel of $\adjoint{\ana}$).
This is in general different from identity, unless the range space of $\ana$ is the whole coefficient space $\CC^\Fdim$.

Even though the list is far from being exhaustive, it illustrates that there are commonly used settings (e.g., the redundant DFT used in some sparsity-based reconstruction algorithms <cit.>) where  and  are not just conceptually but also practically different.
This will be further detailed in the experiments in Section <ref>.

§ ML ESTIMATION BY TREATING THE MISSING SAMPLES AS PARAMETERS: THE  ALGORITHM

As a novel approach, we propose to treat the missing samples as parameters and include them explicitly into the estimation problem, which results in:

\begin{equation}
\hat\W,\hat\H,\FMISSh = \argmin_{\W,\H,\FMISS} -\log p\big( \FOBS, \FMISS \mid \W, \H \big).
\label{eq:mle.2}
\end{equation}

Note that the objective of this novel estimator is a function different from the likelihood in (<ref>), since the parameter space is extended by the inclusion of the missing samples in the problem.
We propose to approach it via AM. This approach consists of two steps – minimization of (<ref>) with respect to the missing samples (signal update) and minimization with respect to the NMF parameters (model update).

§.§ Signal update

Performing the signal update means minimizing the objective in (<ref>) with respect to $\fmiss_n$ while the current estimates of $\W$, $\H$ are fixed.
This is equivalent to finding the mode of the conditional distribution of $_n$ given $_n,,$̋, which, due to the Gaussian assumption, equals its expectation:

\begin{equation}\label{eq:am.signal.update}
\fmissh_n
= \mathbb{E}\left( \fmiss_n \mid \fobs_n, \W, \H \right).
\end{equation}

Note that this expectation is the same quantity as $\Mm_n\syn\coefh_n$ used in Eq. (<ref>).
Using the expression for $\coefh_n$ from Eq. (<ref>) yields:

\begin{equation}\label{eq:am.signal.update}
\fmissh_n
= \Mm_n\syn\D_n\adjoint{\syn}\transp{\Mr}_n
\inv{\left(\Mr_n\syn\D_n\adjoint{\syn}\transp{\Mr}_n\right)}
\fobs_n.
\end{equation}

The whole signal frame, including both the estimated missing samples $\fmissh_n$ and the observed samples $\fobs_n$, can be merged together as:

\begin{equation}\label{eq:am.signal.update.whole}
\fsigh_n = \syn\D_n\adjoint{\syn}\transp{\Mr}_n
\inv{\left(\Mr_n\syn\D_n\adjoint{\syn}\transp{\Mr}_n\right)}
\fobs_n.
\end{equation}

§.§ NMF parameters update

For the NMF model update, we aim at deriving the computation directly from the optimization problem (<ref>).
Since not only the observed samples but also the (estimated) missing ones are fixed in this step,
this is equivalent to minimizing

\begin{equation}
-\log p\left(
\fsig_n
\mid\W,\H
\right)
= \log{\det{\pi\syn\D_n\adjoint{\syn}}}
+ \transp{\fsig_n}
\inv{\left(\syn\D_n\adjoint{\syn}\right)}
\fsig_n
\end{equation}

with respect to $\W,\H$. To simplify the development of the method, let us pose the following assumption.

\begin{assumption}[invertibility of the synthesis]\label{ass:invertible}
The synthesis operator $\syn$ is invertible and the analysis operator is $\ana = \inv{\syn}$.
In particular, this means that $\syn$ is square, i.e., $\Fdim = \Wdim$.
\end{assumption}

Under assumption \ref{ass:invertible}, we see that $π_n = π^^2_n$ and $(_n) = ()_n$.
The optimization problem then reduces to:
\begin{equation}\label{eq:am.invertible}
\argmin_{\W,\H} \; \log\det{\D_n} + \adjoint{\left(\inv{\syn}\fsig_n\right)}\inv{\D_n}\left(\inv{\syn}\fsig_n\right).
\end{equation}
Now recall that $_n = ([v_fn]_f=1,…, )$, thus we can rewrite the objective function as:
\begin{equation}\label{eq:is-nmf}
\log\prod_{f=1}^{\Fdim} v_{fn} + \sum_{f=1}^{F} \left(\inv{\syn}\fsig_n\right)_f^*\frac{1}{v_{fn}}\left(\inv{\syn}\fsig_n\right)_f
\sum_{f=1}^{\Fdim} \log v_{fn} + \sum_{f=1}^{\Fdim}\frac{\abs{\left(\inv{\syn}\fsig_n\right)_f}^2}{v_{fn}}.
\end{equation}
It is now straightforward to show that the minimization of \eqref{eq:is-nmf} is equivalent to the minimization of the Itakura--Saito divergence:
\begin{equation}
\argmin_{v_{fn}} \sum_{f = 1}^{F} d_\textup{IS}\left(\abs{\left(\inv{\syn}\fsig_n\right)_f}^2\mid v_{fn}\right),\quad v_{fn} = \sum_{k} w_{fk}h_{kn}.
\end{equation}
Taking into account all the frames finally leads to the desired result that $,$̋ are obtained by minimizing $D_\mathrm{IS}(\P\mid\W\H)$ where $p_{fn} = \abs{(\inv{\syn}\fsigh_n)_f}^2$ and $\fsigh_n$ is the signal estimate from Eq. (<ref>).
The procedure is summarized in Alg. <ref>.

A simple heuristic possibility for the case of non-invertible $\syn$ is to compute the spectrum of $\fsigh_n$, defined in Eq. (<ref>), as

\begin{equation}
\coefh_n = \ana\fsigh_n = \ana\syn\D_n\adjoint{\syn}\transp{\Mr}_n
\inv{\left(\Mr_n\syn\D_n\adjoint{\syn}\transp{\Mr}_n\right)}
\fobs_n
\end{equation}

and the power spectrogram $p_{fn} = \abs{(\coefh_n)_f}^2$.
Then, we apply the multiplicative rules to minimize $D_\mathrm{IS}(\P\mid\W\H)$.
However, this approach is not justified by the minimization of (<ref>) with respect to $\W,\H$.
The problem is that if we cannot compute the inversion $(_n)$ as $A_nB$ for some matrices $A,B$, we cannot separate the individual diagonal entries of $_n$ to fit it to the IS-NMF problem.
\end{remark}

\begin{remark}
Note that in the setting imposed in Assumption \ref{ass:invertible} (\nameref{ass:invertible}), \EMTF{} is equivalent to \EMT, but the alternating minimization produces a different algorithm, because we do not include any covariance matrix in the power spectra.
Thus, \AM{} might provide a different inpainting solution, as demonstrated by the numerical experiments in Section \ref{sec:experiments}.
\end{remark}

\begin{remark}[computational complexity]\label{rem:computational.complexity}
It is intricate to express the computational complexity of the algorithms \EMTF, \EMT{} and \AM.
However, several claims are evident:
\begin{enumerate}
\item In \EMTF{} on line \ref{alg:sigma} of Alg.\,\ref{alg:emtf}, we do not in fact need to compute the whole matrix $\bSigmah_n$, but solely its diagonal entries.
Thus, if the complicated term $\D_n\adjoint{\syn}\transp{\Mr}_n\inv{\left(\Mr_n\syn \D_n \adjoint{\syn}\transp{\Mr}_n\right)}$ is computed and saved on line \ref{alg:coef} (e.g., using Matlab's \href{https://www.mathworks.com/help/matlab/ref/mrdivide.html}{mrdivide} without the need to perform the matrix inversion), then no more matrix multiplications are needed on line \ref{alg:sigma} since the diagonal entries can be extracted by computing only $F$ scalar products.
\item In \AM{}, the matrix $\D_n\adjoint{\syn}\transp{\Mr}_n\inv{\left(\Mr_n\syn \D_n \adjoint{\syn}\transp{\Mr}_n\right)}$ does not need to be calculated at all -- we can first compute $\inv{\left(\Mr_n\syn \D_n \adjoint{\syn}\transp{\Mr}_n\right)}\fobs_n$ efficiently
and then multiply the resulting vector with the matrix $\D_n\adjoint{\syn}\transp{\Mr}_n$.
The difference to \EMTF{} is that the inversion operates with a vector, not a matrix, which makes \AM{} less computationally demanding per iteration than \EMTF{}.
\item Although a similar strategy as in \EMTF{} can be applied to \EMT, computing even the diagonal entries of $\bSigmah_n^{\textup{alt}}$ as derived in Eq.\,\eqref{eq:e-step.2} is more complicated due to the multiplication with $\ana\syn$.
As a result, \EMT{} has computationally the most expensive iteration from the three algorithms, even in settings when the operations can be implemented efficiently using FFT.
\end{enumerate}
\end{remark}

\section{Experiments}
\label{sec:experiments}

In this section, we evaluate the performance of the proposed estimators for the task of inpainting noise-less musical recordings. The implementation
\edtlang{was} %is
done in Matlab with the use of the LTFAT toolbox \cite{LTFAT}. For the sake of reproducibility, the source code is published online.\footnote{\url{https://github.com/ondrejmokry/InpaintingNMF}}

\edt{We first present}
% The experiments focus on audio inpainting of compact gaps, which is a challenging setup. However, we
% \edtlang{started} %start
% with
\edt{an experiment testing the restoration of missing individual samples} %a preliminary experiment dedicated to the restoration of missing samples
at random location \edt{(see e.g.\ \cite{Lieb2018:Audio.Inpainting} for similar use case)}.
\edt{This setup served as a demonstrative example
for the comparison of the estimators on a single signal only,}
\edtlang{was} %is
mainly motivated by the high computational cost of the algorithms (see Section \ref{sssec:random.missing.samples}).
\edt{Then we considered a more challenging setting of compact short to medium gaps.}
\edt{We evaluated the estimators against each other in the gap-filling experiment (Section \ref{sssec:missing.gaps}),
and finally we compared them to state-of-the-art baselines (Section \ref{ssec:comparison.soa}).}

% PARTLY MOVED TO SEC. 5.2.2
% For the gap-filling experiment, we
% \edtlang{used} %use
% the set of 10 musical recordings from the EBU SQAM dataset \cite{EBUSQAM,EBU_SQAM_Manual}, sampled at 44.1\,kHz and shortened to 7 seconds, as used commonly in recent related publications~\cite{MokryRajmic2020:Inpainting.revisited,TaubockRajbamshiBalasz2021:SPAINMOD}. The proposed estimators
% \edtlang{were} %are
% first evaluated against each other (Section \ref{sssec:missing.gaps}), and then compared to state-of-the-art baselines (Section \ref{ssec:comparison.soa}).

\subsection{\edt{Evaluation metrics}}
\label{ssec:evaluation.metrics}

As a measure of audio quality, we
\edtlang{considered} %consider
the %commonly used
signal-to-noise ratio (SNR) defined as \cite{Adler2012:Audio.inpainting}
\begin{equation}
\textup{SNR}(\sig,\sigh) = 10\log_{10}\frac{\norm{\sig}^2}{\norm{\sig-\sigh}^2},
\end{equation}
and expressed in dB, as well as the perceptually motivated objective measure PEMO-Q \cite{Huber:2006a}.
% and PEAQ \cite{Thiede:2000a,Kabal2002:PEAQ}.
%, since no other standard for perceptual similarity is established for the evaluation of audio inpainting.
All the metrics measure the similarity
\edt{between} %of
the ground-truth signal $$ and its estimate $$; SNR represents the sample-wise similarity
\edt{between} %of
the waveforms, whereas PEMO-Q % and PEAQ
estimates the perceived difference (objective difference grade, ODG) on a scale ranging from $-4$ (very annoying) to $0$ (imperceptible).
Note that since we
\edtlang{assumed} %consider
a noise-less scenario and all the methods considered fit the reliable samples perfectly, we measure the SNR only on the inpainted segments.

The relative solution change, used as a measure of convergence, is computed as $^(i+1)-^(i)/^(i)$, where $^(i)$ is the estimate at the $i$-th iteration folded together from $= [_1,…,_N] = [_1,…,_N]$, where $$ is the output of the proposed algorithms.
Similarly, we
\edtlang{computed} %compute
the relative objective change, where the objective is the negative log-likelihood of Eq.\,\eqref{eq:mle.1} for \EMTF{} and \EMT{} and of Eq.\,\eqref{eq:mle.2} for \AM, evaluated using the current iterates of $$, $$̋ and the corresponding solution.

§.§ Comparison of the proposed estimators

the experiments with the comparison of the behavior of the three estimators, ,  and .
In <ref>, we
a preliminary comparison including the very demanding calculation of the objective function.
Then, we
give a validation of
the results in <ref> for the gap-filling task.

§.§.§ Introductory comparison with random missing samples

First, we illustrate the capabilities of the three estimators ,  and  on a demonstrative example of an inpainting problem.
The original signal is an excerpt of the first 12 seconds of the song Mamavatu by Susheela Raman, containing acoustic guitar and drums, sampled at 16 kHz.
60 % of the signal samples (chosen randomly) and
inpainting using ,  and .
The temporal frames
extracted using sine window of length $\Wdim = 1024$ samples (64 ms), the hop length
was set to
512 samples and we
$K = 20$ components of the NMF.
To distinguish between the individual algorithms,
we used DFT as the transform $\adjoint{\syn}$ with the number of frequency channels $\Fdim$

between $\Wdim$ and $2\Wdim$, which corresponds to examples <ref> and <ref> of subsection <ref>, respectively.
The inner loop for NMF was run for 10 iterations.

As already mentioned in Remark <ref>, choosing the DFT and its inverse transform as $\syn$ and $\adjoint{\syn}$, respectively, causes violation of the independence assumption <ref>.
A natural way of overcoming this issue would be to work with half of the frequency spectrum and only compute the complex-conjugate half when synthesizing the resulting signal.
However, this approach is not practical since all our algorithms require to compute the inverse $\inv{\left(\Mr_n\syn \D_n \adjoint{\syn}\transp{\Mr}_n\right)}$, which is quite involved.
Therefore, we used the whole frequency spectrum throughout the experiments.
To ensure that the resulting signal are real-valued, we initialized $\W$ such that its columns exhibit the desired symmetry.
It is straightforward to show that the symmetry is preserved throughout all the operations during the updates of our algorithms (element-wise addition and multiplication of such symmetrical matrices, right multiplication with a real-valued matrix, etc. - see  <cit.> for more details).

The comparison is visualized in Fig. <ref> by means of several quantities: the negative log-likelihood (i.e., the objective function of the estimation problem (<ref>) or (<ref>)), SNR, relative objective change and relative solution change,
as introduced in Subsection <ref>.
Our key finding is the difference in the observed convergence speed with respect to iteration count.
It is visible in all the plots of Fig. <ref> that
its solution
faster than , which is even more pronounced with respect to CPU time according to remark <ref> (rem:computational.complexity).
However, the quality of the solution was slightly worse and it
decreased after further iterations.

A similar behavior
observed for the case of redundant transform $\syn$.
This redundancy causes that:

* the convergence of  for $\Fdim = 2\Wdim$
slower than with $\Fdim = \Wdim$ while reaching similar reconstruction quality,

* the convergence of  for $\Fdim = 2\Wdim$
faster than both cases of , but as in the case of , it
worse reconstruction quality.

Based on these first observations, a natural question arises: Can we combine the convergence properties of  with the performance of ?
To answer it, we consider a combined algorithm , which consists in initializing  with $5$ iterations of .
As observed in Fig. <ref>, the initialization with
the algorithmic behavior
in the first few iterations and at the same time, the resulting restoration quality
the same as with .
However, the number of iterations needed to reach the peak performance
not reduced by the switching strategy.

Comparison of the performance and convergence properties of ,  and , including the switching variant .
The legend in the middle plot is common for the whole figure.
Note that the first column thus shows three different quantities, since the objective depends on the choice of $\Fdim$ and also on the algorithm.
Especially, the formula for log likelihood switches after the initializing iterations of , which is however disregard on purpose in the plot.

§.§.§ Missing gaps

To validate the
results concerning  and , we
a larger experiment with a set of signals and for the more demanding problem of inpainting short to middle-length gaps (instead of random subsampling).

the set of 10 musical recordings from the EBU SQAM dataset <cit.>, sampled at 44.1 kHz and shortened to 7 seconds, as used commonly in recent related publications <cit.>.

In each of the 10 signals from the EBU SQAM dataset, 10 gaps of given length
artificially introduced, with the gap length ranging from 20 to 80 ms.
not track the objective function in this case, because its computation
was computationally too demanding. We also
on the practical case of $\Fdim = \Wdim$ and invertible transform $\syn$ representing the inverse DFT in each temporal frame, thus
The frame length
was set to
$\Wdim = 4096$ samples (approx. 92 ms), and the temporal frames
extracted using sine window with 50% overlap.
We used $K=20$ components of the low-rank model according to assumption <ref>,
and the algorithms <ref> and <ref> were run with maximum number of iterations set to 100,
while the inner cycle for the NMF used 10 iterations.

As in the previous experiment, we
that the performance difference between  and
not significant, but
its peak faster, as shown in the plots for SNR in Fig. <ref>.
The relative solution change supports this observation, as it has been demonstrated above (see Fig. <ref>) that this measure mostly corresponds to the convergence of the algorithm with respect to its objective value.
A new observation is that this phenomenon depends on the gap length – the longer the gap, the slower the convergence of  is, compared to .

Comparison of the performance of , , and the switching variant  (switching after 5 iterations).
The left column shows the evolution of the SNR over iterations, the right column shows the relative solution change.
Both the metrics are averaged over the dataset and plotted together with 95% confidence interval represented by the light colored areas (note that for the relative solution change plot, this confidence interval is computed from the decimal logarithm of the data).

To provide an example supporting the observation about convergence properties of  and , we present the evolution in the temporal domain in Fig. <ref>.
We have taken a single signal from the gap-filling experiment and studied how the temporal solution changes during iterations inside a particular gap.
The observation is consistent with
that obtained from Fig. <ref>.
We see that  already identified the basic shape of the original signal in iteration 2 and then only improved the estimate in a few iterations.
On the other hand, the first estimate of  in iteration 2 seems to be almost the opposite, compared to the original signal.
As a result,  first returned to the solution being around zero, before finally reaching a quality estimate in later iterations.

The interpretation here is that the inclusion of the covariance matrices $\bSigmah_n$ in  helps to fit the model closer to the original signal, but only when the current estimate of the model (matrices $\W,\H$) is already relevant.
This corresponds to the observation that the initial estimates of the signal via \EMTF{} were observed to be inferior compared to \AM{}.
It is also consistent with Fig.\,\ref{fig:em_am}, which shows that the decrease observed in the performance of \AM{} could be prevented by including the covariance matrices into the estimation (i.e.\ by switching to \EMTF{}) if the estimate is already sufficient to some extent.}

\begin{figure}[ht]
\centering
\scalebox{0.6}{\input{figures/single_gap.tex}}
\caption{\edt{Example of the difference between \EMTF{} and \AM{} during iterations, namely for the guitar recording from our testing set and gap length 30\,ms.
The waveforms are compared to the undegraded original waveform (in black on the plots).
The legend indicates SNR values \emph{for the whole signal} and for each method at the corresponding iteration number.
The iterations are chosen such that the peak solution (in terms of SNR) is displayed both for \EMTF{} and \AM{} (iterations 18 and 4, respectively).
Also note that the solution in the first plot is the same for \EMTF{} and \AM{}, since both algorithms use the same initialization, thus there is no difference in their estimate after 1 iteration (\textit{cf.} line 5 of Alg.\,\ref{alg:emtf} and \ref{alg:am}).
\label{fig:single_gap}
\end{figure}

\subsection{Comparison with the state of the art}
\label{ssec:comparison.soa}

Finally, we
\edtlang{compared} %compare
our NMF-based methods to the following state-of-the-art techniques for the task of inpainting middle-length audio gaps,
\edt{based on recent publications \cite{MokryRajmic2020:Inpainting.revisited,TaubockRajbamshiBalasz2021:SPAINMOD,Rajbamshi2021:Audio.inapinting.ell.1.dictionary.learning}}:

\begin{itemize}
\item Janssen's AR approach \cite{javevr86} (denoted \textbf{Janssen}): This iterative algorithm builds upon a frame-wise AR nature of the clean signal.
At each iteration, it estimates the AR coefficients of the signal estimate (starting from the observed signal with the missing samples initialized as zeros), and then recompute the missing samples using the reliable samples and the current model parameters' estimates.
\item Modified variant of Janssen's algorithm (\textbf{Janssenmod}): This method is similar to the AR approach, but instead of processing the signal in overlapping segments, it treats each gap with its context individually.
\item Cross-faded extrapolation based on AR modeling \cite{Etter1996:Interpolation_AR} (\textbf{LR}): This is an efficient method for inpainting of compact gaps with sufficiently long reliable contexts.
It consists of estimating two AR models for the left and right contexts, extrapolating the contexts into the gap using the fitted model and cross-fading the two candidate solutions.
\item \textbf{SPAINlearned} \cite{TaubockRajbamshiBalasz2021:SPAINMOD}, a variant of the sparsity-based non-convex approach SPAIN \cite{MokryZaviskaRajmicVesely2019:SPAIN}, which treats individual gaps and their contexts (instead of overlapping frames as in SPAIN).
The performance is further enhanced by a dictionary learning step to deform the STFT based on the particular signal such that it allows for a sparser representation than using the STFT.
\item A convex alternative is the (weighted) $\ell_1$ minimization as a relaxation of the non-convex sparsity \cite{MokryRajmic2020:Inpainting.revisited} (\textbf{reweighted}).
\end{itemize}

For \redt{the} % all
methods which use STFT or segmentation (\textbf{Janssen}, \textbf{SPAINlearned}, \textbf{reweighted}, and the NMF-based methods),
\edtlang{we used} %the window is
a sine window of length 4096 samples (approx.\ 92\,ms) with 50\% overlap.
The NMF-based methods
\edtlang{were} %are
applied with $= $, thus \EMT{}
\edtlang{was} %is
equivalent to \EMTF{} and
\edtlang{was therefore} %is
omitted for brevity.
The AR-based methods
\edtlang{used} %use
a model of order 512.
The context of \textbf{Janssenmod} and \textbf{LR}
\edtlang{was} %is
set to 4096 samples, while \textbf{SPAINlearned}
\edtlang{used} %use
a longer context (8192 samples) for the sake of the dictionary learning. These values
\edtlang{were} %are
chosen based on the corresponding studies, where they have shown good performance.
For particular choices of all the parameters of the individual methods, please refer to the published source code.

The performance results averaged over the 10 test signals are shown in Fig.\,\ref{fig:comparison}.
We observe that in terms of SNR the proposed NMF-based methods
\edtlang{outperformed} %outperform
the state-of-the-art for short gaps (up to \redt{40}\,ms) while their performance
\edtlang{dropped} %drops
for longer gaps.
\edt{Note that Janssen % exhibits
\redt{exhibited} a~similar behavior, suggesting that the performance %is
\redt{was} undermined by estimating the missing samples via Wiener filtering (as mentioned in Sec.\,\ref{sec:formulation}) from a limited number of samples (which %is
\redt{was} more pronounced for long gaps with a large number of missing samples per frame).
Note that similar process %is conducted
\redt{appears} within the Janssen algorithm in the signal update step \cite[Sec.\,3.3]{Oudre2018:Janssen.implementation}.
On the other hand, sparsity-based methods tend to loose signal energy within long gaps \cite{MokryRajmic2020:Inpainting.revisited}, which may eventually lead to better SNR results.}

\edt{The perceptually-motivated comparison is realized via PEMO-Q ODG.
Here we}
% The perceptually-motivated comparison includes PEMO-Q and PEAQ ODG.
% However, based on the range of the PEAQ ODG values, we do not find these results very informative, and rather include them for the sake of completeness.
% On the other hand, we
\edtlang{observed} %observe
% from the PEMO-Q ODG values
that \EMTF{} and \AM{}
\edtlang{were} %are
among the top three methods for short to middle gaps \edt{(up to 40\,ms)}.
\edt{For longer gaps, the performance %decreases
\redt{decreased} similarly to the evaluation based on SNR.}

\begin{figure}[ht]
\centering
\scalebox{0.6}{\input{figures/figure3_revised.tex}}
\caption{Comparison with the state-of-the-art algorithms for inpainting short to middle-length gaps.
The legend is common for both plots.}
\label{fig:comparison}
\end{figure}

%\edt{To depict the results in more detail,
% OM: added this sentence
\redt{Since the results vary in general, depending on the particular signal,}
\redt{we present a short signal-specific analysis to illustrate which particular results contribute the most to the improvement over the state-of-the-art methods presented in Fig.\,\ref{fig:comparison}.
To this end,}
Fig.\,\ref{fig:comparison_per_signal} focuses on
\redt{three signals from the dataset on which the NMF-based methods perform particularly well in terms of SNR.}
%the signals for which the NMF-based methods reached a notable improvement over the state-of-the-art methods in terms of SNR.
\redt{While the results in terms of PEMO-Q ODG were mostly consistent with SNR,
we observed a notable difference in the case of the harp signal -- the performance of \AM{} in terms of SNR suddenly dropped for 50\,ms gaps, but the ODG grades stayed on par with \EMTF{}.}
%Interestingly, the observed inconsistency between SNR and PEMO-Q ODG favors the \AM{} algorithm in the case of the harp signal, where the significant drop in SNR is not projected into the ODG scores.

\begin{figure}[ht]
\centering
\scalebox{0.6}{% This file was created by matlab2tikz.
%The latest updates can be retrieved from
%  http://www.mathworks.com/matlabcentral/fileexchange/22022-matlab2tikz-matlab2tikz
%where you can also make suggestions and rate matlab2tikz.
% data:         larger_inpainting_comparison_merged.mat
% generated by: larger_inpainting_comparison.m,
%               larger_inpainting_comparison_NMF.m
% plotted by:   inpainting_comparison_plot.m
\definecolor{mycolor1}{rgb}{0.00000,0.44700,0.74100}%
\definecolor{mycolor2}{rgb}{0.85000,0.32500,0.09800}%
\definecolor{mycolor3}{rgb}{0.92900,0.69400,0.12500}%
\definecolor{mycolor4}{rgb}{0.49400,0.18400,0.55600}%
\definecolor{mycolor5}{rgb}{0.46600,0.67400,0.18800}%
\definecolor{mycolor6}{rgb}{0.30100,0.74500,0.93300}%
\definecolor{mycolor7}{rgb}{0.63500,0.07800,0.18400}%
\def\NMFlinethickness{0.85mm}
\begin{tikzpicture}

\begin{axis}[%
scale only axis,
xlabel style={font=\color{white!15!black}},
xlabel={gap length (ms)},
ylabel style={font=\color{white!15!black}},
ylabel={SNR (dB)},
axis background/.style={fill=white},
title style={font=\bfseries},
legend style={legend cell align=left, align=left, draw=white!15!black}
\addplot [color=mycolor1,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	19.1361172421887\\
30	15.6803392380297\\
40	13.6682344073348\\
50	12.5996994507056\\
60	10.5188936323033\\
70	6.67853560858092\\
80	6.63590756697866\\
% \addlegendentry{\EMTF{}}

\addplot [color=mycolor2,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	19.7414777004135\\
30	16.793687905635\\
40	15.5291559437398\\
50	7.54801705166557\\
60	7.5860562052103\\
70	6.36092914099383\\
80	6.45388004239599\\
% \addlegendentry{\AM{}}

\addplot [color=mycolor3]
table[row sep=crcr]{%
20	13.2227399564113\\
30	11.2909727851921\\
40	9.37328001251518\\
50	8.26036181354689\\
60	6.93454875084045\\
70	5.89628700165694\\
80	5.41941877891664\\
% \addlegendentry{Janssen}

\addplot [color=mycolor4]
table[row sep=crcr]{%
20	13.5629984548588\\
30	11.8166238261482\\
40	9.97652295796737\\
50	8.32662324592138\\
60	8.29226829557175\\
70	8.01436028958267\\
80	7.53761921457195\\
% \addlegendentry{Janssenmod}

\addplot [color=mycolor5]
table[row sep=crcr]{%
20	11.6425182343684\\
30	11.0699894921155\\
40	9.22609570388769\\
50	8.23670036160983\\
60	7.80844036437605\\
70	7.72089408783704\\
80	6.88862187897841\\
% \addlegendentry{LR}

\addplot [color=mycolor6]
table[row sep=crcr]{%
20	17.5510353347415\\
30	14.0333066322075\\
40	12.617247664234\\
50	11.5562705128235\\
60	11.3288911622004\\
70	8.16774278926693\\
80	7.73199790189495\\
% \addlegendentry{SPAINlearned}

\addplot [color=mycolor7]
table[row sep=crcr]{%
20	13.5975973169458\\
30	11.382445902619\\
40	9.77754489992811\\
50	9.64295288628557\\
60	8.36236057748718\\
70	7.28481780016079\\
80	5.18265398659583\\
% \addlegendentry{reweighted}

\end{axis}

\begin{axis}[%
scale only axis,
xlabel style={font=\color{white!15!black}},
xlabel={gap length (ms)},
ylabel style={font=\color{white!15!black}},
ylabel={SNR (dB)},
axis background/.style={fill=white},
title style={font=\bfseries},
legend style={legend cell align=left, align=left, draw=white!15!black}
\addplot [color=mycolor1,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	18.1485742675895\\
30	14.2704783263745\\
40	12.0644040356025\\
50	9.36867245074381\\
60	7.3404741150593\\
70	3.51008511476938\\
80	2.9314623795503\\
% \addlegendentry{\EMTF{}}

\addplot [color=mycolor2,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	17.0474898226938\\
30	14.7166200809083\\
40	12.9550605585822\\
50	10.3702063500069\\
60	7.82366459434755\\
70	4.2547347992531\\
80	3.37193053598567\\
% \addlegendentry{\AM{}}

\addplot [color=mycolor3]
table[row sep=crcr]{%
20	15.7481839832033\\
30	12.3474309574973\\
40	10.559396660095\\
50	7.84160924011624\\
60	3.60250875996043\\
70	3.03009553613289\\
80	2.58388021816696\\
% \addlegendentry{Janssen}

\addplot [color=mycolor4]
table[row sep=crcr]{%
20	15.6165904711584\\
30	13.0260311512341\\
40	11.2082721526508\\
50	8.25721526184968\\
60	5.91113503402587\\
70	4.99530345609224\\
80	3.61884434236907\\
% \addlegendentry{Janssenmod}

\addplot [color=mycolor5]
table[row sep=crcr]{%
20	9.97227026978335\\
30	7.45522147358752\\
40	5.58773152561338\\
50	4.51933996112026\\
60	3.96501759760622\\
70	3.2513357483852\\
80	2.79561710896413\\
% \addlegendentry{LR}

\addplot [color=mycolor6]
table[row sep=crcr]{%
20	12.0344902322419\\
30	11.2386238211\\
40	8.54511396635549\\
50	7.32823191814608\\
60	5.21292587405024\\
70	4.98218940599571\\
80	3.62479966584014\\
% \addlegendentry{SPAINlearned}

\addplot [color=mycolor7]
table[row sep=crcr]{%
20	12.847110962485\\
30	8.63398534710363\\
40	5.45988634805048\\
50	4.14849470988096\\
60	3.26534914607855\\
70	2.78252231784103\\
80	2.11360451778408\\
% \addlegendentry{reweighted}

\end{axis}

\begin{axis}[%
scale only axis,
xlabel style={font=\color{white!15!black}},
xlabel={gap length (ms)},
ylabel style={font=\color{white!15!black}},
ylabel={SNR (dB)},
axis background/.style={fill=white},
title style={font=\bfseries},
legend style={legend cell align=left, align=left, draw=white!15!black}
\addplot [color=mycolor1,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	8.28240897050337\\
30	7.00994850417444\\
40	5.54042492869449\\
50	3.78761525714733\\
60	2.55803367860929\\
70	2.4415559209961\\
80	2.14624028471493\\
\addlegendentry{\EMTF{}}

\addplot [color=mycolor2,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	8.54934870739845\\
30	5.32828901585732\\
40	4.06304154647852\\
50	3.89644432879172\\
60	2.41392970669487\\
70	2.38218020152002\\
80	1.77818985019129\\
\addlegendentry{\AM{}}

\addplot [color=mycolor3]
table[row sep=crcr]{%
20	3.24519126017131\\
30	2.88303826478977\\
40	2.30874891316369\\
50	1.51787624304036\\
60	1.10118076445445\\
70	0.787342730627662\\
80	0.784559205843556\\
\addlegendentry{Janssen}

\addplot [color=mycolor4]
table[row sep=crcr]{%
20	3.49230579371693\\
30	3.40798829420266\\
40	3.17116240092692\\
50	2.23402159890218\\
60	2.15788845193063\\
70	2.08258905633536\\
80	1.62848981948071\\
\addlegendentry{Janssenmod}

\addplot [color=mycolor5]
table[row sep=crcr]{%
20	3.0121766417137\\
30	2.79228397794626\\
40	2.35929867164166\\
50	1.86554890055172\\
60	1.54716444150247\\
70	1.47512446626154\\
80	1.18660699569479\\
\addlegendentry{LR}

\addplot [color=mycolor6]
table[row sep=crcr]{%
20	4.4164329578349\\
30	3.8144306043834\\
40	5.01973495244824\\
50	5.07142743812257\\
60	4.67546865171506\\
70	5.49701528681812\\
80	4.34052141370421\\
\addlegendentry{SPAINlearned}

\addplot [color=mycolor7]
table[row sep=crcr]{%
20	3.3426284409104\\
30	3.22023335319325\\
40	2.81587280898705\\
50	2.59388472369284\\
60	1.98153437033287\\
70	1.72351442878818\\
80	0.641471506798929\\
\addlegendentry{reweighted}

\end{axis}

\begin{axis}[%
scale only axis,
xlabel style={font=\color{white!15!black}},
xlabel={gap length (ms)},
ylabel style={font=\color{white!15!black}},
ylabel={PEMO-Q ODG},
axis background/.style={fill=white},
title style={font=\bfseries},
legend style={legend cell align=left, align=left, draw=white!15!black}
\addplot [color=mycolor1,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	-0.302031524108866\\
30	-0.519691483338935\\
40	-0.869670251917254\\
50	-1.5091650435296\\
60	-2.26502528797598\\
70	-3.12879104883682\\
80	-3.48956616558837\\
% \addlegendentry{\EMTF{}}

\addplot [color=mycolor2,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	-0.301210621030027\\
30	-0.59006032744594\\
40	-0.904371165946728\\
50	-1.43288710887187\\
60	-2.60138603358014\\
70	-2.87399223177433\\
80	-3.18018994811886\\
% \addlegendentry{\AM{}}

\addplot [color=mycolor3]
table[row sep=crcr]{%
20	-0.352137957825633\\
30	-0.707399607423694\\
40	-1.15525035829876\\
50	-1.68021410741906\\
60	-2.40198986928895\\
70	-3.05667490577138\\
80	-3.22492393027218\\
% \addlegendentry{Janssen}

\addplot [color=mycolor4]
table[row sep=crcr]{%
20	-0.702996390968165\\
30	-1.27561861337595\\
40	-1.6181924511861\\
50	-2.22527952941075\\
60	-2.31203388188366\\
70	-2.82137612845688\\
80	-3.02891772959478\\
% \addlegendentry{Janssenmod}

\addplot [color=mycolor5]
table[row sep=crcr]{%
20	-1.42044503641207\\
30	-2.15506249926501\\
40	-2.56709741758789\\
50	-2.98993946416385\\
60	-3.04754339108381\\
70	-3.14517258145656\\
80	-3.29465165709315\\
% \addlegendentry{LR}

\addplot [color=mycolor6]
table[row sep=crcr]{%
20	-0.542268451674895\\
30	-1.16047816559918\\
40	-1.35284554361613\\
50	-2.23501315200126\\
60	-2.58091845464103\\
70	-2.63699422097931\\
80	-2.88900287260838\\
% \addlegendentry{SPAINlearned}

\addplot [color=mycolor7]
table[row sep=crcr]{%
20	-0.726920506872897\\
30	-1.53781017733781\\
40	-2.22660946984008\\
50	-2.42834532899414\\
60	-2.79082798521607\\
70	-3.07118087934658\\
80	-3.40347560125479\\
% \addlegendentry{reweighted}

\end{axis}

\begin{axis}[%
scale only axis,
xlabel style={font=\color{white!15!black}},
xlabel={gap length (ms)},
ylabel style={font=\color{white!15!black}},
ylabel={PEMO-Q ODG},
axis background/.style={fill=white},
title style={font=\bfseries},
legend style={legend cell align=left, align=left, draw=white!15!black}
\addplot [color=mycolor1,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	-0.646022303045996\\
30	-1.42075980260925\\
40	-2.35647740234348\\
50	-2.98141883118414\\
60	-3.22441648511898\\
70	-3.47247640252925\\
80	-3.53374111213868\\
% \addlegendentry{\EMTF{}}

\addplot [color=mycolor2,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	-0.776402066276537\\
30	-1.32639461796883\\
40	-2.00187122287439\\
50	-2.59717241120742\\
60	-3.27756620474786\\
70	-3.43374320346454\\
80	-3.48786146992494\\
% \addlegendentry{\AM{}}

\addplot [color=mycolor3]
table[row sep=crcr]{%
20	-0.23396510444001\\
30	-0.486636262098871\\
40	-0.996856487838537\\
50	-1.80601208429609\\
60	-2.66114466939698\\
70	-3.27133242682999\\
80	-3.37364209376629\\
% \addlegendentry{Janssen}

\addplot [color=mycolor4]
table[row sep=crcr]{%
20	-0.442189403454966\\
30	-0.843968692379102\\
40	-1.57733894662096\\
50	-2.08593995901462\\
60	-2.64998312490493\\
70	-2.92910499100512\\
80	-3.0232406161433\\
% \addlegendentry{Janssenmod}

\addplot [color=mycolor5]
table[row sep=crcr]{%
20	-0.999595108308892\\
30	-1.87483074241384\\
40	-2.71243712200472\\
50	-3.05110565574591\\
60	-3.19250977447412\\
70	-3.34143765636918\\
80	-3.43130762206045\\
% \addlegendentry{LR}

\addplot [color=mycolor6]
table[row sep=crcr]{%
20	-0.374755737460941\\
30	-0.991105887807402\\
40	-1.87826274209843\\
50	-2.53553354916069\\
60	-2.83455589899156\\
70	-3.09997226957493\\
80	-3.30606336629998\\
% \addlegendentry{SPAINlearned}

\addplot [color=mycolor7]
table[row sep=crcr]{%
20	-0.926218718599799\\
30	-1.86825554272934\\
40	-2.97027871023008\\
50	-3.27296822471774\\
60	-3.38787710709368\\
70	-3.45661616987944\\
80	-3.57997348392005\\
% \addlegendentry{reweighted}

\end{axis}

\begin{axis}[%
scale only axis,
xlabel style={font=\color{white!15!black}},
xlabel={gap length (ms)},
ylabel style={font=\color{white!15!black}},
ylabel={PEMO-Q ODG},
axis background/.style={fill=white},
title style={font=\bfseries},
legend style={legend cell align=left, align=left, draw=white!15!black}
\addplot [color=mycolor1,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	-0.480367493982087\\
30	-1.06329629922404\\
40	-1.86956723450374\\
50	-2.66778305327009\\
60	-3.08754991960833\\
70	-3.26673467531541\\
80	-3.46522605155436\\
% \addlegendentry{\EMTF{}}

\addplot [color=mycolor2,line width=\NMFlinethickness]
table[row sep=crcr]{%
20	-0.496467021285349\\
30	-1.08377450746851\\
40	-1.73578746028709\\
50	-2.7522913607735\\
60	-3.14642608343068\\
70	-3.26474808057221\\
80	-3.33808918928689\\
% \addlegendentry{\AM{}}

\addplot [color=mycolor3]
table[row sep=crcr]{%
20	-0.888747699947571\\
30	-1.74929839339753\\
40	-2.35918821489418\\
50	-2.81131161013796\\
60	-3.04051214115639\\
70	-3.31094154004197\\
80	-3.41177863220069\\
% \addlegendentry{Janssen}

\addplot [color=mycolor4]
table[row sep=crcr]{%
20	-1.32137678241264\\
30	-1.90002312470463\\
40	-2.47653538860882\\
50	-2.96550089446814\\
60	-3.02310316109408\\
70	-3.22869755684609\\
80	-3.34478089054316\\
% \addlegendentry{Janssenmod}

\addplot [color=mycolor5]
table[row sep=crcr]{%
20	-1.96931487951185\\
30	-2.74790104414906\\
40	-3.03655547871729\\
50	-3.15137313435243\\
60	-3.28379199519981\\
70	-3.42192245319703\\
80	-3.49378174889751\\
% \addlegendentry{LR}

\addplot [color=mycolor6]
table[row sep=crcr]{%
20	-0.35816031164331\\
30	-0.662436132520574\\
40	-1.00836440941263\\
50	-1.4084286955112\\
60	-2.18802991482414\\
70	-2.51075639359778\\
80	-2.78310659965622\\
% \addlegendentry{SPAINlearned}

\addplot [color=mycolor7]
table[row sep=crcr]{%
20	-1.22925557429904\\
30	-1.92233302905355\\
40	-2.33282379399973\\
50	-2.84539378247752\\
60	-3.24796657430684\\
70	-3.43370457863927\\
80	-3.56008357112456\\
% \addlegendentry{reweighted}

\end{axis}

\end{tikzpicture}%}
\caption{\edt{Comparison with the state-of-the-art algorithms for inpainting short to middle-length gaps on some selected signals.
The legend is common for all the plots.}}
\label{fig:comparison_per_signal}
\end{figure}

\edt{To summarize the results,
\AM{} reaches a similar performance to \EMTF{} but with fewer iterations (and reduced CPU time given the lower complexity of \AM, see remark \ref{rem:computational.complexity}),
and both methods outperform the state-of-the-art in terms of SNR for small gaps (and rank second in terms of ODG), but they %are less competitive
\redt{underperformed the other algorithms} for medium to large gaps.}

\section{Conclusion}
\label{sec:conclusion}

In this paper, we derived new estimators for NMF-based audio inpainting.
We formulated inpainting as an optimization problem where the goal is to estimate the signal's power spectrum, which is structured using NMF.
To that end, we have derived three algorithms, among which two are new, which encompasses and extends previous related works~\cite{BilenOzerovPerez2015:declipping.via.NMF,Bilen2018:NTF_audio_inverse_problems}.
Even though the proposed estimators build upon the same low-rank assumption about the signal's TF spectrum, we have shown that there are both theoretical and practical differences between them.
In particular,
they all exhibit a different behavior and lead to different solutions to the audio inpainting problem.
Importantly, the novel approaches (\EMT{} and \AM{}) improve the convergence rate compared to \EMTF{}, while reaching similar reconstruction quality.
% They have also been demonstrated competitive with state-of-the-art audio inpainting methods.
\edt{%The
\redt{Both} NMF-based methods \redt{considered in the comparison} are in general efficient for short to middle gaps (up to 40\,ms), where they %improve the state-of-the-art.
\redt{provided better results than the state-of-the-art algorithms}.}

Throughout the derivation, we assumed independence of the temporal frames.
A natural extension of the model would be to employ temporal Markov NMF models, as presented e.g., in \cite{Smaragdis2014:Static.dnymic.source.separation.NMF}.
Future research will also study the possibility to leverage psychoacoustics in audio inpainting, e.g., by using perceptually-motivated TF representations \cite{Lieb2018:Audio.Inpainting, Necciari2013:ERBlet}.

\bibliography{literatura.bib}

\end{document}


# Intrinsic Nonlinear Spin Magnetoelectricity in Centrosymmetric Magnets

Cong Xiao  Department of Physics, The University of Hong Kong, Hong Kong,
China  HKU-UCAS Joint Institute of Theoretical and Computational Physics at
Hong Kong, China  Huiying Liu<EMAIL_ADDRESS>Research Laboratory for
Quantum Materials, Singapore University of Technology and Design, Singapore
487372, Singapore  Weikang Wu<EMAIL_ADDRESS>Research Laboratory for
Quantum Materials, Singapore University of Technology and Design, Singapore
487372, Singapore  Division of Physics and Applied Physics, School of Physical
and Mathematical Sciences, Nanyang Technological University, Singapore 637371,
Singapore  Hui Wang  Research Laboratory for Quantum Materials, Singapore
University of Technology and Design, Singapore 487372, Singapore  Qian Niu
School of Physics, University of Science and Technology of China, Hefei, Anhui
230026, China  Shengyuan A. Yang  Research Laboratory for Quantum Materials,
Singapore University of Technology and Design, Singapore 487372, Singapore

###### Abstract

We propose an intrinsic nonlinear spin magnetoelectric effect in magnetic
materials, offering the potential of all-electric control of spin degree of
freedom in centrosymmetric magnets, which reside outside of the current
paradigm based on linear spin response. We reveal the band geometric origin of
this effect in the momentum and magnetization space Berry connection
polarizabilities, and clarify its symmetry characters. As an intrinsic effect,
it is determined solely by the material’s band structure and represents a
material characteristic. Combining our theory with first-principles
calculations, we predict sizable nonlinear spin magnetoelectricity in single-
layer MnBi2Te4, which can be detected in experiment. Our theory paves the way
for exploring rich nonlinear spintronic effects and novel device concepts
based on them.

Electric control of spin degree of freedom is a central topic of spintronics.
Much theoretical and experimental effort Pikus1978 ; Aronov1989 ; Edelstein ;
Kato2004 ; Geller2009 ; Vyborny2011 ; Manchon2019  has been devoted to
generating spin polarization $\delta\bm{s}$ by an applied electric field,
characterized by a _linear_ spin magnetoelectric response tensor $\alpha_{ij}$
with $\delta s_{i}=\alpha_{ij}E_{j}$, where $i$, $j$ are Cartesian indices and
the Einstein summation convention is adopted. In a magnet, the generated spin
polarization may further induce torques on the magnetization and cause
magnetic reversal Manchon2019 . Importantly, since spin is even under space
inversion but electric field is odd, the linear effect is constrained to
systems with broken inversion symmetry Pikus1978 ; Aronov1989 ; Edelstein ;
Kato2004 ; Geller2009 ; Vyborny2011 ; Culcer2007 ; Manchon2019 ; Garate2009 ;
Franz2010 ; Kurebayashi2014 ; Zelezny2014 ; Freimuth2014 . For the class of
centrosymmetric magnets, the linear spin response is strictly forbidden in the
bulk and may only occur at interfaces when forming hetero-junctions with other
materials Miron2010 ; Miron2011 ; Liu2012 .

In this Letter, we unveil that sizable intrinsic nonlinear spin
magnetoelectric effect can exist in centrosymmetric magnets, thus
substantially extending the playing field of spin-charge conversion and
magnetoelectricity. Since the linear response is forbidden, the leading
contribution is of the second order:

$\delta s_{i}=\alpha_{ij\ell}E_{j}E_{\ell},$ (1)

characterized by a nonlinear response tensor $\alpha_{ij\ell}$. We show that
$\alpha_{ij\ell}$ contains an _intrinsic_ part determined solely by the band
structure. For insulators, $\alpha_{ij\ell}$ is closely connected to the
electric polarization of Bloch electrons and can be expressed by an important
band geometric quantity, the momentum space Berry connection polarizability
(BCP). For metals, there is an extra Fermi surface contribution, which
involves a magnetization space BCP. We clarify the symmetry properties of the
nonlinear intrinsic spin magnetoelectric effect. The intrinsic response has
the advantage of allowing a quantitative evaluation. By combining our theory
with first-principles calculations, we study the effect in single-layer
MnBi2Te4 and find sizable result that can be detected in experiment. Our work
develops the first theory for nonlinear electric spin generation, uncovers the
important roles of BCPs in spintronic effects, and opens the door to new
nonlinear spintronic device concepts. The approach here also offers a general
recipe for investigating other intrinsic nonlinear response properties of
Bloch electrons.

Thermodynamic argument for insulators. We first present a thermodynamic
argument, which applies to insulating cases and captures both linear and
nonlinear spin magnetoelectric effects. It also helps to expose the role of
electric polarization in the effect.

To evaluate the spin response, a conventional way is to introduce a fictitious
(homogeneous) Zeeman-like field $\bm{m}$ that couples to spin in the form of
$-\hat{\bm{{s}}}\cdot\bm{m}$, with $\hat{\bm{s}}$ the spin operator. This
auxiliary field is to be distinguished from the genuine magnetization of the
system and is set to zero at the end of the calculation Dong2020 . Under the
$\bm{m}$ field and the electric field, the electronic enthalpy of the magnetic
insulator system follows the relation $d\mathcal{H}=\mathcal{-}\bm{s}\cdot
d\bm{m}-\bm{P}\cdot d\bm{E}$, where $\bm{s}$ and $\bm{P}$ denote the spin
magnetization and the electric polarization of electrons, respectively.
According to the Maxwell relation, we have

$\frac{\partial s_{i}}{\partial E_{j}}=\frac{\partial P_{j}}{\partial m_{i}}.$
(2)

Therefore, the spin magnetoelectric response can be extracted from studying
the electric polarization.

To include the second order spin generation, we need the polarization expanded
to the first order in the electric field:
$\bm{P}=\bm{P}^{\left(0\right)}+\bm{P}^{\left(1\right)}$. The zero-field part
$\bm{P}^{(0)}=-\int[d\bm{k}]\bm{\mathcal{A}}(\bm{k})$ is well known KS1993 ;
Resta1994 ; Xiao2010 ; note  (we set $e=\hbar=1$),
$\bm{\mathcal{A}}\left(\bm{k}\right)=\langle
u_{n}(\bm{k})|i\partial_{\bm{k}}|u_{n}(\bm{k})\rangle$ is the intraband Berry
connection for a band eigenstate $|u_{n}(\bm{k})\rangle$, $[d\bm{k}]$ is a
short-hand notation for $\sum_{n}d\bm{k}/(2\pi)^{d}$ with $d$ being the
dimension of the system, and the summation is over all occupied bands. Here
and hereafter, for simple notations, we drop the band index $n$ wherever
appropriate. $\bm{P}^{(0)}$ gives the linear spin magnetoelectric response
with $\delta\bm{s}=\partial_{\bm{m}}P^{(0)}_{j}E_{j}$, reproducing the result
in previous works Yuriy2017 ; Xiao2021 .

On the other hand, the second order spin response is contained in
$\bm{P}^{(1)}$, the polarization linear in ${\bm{E}}$ field Sipe1995 ;
Nunes2001 ; Vanderbilt2002 . It can be expressed as
$\bm{P}^{(1)}=-\int[d\bm{k}]\bm{a}(\bm{k})$ Gao2014  in terms of the field-
induced Berry connection $a_{i}=G_{ij}E_{j}$, where

$G_{ij}=2\text{Re}\sum_{n^{\prime}\neq
n}\frac{(v_{i})_{nn^{\prime}}(v_{j})_{n^{\prime}n}}{(\varepsilon_{n}-\varepsilon_{n^{\prime}})^{3}}$
(3)

is known as the BCP for the state $|u_{n}(\bm{k})\rangle$ Gao2014 ,
$(v_{i})_{nn^{\prime}}$ is the interband velocity matrix element, and
$\varepsilon_{n}$ is the unperturbed energy for $|u_{n}(\bm{k})\rangle$.
Combining this result with Eq. (2), we immediately find the second order spin
polarization

$\delta\bm{s}=\partial_{\bm{m}}\Big{[}\frac{1}{2}\bm{E}\cdot\bm{P}^{(1)}\Big{]}$
(4)

with the nonlinear response tensor

$\alpha_{ij\ell}=-\frac{1}{2}\partial_{m_{i}}\int[d\bm{k}]G_{j\ell}(\bm{k})$
(5)

expressed nicely in terms of the momentum space BCP of occupied states. Recent
studies have highlighted the role of this BCP in nonlinear transport phenomena
Gao2014 ; Gao2017 ; Gao2018 ; Xiao2021adiabatic ; Lai2021 ; Wang2021 ; Liu2021
; Liu2022 , whereas our result here unveils its significance in
magnetoelectric effects.

The thermodynamic argument reveals the important role of electric polarization
in the spin magnetoelectricity. Particularly, the term in the bracket of Eq.
(4) is just the material-dependent part of the electric energy density in a
dielectric that is of $E^{2}$ order. For centrosymmetric systems, the zero-
field polarization $\bm{P}^{(0)}$ vanishes, so the argument confirms that the
linear spin response must also vanish. Meanwhile, $\bm{P}^{(1)}$ can be
nonzero regardless of the inversion symmetry, so that the nonlinear spin
response would become dominant in centrosymmetric materials.

It is also important to note that the argument above applies only to
insulators. For metals, the electric polarization ceases to be well defined,
hence we need a more general approach to the problem. As we shall see, in a
magnetic metal, there will be additional nonlinear contributions from the
Fermi surface.

Intrinsic nonlinear spin magnetoelectricity. To establish a general result
which is applicable also to metallic cases, we develop a semiclassical theory
for Bloch electrons in the nonlinear response regime. In Refs. Gao2014 ;
Gao2015 , Gao _et al._ extended the semiclassical theory to second order
accuracy. Nonetheless, the formulation there is focusing on the charge degree
of freedom, but does not explicitly handle spin. Here, we add this missing
piece. As is detailed in the Supplemental Material supp , within the extended
semiclassical framework, we derive the following spin expectation value
corrected to second order of $E$ field for an electron wave packet centered at
$|u_{n}(\bm{k})\rangle$:

$\bm{s}_{n}(\bm{k})=-\partial_{\bm{m}}\tilde{\varepsilon}_{n}+\bm{{\Omega}}_{\bm{m}\bm{k}}\cdot\bm{E}.$
(6)

Here, $\tilde{\varepsilon}_{n}=\varepsilon_{n}-(1/2)G_{ij}E_{i}E_{j}$ is the
field-corrected band energy, and
$(\bm{{\Omega}}_{\bm{m}\bm{k}})_{ij}=\partial_{m_{i}}(\mathcal{A}_{j}+a_{j})-\partial_{k_{j}}(\mathfrak{A}_{i}+\mathfrak{a}_{i})$
is the field-corrected Berry curvature in the hybrid $k$-$m$ space. The
definitions of the Berry connections $\mathfrak{\bm{A}}$ and
$\mathfrak{\bm{a}}$ are analogous to their counterparts $\mathcal{\bm{A}}$ and
${\bm{a}}$ in $k$ space. Specifically, $\mathfrak{A}_{i}=\langle
u_{n}(\bm{k})|i\partial_{m_{i}}|u_{n}(\bm{k})\rangle$, and
$\mathfrak{a}_{i}=\mathfrak{G}_{ij}E_{j}$ can be expressed using a $m$-space
BCP

$\mathfrak{G}_{ij}=-2\text{Re}\sum_{n^{\prime}\neq
n}\frac{(s_{i})_{nn^{\prime}}(v_{j})_{n^{\prime}n}}{(\varepsilon_{n}-\varepsilon_{n^{\prime}})^{3}},$
(7)

where the numerator involves the interband matrix elements of spin and
velocity operators. Like $G_{ij}$, $\mathfrak{G}_{ij}$ is gauge invariant, so
it is also an intrinsic band geometric property.

Table 1: Constraints on the intrinsic nonlinear spin magnetoelectric response tensor elements from magnetic point group symmetries. “$\checkmark$” (“$\times$”) means that the element is symmetry allowed (forbidden). Here, we choose to symmetrize the second and the third tensor indices, by defining $\alpha_{i(xy)}\equiv\frac{1}{2}(\alpha_{ixy}+\alpha_{iyx})$, and we omit the superscript ‘int’ in the table. Symmetry operations $\mathcal{T}$, $\mathcal{PT}$, $C_{3}\mathcal{T}$, and $S_{6}\mathcal{T}$ forbid all the elements here, hence are not listed. |  $\mathcal{P}$ |  $C_{2}^{z}$ |  $C_{3}^{z}$, |  $C_{4,6}^{z}$, |  $C_{2,4,6}^{x}$, |  $C_{3}^{x}$, |  $\sigma_{z}$ |  $\sigma_{x}$ |  $C_{2}^{z}\mathcal{T}$ |  $C_{4}^{z}\mathcal{T}$, |  $C_{6}^{z}\mathcal{T}$ |  $C_{2}^{x}\mathcal{T}$ |  $C_{4}^{x}\mathcal{T}$, |  $C_{6}^{x}\mathcal{T}$ |  $\sigma_{z}\mathcal{T}$ |  $\sigma_{x}\mathcal{T}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  $S_{6}^{z}$ |  $S_{4}^{z}$ |  $S_{4}^{x}$ |  $S_{6}^{x}$ |  |  |  |  $S_{4}^{z}\mathcal{T}$ |  |  |  $S_{4}^{x}\mathcal{T}$ |  |  |
$\alpha_{xxx}$ |  $\checkmark$ |  $\times$ |  $-\alpha_{xyy}$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $-\alpha_{xyy}$ |  $\times$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\times$
$\alpha_{x(xy)}$ |  $\checkmark$ |  $\times$ |  $\alpha_{yxx}$ |  $\times$ |  $\times$ |  $\times$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\alpha_{yxx}$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\checkmark$
$\alpha_{xyy}$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$
$\alpha_{yxx}$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\times$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\checkmark$
$\alpha_{y(xy)}$ |  $\checkmark$ |  $\times$ |  $\alpha_{xyy}$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\alpha_{xyy}$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$
$\alpha_{yyy}$ |  $\checkmark$ |  $\times$ |  $-\alpha_{yxx}$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\times$ |  $-\alpha_{yxx}$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
$\alpha_{zxx}$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\times$ |  $\checkmark$
$\alpha_{z(xy)}$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\times$
$\alpha_{zyy}$ |  $\checkmark$ |  $\checkmark$ |  $\alpha_{zxx}$ |  $\alpha_{zxx}$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\times$ |  $-\alpha_{zxx}$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\times$ |  $\checkmark$

With the spin polarization for each state, the total spin polarization in the
system can be obtained as

$\bm{s}=\int[d\bm{k}]\bm{s}_{n}(\bm{k})f_{n}(\bm{k}),$ (8)

where $f$ is the electron distribution function. The second-order spin
response is obtained by inserting the expression in (6) and retaining terms
that are of $E^{2}$ order. Here, we are particularly interested in the
intrinsic contribution that involves only the equilibrium Fermi distribution
function $f_{0}({\varepsilon}_{n})$ of the unperturbed band structure. The
intrinsic nonlinear response tensor is obtained as

$\begin{split}\alpha_{ij\ell}^{\text{int}}=-\frac{1}{2}\partial_{m_{i}}\int[d\bm{k}]G_{j\ell}f_{0}-\int[d\bm{k}](s_{i}G_{j\ell}+v_{j}\mathfrak{G}_{i\ell})f_{0}^{\prime},\end{split}$
(9)

where $s_{i}$ ($v_{j}$) are the intraband spin (velocity) matrix elements for
the state $|u_{n}(\bm{k})\rangle$.

Equation (9) is the key result of this work. First, it applies to both
insulators and metals. Compared to (5), Eq. (9) contains an additional Fermi
surface term (the second term). For the case of an insulator, the Fermi
surface term vanishes, and the result recovers (5), confirming the consistency
between the two approaches. Second, $\alpha_{ij\ell}^{\text{int}}$ is
suppressed by time reversal symmetry, and is nonzero only for magnetic
systems, as can be verified directly from Eq. (9). Third, as an intrinsic
contribution, Eq. (9) is a genuine material property, determined solely by the
material’s band structure. We mention that the $m$-derivative in the first
term of (9) can be done straightforwardly to obtain an expression involving
only the spin and velocity matrix elements [see Eq. (S6) in supp ] of the band
structure. Hence, the response can be readily evaluated in first-principles
calculations.

Symmetry property. As we have discussed, the intrinsic nonlinear response is
not suppressed by the inversion symmetry. Meanwhile, other crystalline
symmetries also put constraints on the form of $\alpha_{ij\ell}^{\text{int}}$,
which we analyze here.

Since spin is a time reversal ($\mathcal{T}$) odd pseudovector and the
electric field is a $\mathcal{T}$ even vector, $\alpha_{ij\ell}^{\text{int}}$
transforms as a third-rank $\mathcal{T}$ odd pseudotensor, which respects

$\alpha^{\text{int}}_{i^{\prime}j^{\prime}\ell^{\prime}}=\eta_{T}\text{det}(O)O_{i^{\prime}i}O_{j^{\prime}j}O_{\ell^{\prime}\ell}\alpha^{\text{int}}_{ij\ell}.$
(10)

Here $O$ is a point group operation, and the factor $\eta_{T}=\pm$ is
connected with the character of $\alpha_{ij\ell}^{\text{int}}$ being
$\mathcal{T}$ odd: $\eta_{T}=-1$ for primed operations, i.e., the magnetic
symmetry operations of the form $R\mathcal{T}$ with $R$ a spatial operation;
and $\eta_{T}=+1$ for nonprimed operations.

Assuming the applied electric field is in the $x$-$y$ plane, the constraints
from different magnetic point group symmetries are listed in Table 1. This
offers useful guidance for analyzing the nonlinear spin response for a
particular material. For example, consider a ferromagnet which preserves the
inversion and a horizontal mirror $\sigma_{z}$, with its magnetization along
the $z$ direction. Then, Table 1 tells us that for an applied in-plane
electric field, the generated nonlinear spin polarization must be out-of-
plane, along the magnetization direction. For the purpose to induce magnetic
reversal, one may want the induced spin polarization to have a component
normal to the magnetization, such that it can generate a torque, which means
that the desired material should not have a horizontal mirror plane.

Figure 1: (a) Top and (b) side views of the lattice structure of single-layer
MnBi2Te4. (c) shows the Brillouin zone. (d) Calculated band structure (spin-
orbit coupling is included). (e) Calculated nonlinear response tensor element
$\alpha_{yxx}^{\text{int}}$ versus the chemical potential $\mu$. (f)
Distribution of $\alpha_{yxx}^{\text{int}}(\bm{k})$, i.e., the integrand of
Eq. (9), in the momentum space for chemical potential at $-0.1$ eV [marked by
the red arrow in (e)]. The unit is $\mu_{B}/\text{V}^{2}$ per unit cell. In
the calculation, the temperature is set to 8 K.

A material example. We demonstrate the implementation of our theory in first-
principles calculations to study a concrete material. Guided by the symmetry
constraints in Table 1, we consider the effect in single-layer MnBi2Te4.
MnBi2Te4 in its bulk and two dimensional few-layer forms has attracted
considerable research interest recently, because it provides a platform for
realizing various types of topological states Gong2019 ; Chulkov2019 ;
Deng2020 ; Zhang2019 ; Xu2019 ; Otrokov2019 ; Zeugner2019 ; Yan2019 ; Lee2019
; Cui2019 ; Liu2020 . Our focus here is on its single layer, which has been
successfully fabricated in experiment, either by exfoliation from the bulk or
by molecular beam epitaxy growth Deng2020 ; Liu2020 ; Gong2019 . Its crystal
structure is shown in Fig. 1(a) and 1(b), characterized by a hexagonal lattice
with the space group $P\bar{3}m1$ (No. 164) and the point group $D_{3d}$. It
consists of seven atomic layers, stacked in the sequence of Te-Bi-Te-Mn-Te-Bi-
Te. Previous works have established that the ground state of single-layer
MnBi2Te4 is a topologically trivial ferromagnetic semiconductor with out-of-
plane magnetization, and the Curie temperature is about 12 K Chulkov2019 . The
ground-state magnetic configuration possesses a magnetic point group of
$\bar{3}m^{\prime}$. Importantly, the inversion symmetry is preserved, which
forbids the linear spin magnetoelectric response. According to Table 1, the
symmetries $C_{3}^{z}$, $C_{2}^{x}\mathcal{T}$, and $\sigma_{x}\mathcal{T}$
further enforce the following relations among elements of the nonlinear
response tensor:
$\alpha_{yxx}^{\text{int}}=\alpha_{xxy}^{\text{int}}=\alpha_{xyx}^{\text{int}}=-\alpha_{yyy}^{\text{int}}$
and $\alpha_{zxx}^{\text{int}}=\alpha_{zyy}^{\text{int}}$.

It follows that the nonlinear spin magnetoelectric response of single-layer
MnBi2Te4 is specified by only two independent elements
$\alpha_{yxx}^{\text{int}}$ and $\alpha_{zxx}^{\text{int}}$. To see this more
clearly, we assume the electric field is along an in-plane direction that
makes an angle $\theta$ from the $x$ axis, i.e.,
$\bm{E}=E(\cos\theta,\sin\theta,0)$. The induced out-of-plane nonlinear spin
polarization takes the form of $\delta s_{z}=\alpha_{zxx}^{\text{int}}E^{2}$,
which is independent of the field direction. Meanwhile, the induced in-plane
spin polarization can be expressed as supp

$(\delta s_{\parallel},\delta s_{\perp})=\alpha_{yxx}^{\text{int}}(\sin
3\theta,\cos 3\theta)E^{2},$ (11)

where $\delta s_{\parallel}$ and $\delta s_{\perp}$ denote the components
parallel and perpendicular to the electric field, respectively. It is
interesting to note that the in-plane nonlinear response is characterized by a
single element $\alpha_{yxx}^{\text{int}}$, and exhibits an angular dependence
with $2\pi/3$ periodicity. Since we are more interested in the spin
polarization that is normal to the magnetization direction, we shall focus on
$\alpha_{yxx}^{\text{int}}$ in the following discussion.

Next, we evaluate the intrinsic nonlinear spin magnetoelectric response tensor
by combining our theory with first-principles calculations (the calculation
details are presented in supp ). Figure 1(d) shows the calculated band
structure of single-layer MnBi2Te4. The system is a ferromagnetic
semiconductor with an indirect gap of 337 meV, which agrees with previous
result Chulkov2019 . We have computed all the relevant tensor elements for
$\alpha_{ij\ell}^{\text{int}}$ according to Eq. (9). The results comply with
the symmetry constraints discussed above. As mentioned, we focus on the in-
plane spin generation. The obtained $\alpha_{yxx}^{\text{int}}$ as a function
of the chemical potential is plotted in Fig. 1(e). Within the large band gap,
the value of $\alpha_{yxx}^{\text{int}}$ is small but nonzero. It is
$\sim-0.036$ $\mu_{\text{B}}$/V2, with $\mu_{\text{B}}$ as the Bohr magneton.
The response is greatly enhanced by hole doping, especially when the chemical
potential is shifted to band near degeneracy regions in the valence bands.
Because Eqs. (3) and (7) show that the BCPs are generally large around band
near degeneracies, it follows that $\alpha_{ij\ell}^{\text{int}}$, involving
integrals of BCPs, must also be peaked when the chemical potential is aligned
in such regions. Particularly, in Fig. 1(e), a peak of $375$
$\mu_{\text{B}}$/V2 is observed around $-0.1$ eV, which can be attributed to
the small gap regions marked by red arrows in Fig. 1(d). At the peak, we plot
the $k$-resolved contribution to $\alpha_{yxx}^{\text{int}}$, i.e., the
integrand in Eq. (9), in Fig. 1(f). The distribution shows an even function
with respect to both the $x$ and $y$ axis, and is peaked around the small-gap
regions.

Consider the hole doped case with the chemical potential $\sim-0.1$ eV and a
moderate applied electric field of 1 kV/cm which is readily achievable in
experiment. The induced in-plane spin magnetization in single-layer MnBi2Te4
is $\sim 0.4\times 10^{-5}$ $\mu_{\text{B}}$/nm2, which is comparable to the
reported linear spin generation in typical noncentrosymmetric systems
Mertig2016 ; Rossi2017 ; Manchon2019 . Hence, it should be detectable in
experiment, e.g., by magneto-optical Kerr spectroscopy, and it can produce
considerable spin-orbit torque effects. The effect can be further enhanced by
more than an order of magnitude at higher doping levels $\sim-0.16$ eV (shown
in supp ). In practice, the doping can be readily controlled for 2D materials
by gating.

Discussion. We have presented the first theory of nonlinear spin
magnetoelectric effect, which is the leading response in centrosymmetric
magnets. The focus here is on the intrinsic contribution, which can be
quantitatively evaluated for each material. For insulators, it captures the
total response, whereas for metals, there are additional extrinsic
contributions from the nonequilibrium distribution at the Fermi surface. The
extrinsic contributions are in principle also contained in Eq. (8), and can be
extracted by solving the distribution function, e.g., from the Boltzmann
equation. They are connected with carrier scattering and will involve the
relaxation time parameter. A systematic study of the extrinsic effect is an
interesting topic to explore in future works. In practice, the intrinsic and
extrinsic parts can be separated by their different scaling with the
relaxation time and distinct symmetry constraints, analogous to cases in
nonlinear charge transport Lai2021 ; Wang2021 ; Liu2021 ; Fu2015 ; Kang2019 ;
Du2019 .

We have demonstrated the implementation of our theory with first-principles
calculations. This will guide the experimental study and facilitate the search
for nonlinear spintronic material platforms. The effect should exist in
conventional ferromagnets like fcc Ni and Co, which preserve the inversion
symmetry. We also expect the recently fabricated 2D centrosymmetric magnets,
such as 1T-MnSe2 Kawakami2018 , CrI3 Xu2017 ; Xu2018 , and 1T-VSe2 Bonilla2018
, would be good candidates, due to their great tunability.

Finally, we note that intrinsic second order responses of other observables
that correspond to local operators, such as charge current Gao2014 , spin
current, and pseudospin, admit a similar formulation as the theory developed
here (see the Supplemental Material supp ). Therefore, our finding not only
serves as a building block for the emerging field of nonlinear spintronics,
but also forms the basis for exploring rich intrinsic nonlinear response
properties of Bloch electrons.

## References

  * (1) E. L. Ivchenko and G. E. Pikus, JETP Lett. 27, 604 (1978).
  * (2) A. G. Aronov and Y. Lyanda-Geller, JETP Lett. 50, 431 (1989).
  * (3) V. M. Edelstein, Solid State Commun. 73, 233 (1990).
  * (4) Y. K. Kato, R. C. Myers, A. C. Gossard, and D. D. Awschalom, Phys. Rev. Lett. 93, 176601 (2004).
  * (5) A. Chernyshov, M. Overby, X. Liu, J. K. Furdyna, Y. Lyanda-Geller, and L. P. Rokhinson, Nat. Phys. 5, 656 (2009).
  * (6) D. Fang, H. Kurebayashi, J. Wunderlich, K. Vyborny, L. P. Zarbo, R. P. Campion, A. Casiraghi, B. L. Gallagher, T. Jungwirth, and A. J. Ferguson, Nat. Nanotechnol. 6, 413 (2011).
  * (7) A. Manchon, J. Zelezny, I. M. Miron, T. Jungwirth, J. Sinova, A. Thiaville, K. Garello, and P. Gambardella, Rev. Mod. Phys. 91, 035004 (2019).
  * (8) D. Culcer and R. Winkler, Phys. Rev. Lett. 99, 226601 (2007).
  * (9) I. Garate and A. H. MacDonald, Phys. Rev. B 80, 134403 (2009).
  * (10) I. Garate and M. Franz, Phys. Rev. Lett. 104, 146802 (2010).
  * (11) H. Kurebayashi, J. Sinova, D. Fang, A. C. Irvine, T. D. Skinner, J. Wunderlich, V. Novak, R. P. Campion, B. L. Gallagher, E. K. Vehstedt, L. P. Zarobo, K. Vyborny, A. J. Ferguson, and T. Jungwirth, Nat. Nanotechnol. 9, 211 (2014).
  * (12) J. Zelezny, H. Gao, K. Vyborny, J. Zemen, J. Masek, A. Manchon, J. Wunderlich, J. Sinova, and T. Jungwirth, Phys. Rev. Lett. 113, 157201 (2014).
  * (13) F. Freimuth, S. Blugel, and Y. Mokrousov, Phys. Rev. B 90, 174423 (2014).
  * (14) I. M. Miron, G. Gaudin, S. Auffret, B. Rodmacq, A. Schuhl, S. Pizzini, J. Vogel, and P. Gambardella, Nat. Mater. 9, 230 (2010).
  * (15) I. M. Miron, K. Garello, G. Gaudin, P.-J. Zermatten, M. V. Costache, S. Auffret, S. Bandiera, B. Rodmacq, A. Schuhl, and P. Gambardella, Nature (London) 476, 189 (2011).
  * (16) L. Liu, C.-F. Pai, Y. Li, H.W. Tseng, D. C. Ralph, and R. A. Buhrman, Science 336, 555 (2012).
  * (17) L. Dong, C. Xiao, B. Xiong, and Q. Niu, Phys. Rev. Lett. 124, 066601 (2020).
  * (18) R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  * (19) R. Resta, Rev. Mod. Phys. 66, 899 (1994).
  * (20) D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * (21) Here we assume the system has no nontrivial Chern number.
  * (22) J.-P. Hanke, F. Freimuth, C. Niu, S. Blugel, and Y. Mokrousov, Nat. Commun. 8, 1479 (2017).
  * (23) C. Xiao, B. Xiong, and Q. Niu, Phys. Rev. B 104, 064433 (2021).
  * (24) C. Aversa and J. E. Sipe, Phys. Rev. B 52, 14636 (1995).
  * (25) R. W. Nunes and X. Gonze, Phys. Rev. B 63, 155107 (2001).
  * (26) I. Souza, J. Iniguez, and D. Vanderbilt, Phys. Rev. Lett. 89, 117602 (2002).
  * (27) Y. Gao, S. A. Yang, and Q. Niu, Phys. Rev. Lett. 112, 166601 (2014).
  * (28) Y. Gao, S. A. Yang, and Q. Niu, Phys. Rev. B 95, 165135 (2017).
  * (29) Y. Gao and D. Xiao, Phys. Rev. B 98, 060402(R) (2018).
  * (30) C. Xiao, Y. Ren, and B. Xiong, Phys. Rev. B 103, 115432 (2021).
  * (31) S. Lai, H. Liu, Z. Zhang, J. Zhao, X. Feng, N. Wang, C. Tang, Y. Liu, K. S. Novoselov, S. A. Yang, and W. B. Gao, Nat. Nanotechnol. 16, 869 (2021).
  * (32) C. Wang, Y. Gao, and D. Xiao, Phys. Rev. Lett. 127, 277201 (2021).
  * (33) H. Liu, J. Zhao, Y.-X. Huang, W. Wu, X.-L. Sheng, C. Xiao, and S. A. Yang, Phys. Rev. Lett. 127, 277202 (2021).
  * (34) H. Liu, J. Zhao, Y.-X. Huang, X. Feng, C. Xiao, W. Wu, S. Lai, W. B. Gao, S. A. Yang, arXiv:2106.04931
  * (35) Y. Gao, S. A. Yang, and Q. Niu, Phys. Rev. B 91, 214405 (2015).
  * (36) See Supplemental Material for details about the semiclassical nonlinear response formulation applicable not only to spin, the derivation of Eq. (11), the computational method and the spin response within a wider range of chemical potential in MnBi2Te4.
  * (37) Y. Gong, J. Guo, J. Li, K. Zhu, M. Liao, X. Liu, Q. Zhang, L. Gu, L. Tang, X. Feng, D. Zhang, W. Li, C. Song, L. Wang, P. Yu, X. Chen, Y. Wang, H. Yao, W. Duan, Y. Xu $et$ $al.$, Chin. Phys. Lett. 36, 076801 (2019).
  * (38) M. M. Otrokov, I. P. Rusinov, M. Blanco-Rey, M. Hoffmann, A. Yu. Vyazovskaya, S. V. Eremeev, A. Ernst, P. M. Echenique, A. Arnau, and E. V. Chulkov, Phys. Rev. Lett. 122, 107202 (2019).
  * (39) D. Zhang, M. Shi, T. Zhu, D. Xing, H. Zhang, and J. Wang, Phys. Rev. Lett. 122, 206401 (2019).
  * (40) J.-H. Li, Y. Li, S. Du, Z. Wang, B.-L. Gu, S.-C. Zhang, K. He, W. Duan, and Y. Xu, Sci. Adv. 5, eaaw5685 (2019).
  * (41) M. M. Otrokov $et$ $al.$, Nature (London) 576, 416 (2019).
  * (42) A. Zeugner $et$ $al.$, Chem. Mater. 31, 2795 (2019).
  * (43) J.-Q. Yan, Q. Zhang, T. Heitmann, Z. Huang, K. Y. Chen, J.-G. Cheng, W. Wu, D. Vaknin, B. C. Sales, and R. J. McQueeney, Phys. Rev. Mater. 3, 064202 (2019).
  * (44) S. H. Lee, Y. Zhu, Y. Wang, L. Miao, T. Pillsbury, H. Yi, S. Kempinger, J. Hu, C. A. Heikes, P. Quarterman, W. Ratcliff, J. A. Borchers, H. Zhang, X. Ke, D. Graf, N. Alem, C.-Z. Chang, N. Samarth, and Z. Mao, Phys. Rev. Res. 1, 012011(R) (2019).
  * (45) J. Cui, M. Shi, H. Wang, F. Yu, T. Wu, X. Luo, J. Ying, and X. Chen, Phys. Rev. B 99, 155125 (2019).
  * (46) Y. Deng, Y. Yu, M. Z. Shi, Z. Guo, Z. Xu, J. Wang, X. H. Chen, and Y. Zhang, Science 367, 895 (2020).
  * (47) C. Liu, Y. Wang, H. Li, Y. Wu, Y. Li, J. Li, K. He, Y. Xu, J. Zhang, and Y. Wang, Nat. Mater. 19, 522 (2020).
  * (48) A. Johansson, J. Henk, and I. Mertig, Phys. Rev. B 93, 195440 (2016).
  * (49) M. Rodriguez-Vega, G. Schwiete, J. Sinova, and E. Rossi, Phys. Rev. B 96, 235419 (2017).
  * (50) I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * (51) K. Kang, T. Li, E. Sohn, J. Shan, and K. F. Mak, Nat. Mater. 18, 324 (2019).
  * (52) Z. Z. Du, C. M. Wang, S. Li, H.-Z. Lu, and X. C. Xie, Nat. Commun. 10, 3047 (2019).
  * (53) D. J. O’Hara, T. Zhu, A. H. Trout, A. S. Ahmed, Y. K. Luo, C. H. Lee, M. R. Brenner, S. Rajan, J. A. Gupta, D. W. McComb, and R. K. Kawakami, Nano Lett. 18, 3125 (2018).
  * (54) B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Nature (London) 546, 270 (2017).
  * (55) K. L. Seyler, D. Zhong, D. R. Klein, S. Gao, X. Zhang, B. Huang, E. Navarro-Moratalla, L. Yang, D. H. Cobden, M. A. McGuire, W. Yao, D. Xiao, P Jarillo-Herrero, and X. Xu, Nat. Phys. 14, 277 (2018).
  * (56) M. Bonilla, S. Kolekar, Y. Ma, H. C. Diaz, V. Kalappattil, R. Das, T. Eggers, H. R. Gutierrez, M.-H. Phan and M. Batzill, Nat. Nanotechnol. 13, 289 (2018).