# Inference in generalized linear models with robustness to misspecified
variances

Riccardo De Santis111To contact<EMAIL_ADDRESS>University
of Padova, Italy  Jelle J. Goeman  Leiden University Medical Center, The
Netherlands  Jesse Hemerik  Wageningen University and Research, The
Netherlands  Livio Finos  University of Padova, Italy

###### Abstract

Generalized linear models usually assume a common dispersion parameter, an
assumption that is seldom true in practice. Consequently, standard parametric
methods may suffer appreciable loss of type I error control. As an
alternative, we present a semi-parametric group invariance method based on
sign flipping of score contributions. Our method requires only the correct
specification of the mean model, but is robust against any misspecification of
the variance. We present tests for single as well as multiple regression
coefficients. The method is available in the R library flipscores.

##  1 Introduction

When testing for equality of means between two groups, statisticians often
worry that there may be a difference in variance between the two groups. The
task of properly taking this variance into account, known as the Behrens-
Fisher problem (Fisher, 1935, 1941), has generated a huge statistical
literature (Kim and Cohen, 1998; Chang and Pal, 2008). In the linear
regression model, which generalizes the two-group model, the assumption that
all error terms have equal variance is less often questioned. While goodness-
of-fit tests and other diagnostics for heteroskedasticity exist (Goldfeld and
Quandt, 1965; Glejser, 1969; Breusch and Pagan, 1979; Jarque and Bera, 1980;
Cook and Weisberg, 1983; Long and Ervin, 2000), there is no simple and general
way to follow up on obvious lack of fit (Rochon et al., 2012). Still, the
problems arising from misspecified variances in regression can be as severe in
regression as they are the two-group model, especially when the variance
depends on the covariates in the model in an unknown way: standard errors are
too large or too small, and statistical tests can become severely conservative
or anti-conservative.

The situation is worse if we broaden the perspective to generalized linear
models (GLMs). Overdispersed GLMs, that allow for additional variation between
subjects, have become important in many fields; e.g. negative binomial or
quasi-Poisson models in RNA sequencing (Love et al., 2014). GLMs and
overdispersed GLMs generally use a single dispersion parameter for all
subjects. This assumption is similar to the common variance assumption in the
linear model, and violation of the assumption can lead to loss of Type I error
control in the same way. The diagnostics to detect a failure of this variance
model are more complex in GLMs than in the linear model, since the variance of
residuals differs between observations even in a correctly specified model. A
general approach that is robust to variance misspecification is to use the
Sandwich estimator (Eicker, 1967; Huber, 1967; White, 1982). However, as we
will see in the simulations of section 7, this is not always an optimal
choice, since it often performs poorly for small sample sizes (Boos, 1992;
Kauermann and Carroll, 2000; Maas and Hox, 2004; Freedman, 2006).

In this paper we will develop a novel method for inference on regression
coefficients in GLMs that is robust to misspecification of the variance. Our
method is an extension of the method of Hemerik et al. (2020), which compares
the score to a semi-parametric reference distribution calculated by randomly
sign-flipping each subject’s contribution to the score. This method was known
to have robustness to misspecification of the variance by a common
multiplicative constant only. We extend this method in two ways. First, we
standardize each sign flipping of the scores to achieve faster convergence to
the nominal Type I error level. Second, we show that the new standardized test
remains asymptotically valid under any misspecification of the variance. The
resulting test gives valid inference without concern for heteroskedasticity.
We find that the convergence to the nominal level is comparable, and sometimes
even faster than a parametric test when both methods are asympotically valid,
and power is also comparable. Our results actually imply that the same
robustness holds for the original test of Hemerik et al. (2020), but the new
test achieves better error control than the original in small samples.

In the more limited context of linear models, other semi-parametric approaches
have been proposed before. Winkler et al. (2014) gives a systematic comparison
of such methods for linear regression in a realistic setting. Our novel method
has a connection to the Wild Bootstrap (Davidson and Flachaire, 2008), which
addresses the problem of heteroskedasticity only in linear regression. Our
proposed method can be seen as an extension of the Wild Bootstrap to
generalized linear models and to situations with nuisance parameters.

##  2 Generalized linear models

Assume that we observe $n$ independent observations
$y=\\{y_{1},\dots,y_{n}\\}^{T}$ from the exponential dispersion family, i.e. a
density of the form (Agresti, 2015)

$f(y_{i};\theta_{i},\phi_{i})=\exp\left\\{\frac{y_{i}\theta_{i}-b(\theta_{i})}{a(\phi_{i})}+c(y_{i},\phi_{i})\right\\},$

where $\theta_{i}$ and $\phi_{i}$ are respectively the canonical and the
dispersion parameter. We will assume throughout the paper that the model
chosen fulfills the usual regularity conditions (Azzalini, 1996, Chapter 3).
We derive mean and variance of the observed outcome, respectively, as

$\mu_{i}=E(y_{i})=b^{\prime}(\theta_{i});\qquad\mbox{var}(y_{i})=b^{\prime\prime}(\theta_{i})a(\phi_{i}),$

where primes denote derivatives. We will consider through the paper, without
loss of generality, $a(\phi_{i})=\phi_{i}$. We assume that the mean of $y_{i}$
depends on observed covariates $(x_{i},z_{i})$ through the following relation,
written in vector-matrix form as

$g(\mu)=\eta=X\beta+Z\gamma$

where $\mu=(\mu_{1},\dots,\mu_{n})^{T}$ denotes the mean vector, $g(\cdot)$ is
the link function, taken to operate elementwise on a vector, $(X,Z)$ is the
full rank design matrix with $\mathrm{dim}(X)=n\times 1$,
$\mathrm{dim}(Z)=n\times q$, where $q$ does not depend on $n$, and
$(\beta,\gamma)$ are unknown parameters. We consider $\beta$ the parameter of
interest, and $\gamma$ a nuisance parameter.

Let $\ell(\beta,\gamma)$ denote the log-likelihood function, from which we can
derive the score vector with elements

$\frac{\partial\ell(\beta,\gamma)}{\partial\beta}=s_{\beta}=X^{T}DV^{-1}(y-\mu);\qquad\frac{\partial\ell(\beta,\gamma)}{\partial\gamma}=s_{\gamma}=Z^{T}DV^{-1}(y-\mu),$

where, in a compact matrix form, we have

$D=\mathrm{diag}\left\\{\frac{\partial\mu_{i}}{\partial\eta_{i}}\right\\};\qquad
V=\mathrm{diag}\\{\mathrm{\mbox{var}}(y_{i})\\}.$

Taking the derivatives of the score vector we obtain the Fisher information
matrix

$\mathcal{I}=\begin{pmatrix}\mathcal{I}_{\beta,\beta}&\mathcal{I}_{\beta,\gamma}\\\
\mathcal{I}_{\gamma,\beta}&\mathcal{I}_{\gamma,\gamma}\end{pmatrix}=\begin{pmatrix}X^{T}WX&X^{T}WZ\\\
Z^{T}WX&Z^{T}WZ\end{pmatrix}$

where $W=DV^{-1}D$. Denote by $d_{i}$, $v_{i}$, and $w_{i}$ the $i$-th
diagonal elements of $D$, $V$ and $W$, respectively. Note that all the
matrices defined here and in the rest of the paper depend on $n$. We will
suppress this dependence in the notation. The reader may assume that all
quantities depend on $n$, except when explicitly stated otherwise.

The model as defined above has separate dispersion parameters $\phi_{i}$ for
$i=1,\ldots,n$. These parameters cannot be consistently estimated unless they
satisfy some restrictions (Neyman and Scott, 1948). Thus, in practice, the
analyst will use some putative model, which imposes such restrictions. This
leads to computable estimates ${\hat{\phi}}_{1},\dots,{\hat{\phi}}_{n}$, which
will generally be inconsistent. Remarkably, the incorrect estimation of the
dispersion parameters does not affect the consistency of the estimates of the
regression coefficients, as long as the mean (and hence the link function) are
well-specified (Agresti, 2015). The corresponding standard errors, however,
are no longer reliable.

As an example of the set-up, consider the situation that the data are
generated according to a negative binomial model with a log link function and
a different and unknown dispersion parameter for every observation. The
analyst could estimate these dispersion parameters under the additional
constraint that $\phi_{i}=\phi$ for all $i$, or even that $\phi_{i}=1$ for all
$i$. Alternatively, the dispersion parameters could be modeled as a function
of the covariates (McCullagh and Nelder, 1989, Chapter 10). Such strategies
will lead to inconsistent estimates of the dispersion parameters unless the
true model happens to fulfil the chosen constraints (White, 1982).

About the model and the estimation strategy we make the following assumptions.
First, we demand that the regression parameters $\beta$ and $\gamma$ are
consistently estimated, and the true variances to be uniformly bounded.
Secondly, we assume that the estimates of $\phi_{1},\phi_{2},\ldots$, though
not consistent, will converge, as stated more formally below.

######  Assumption 1.

For $i=1,2,\ldots$, we have $\hat{\phi}_{i}-\tilde{\phi}_{i}\to 0$ in
probability as $n\to\infty$, where $K_{1}<\tilde{\phi}_{i}\leq K_{2}$, and
$K_{1},K_{2}$ are strictly positive constants not depending on $n$.

It is known (Huber, 1967) that the use of maximum likelihood estimation in a
misspecified model leads under minimal conditions to a well-defined limit of
the estimator. In this sense, it is possible to define a “true” density, even
in a misspecified model, which is intended to be the closest to the true
density which generates the data in terms of Kullback-Leibler distance (White,
1982). This implies that Assumption 1 generally holds in situations where we
estimate $\phi_{1},\ldots,\phi_{n}$ using a restricted model.

We place a tilde symbol on all the quantities obtained when plugging the
limits $\tilde{\phi}_{1},\ldots,\tilde{\phi}_{n}$ into the model. In
particular, let $\tilde{V},\tilde{W}$ be the variance and weight matrices
obtained by fixing $\phi_{i}=\tilde{\phi}_{i}$. Note that the putative model
only misspecifies the variance, but not the mean: the mean vector and the link
function remain correctly specified and can still be consistently estimated
(Agresti, 2015).

We further assume the following assumption related to the Fisher information
given by the quantities $\tilde{\phi}_{1},\dots,\tilde{\phi}_{n}$. This is a
standard condition in regular models (Van der Vaart, 1998).

######  Assumption 2.

Let $\tilde{\mathcal{I}}$ be the Fisher information matrix for $W=\tilde{W}$.
The $\lim_{n\xrightarrow{}\infty}n^{-1}\tilde{\mathcal{I}}_{\beta,\beta}$
converges to some positive constant.

Further, the following mild condition is needed to apply the multivariate
central limit theorem (Billingsley, 1986, Chapter 5) within the framework of
Hemerik et al. (2020), that we will apply. It is needed to avoid pathological
cases, such as vanishing or dominating observations.

######  Assumption 3.

Let

$\tilde{\nu}_{i,\beta}=\frac{(y_{i}-\hat{\mu}_{i})x_{i}d_{i}}{\tilde{v}_{i}};\qquad\tilde{\nu}_{i,\gamma}=\frac{(y_{i}-\hat{\mu}_{i})z_{i}d_{i}}{\tilde{v}_{i}};$

define the element-wise score contribution and

$\tilde{\nu}^{*}_{i,\beta}=\tilde{\nu}_{i,\beta}-\tilde{\mathcal{I}}_{\beta,\gamma}\tilde{\mathcal{I}}_{\gamma,\gamma}\tilde{\nu}_{i,\gamma}.$

We require, for all $\epsilon>0$,

$\lim_{n\xrightarrow{}\infty}\frac{1}{n}\sum_{i=1}^{n}E[(\tilde{\nu}_{i,\beta}^{*})^{2}\mathbf{1}_{\\{|\tilde{\nu}_{i,\beta}^{*}|/\sqrt{n}>\epsilon\\}}]\xrightarrow{}0$

where $\mathbf{1}_{\\{\dots\\}}$ is the indicator function, and
$n^{-1}\sum_{i=1}^{n}\mbox{var}(\tilde{\nu}_{i,\beta}^{*})$ to converge to
some positive constant.

When the model is correctly specified the last condition implies Assumption 2.
In case of misspecified models we do not have this relation, since the so
called “information identity” does not hold (Azzalini, 1996).

Throughout the paper, we will make a distinction between two situations: the
situation that the model is correctly specified, i.e. $\tilde{W}=W$, and the
putative model is true; or the general case that $\tilde{W}\neq W$.

##  3 Sign-flipping effective score test

In the model described in the previous section we are interested in testing
the following null hypothesis

$\displaystyle H_{0}:\ \beta=\beta_{0}\ |\
(\gamma,\phi_{1},\dots,\phi_{n})\in\Gamma\times\Phi\times\dots\times\Phi$ (1)

against a one or two-sided alternative, where $\Gamma\subseteq\mathbb{R}^{q}$
and $\Phi\subseteq(0,\infty)$. The formulation of the null hypothesis above
makes explicit that only the target parameter $\beta$ is fixed under $H_{0}$,
but the nuisance parameters are unconstrained.

Since $\phi_{1},\dots,\phi_{n}$, as remarked, are difficult to estimate, we
require a test that is robust to misspecification of these parameters. Hemerik
et al. (2020) proposed a general semi-parametric test that has robustness to
misspecification of the variance by a constant, and presented promising
simulation results suggesting more general robustness properties. We will
start from this test.

Hemerik et al. (2020) uses the effective score for $\beta$ as the test
statistic. In the model (2) the effective score is

$S=s_{\beta}-\mathcal{I}_{\beta,\gamma}\mathcal{I}_{\gamma,\gamma}^{-1}s_{\gamma}.$

In the context of generalized linear models, the statistic can be written as

$S=n^{-1/2}X^{T}W^{1/2}(I-H)V^{-1/2}(y-{\hat{\mu}})$

where

$H=W^{1/2}Z(Z^{T}WZ)^{-1}Z^{T}W^{1/2},$ (2)

is the hat matrix, and $\hat{\mu}$ is the vector of fitted values of the model
under the null hypothesis. Note that, since $S$ is an inner product of the two
$n$-vectors $n^{-1/2}V^{-1/2}(I-H)W^{1/2}X$ and $y-{\hat{\mu}}$, it can be
written as a sum of $n$ terms, which we call the effective score
contributions.

To calculate the critical value, Hemerik et al. (2020) proposed to use sign
flips of the effective score contributions, randomly multiplying each score
contribution by $-1$ or $1$. In matrix notation these sign flips can be
represented by a random diagonal matrix $\mathcal{F}$ of dimension $n$, whose
non-zero elements are independent random variables that take values $-1$ and 1
with equal probability. Consequently, the effective sign-flip score statistic
for a given flip matrix $\mathcal{F}=F$ is defined as

$S(F)=n^{-1/2}X^{T}W^{1/2}(I-H)V^{-1/2}F(y-{\hat{\mu}}).$ (3)

Note that for $\mathcal{F}=\mathbf{I}$ (the identity matrix) we recover the
observed effective score.

An asymptotic $\alpha$ level test is then derived as follows. First, Hemerik
et al. (2020) prove (asymptotic) invariance of the first two moments of the
test statistic under the action of $\mathcal{F}$, i.e. that
$E\\{S(\mathcal{F})\\}-E\\{S(\mathbf{I})\\}=0$ and
$\mbox{var}\\{S(\mathcal{F})\\}-\mbox{var}\\{S(\mathbf{I})\\}\to 0$ as
$n\to\infty$. Next, they apply the Lindberg-Feller multivariate central limit
theorem to show that, for independently drawn flip matrices
$F_{2},\ldots,F_{g}$, where $g$ does not depend on $n$, the vector
$S(\mathbf{I}),S(F_{2}),\ldots,S(F_{g})$ converges to a vector of independent
and identically distributed random variables. By Lemma 1 of Hemerik et al.
(2020) is then possible to obtain an asymptotic $\alpha$ level test for the
null hypothesis (1) against a one or two-sided alternative. Assume we observe
values $S_{1},\dots,S_{g}$, where $S_{1}=S(\mathbf{I})$ is the observed test
statistic, while the sorted values are $S_{(1)}\leq\dots\leq S_{(g)}$. Without
loss of generality, consider testing (1) against $H_{1}:\beta>\beta_{0}$. The
test that they consider rejects the null hypothesis if

$S_{1}>S_{(\lceil(1-\alpha)g\rceil)}$ (4)

where $\lceil\cdot\rceil$ represents the ceiling function. Analogous
procedures can be defined for $H_{1}:\beta<\beta_{0}$ or $\beta\neq\beta_{0}$
(Hemerik et al., 2020).

The definition of the test involves $W$, which is unknown. In practice we only
have an estimate of $W$ available, which converges to $\tilde{W}$ by
Assumption 1. For the theoretical results of the remainder of this paper, we
will treat $\tilde{W}$, though not $W$, as known. To motivate this, we note,
in the first place, that any error terms relating to estimation of $\tilde{W}$
are of lower asymptotic order (Barndorff-Nielsen and Cox, 1994). Secondly, as
we will show below in Section 5, the validity of the test is invariant to
misspecification of $\tilde{W}$, so it is not important to estimate it
correctly.

##  4 Standardized sign-flips score test

The test based on the effective score is proven to be asymptotically exact by
Hemerik et al. (2020), but simulations in that paper show an anti-conservative
behavior for small sample sizes. This problem relates to a lack of null-
invariance of the distribution of the test statistic. In this section we
revisit the concept of null-invariance in some detail. This discussion will
motivate our modification of the test of Hemerik et al. (2020).

Tests based on random data transformations such as sign flips or permutations
can have exact control of the type I error rate if the observed test statistic
$S(\mathbf{I})$ and the randomly flipped $S(\mathcal{F})$ have the same
distribution. If this property holds, we say that the transformations are
null-invariant. Hemerik and Goeman (2018) showed how null-invariance can be
used to construct tests with exact control of the type I error. Null-
invariance can be achieved in linear models (Solari et al., 2014; Huh and
Jhun, 2001; Kherad-Pajouh and Renaud, 2010) or in GLMs with simple
experimental designs. With more complex designs, ad hoc solutions sometimes
exist, e.g. when all covariates are discrete (Pesarin, 2001). For general
GLMs, null-invariant transformations do not exist.

In the case of the test of Hemerik et al. (2020), null-invariance is therefore
replaced by an equality in distribution that holds only asymptotically. We
restate the result concerning the validity of that test below for the special
case of GLMs. Our alternative proof emphasizes the importance of asymptotic
null-invariance.

######  Theorem 1.

Assume that the variances are correctly specified, that is, $\tilde{V}=V$, and
that Assumption 1-3 hold. For $n\xrightarrow{}\infty$, the test that rejects
$H_{0}$ if (4) holds is an asymptotically $\alpha$ level test.

###### Proof.

This, and all the other proofs, are given in the Appendix. ∎

If null-invariance holds only asymptotically, the distributions of
$S(\mathbf{I})$ and $S(\mathcal{F})$ are not identical, but both converge to
the same distribution. It is natural to suppose that the more aspects of the
distributions of $S(\mathbf{I})$ and $S(\mathcal{F})$ are equal to each other
in all finite samples, the closer the resulting test is to an exact test. This
intuition motivated (Commenges, 2003; Solari et al., 2014) to formulate the
concept of second-moment null-invariance. A random transformation
$\mathcal{F}$ is second-moment null-invariant if $S(\mathbf{I})$ and
$S(\mathcal{F})$ have identical first and second moment. For second-moment
null-invariant transformations, asymptotics is only needed for the convergence
of the higher-order moments.

The test of Hemerik et al. (2020) does not have the second-moment null-
invariance property: though the first moment of $S(\mathbf{I})$ and
$S(\mathcal{F})$ is equal in finite samples, the variance is not. In fact, we
prove in Theorem 2 below that for finite samples $S(\mathbf{I})$ always has a
larger variance than $S(\mathcal{F})$, at least in the linear model. As a
consequence, the test shows a tendency to anti-conservativeness in finite
samples. Since the variance of the observed test statistic is always larger
than its flipped counterpart, extreme values in $S(\mathbf{I})$ are more
probable than in the reference distribution of $S(\mathcal{F})$. Consequently,
the test tends to reject the null hypothesis too often.

######  Proposition 1.

Consider a normal regression model with identity link. Assume that the
variances are correctly specified, that is, $\tilde{V}=V$, and that Assumption
1-3 hold. For finite sample size, the effective sign-flip score statistic
defined as in (3) has
$\mbox{var}\\{S(\mathbf{I})\\}>\mbox{var}\\{S(\mathcal{F})\\}$.

Proposition 1 is formulated for the linear model only. For the general case
the same positive difference between the variance of $S(\mathbf{I})$ and
$S(\mathcal{F})$ is present, but it is not the only term in the asymptotic
expansion. However, we see a tendency to small sample anti-conservativeness in
all simulations, as was also observed by Hemerik et al. (2020).

The concept of second-moment null-invariance suggests a simple but efficient
way to remedy the anti-conservativeness of the test of Hemerik et al. (2020)
by studentization of the flipped effective score test statistic. We call the
resulting test the standardized sign-flip score statistic, defined for a given
flip $\mathcal{F}=F$ as

$S^{*}(F)=S(F)/\mbox{var}\\{S(F)\\}^{1/2}.$ (5)

The expression

$\mbox{var}\\{S(F)\\}=n^{-1}X^{T}W^{1/2}(I-H)F(I-H)F(I-H)W^{1/2}X+o_{p}(1)$

is easy to estimate. This estimate can be calculated for each $F$ in linear
time in $n$, as we show in Lemma 5 in the Appendix.

This simple transformation has great impact and it substantially improves the
type I error control of the previous test. Analogously with the test defined
in (4), consider without loss of generality testing (1) against
$H_{1}:\beta>\beta_{0}$. Take the test that rejects the null hypothesis if

$S^{*}_{1}>S^{*}_{(\lceil(1-\alpha)g\rceil)}.$ (6)

The following Proposition concerns the validity of the test.

######  Proposition 2.

Assume that the variances are correctly specified, that is, $\tilde{V}=V$, and
that Assumptions 1-3 hold. The standardized sign-flip score statistic is
finite sample second-moment null-invariant. The test that rejects $H_{0}$ if
(6) holds is i) an exact $\alpha$ level test in case of normal response with
identity link and ii) asymptotically exact in all other cases.

##  5 Robustness of the test

Hemerik et al. (2020) showed promising simulation results suggesting
robustness against variance misspecification, but a formal proof was given
only for the special case that all variances were misspecified by the same
multiplicative constant. In this section we prove general robustness of the
test to variance misspecification, both for the novel test derived in the
previous section and for the original version. This is the most important
result of this paper.

Wrong specification of the model variance means that

$E\left\\{(y-\mu)(y-\mu)^{T}\right\\}=V\neq\tilde{V}.$

We first revisit the simple case considered by Hemerik et al. (2020), in which
the variance is misspecified by a multiplicative constant. In this case, the
properties of the standardized sign-flip score test are not affected by this
misspecification. In particular, we have second-moment null-invariance.

######  Proposition 3.

Assume that Assumptions 1-3 hold. If the variances are misspecified by any
finite constant $c>0$, that is, $V=c\tilde{V}$, the standardized sign-flip
score statistic is second-moment null-invariant. The test that rejects $H_{0}$
if (6) holds is i) an exact $\alpha$ level test in case of normal response
with identity link and ii) asymptotically exact in all other cases.

The proposition uses the property that the test is invariant to multiplication
by a constant. This result is relevant, for instance, in models with a common
unknown dispersion parameter, such as normal regression model. In such cases
the test can be performed, and retains second-moment null-invariance, without
the need to estimate the common dispersion parameter. We can save ourselves
the effort of estimating it, taking simply $\hat{\phi}=1$. This special
situation coincides with the standard parametric framework based on the quasi-
likelihood approach.

The main result of this paper concerns the case of general variance
misspecification. The following theorem, which is a strong improvement over
the properties shown by Hemerik et al. (2020), shows that we can still get an
asymptotic exact test.

######  Theorem 2.

Assume that the variances are misspecified, that is, $V\neq\tilde{V}$ and that
Assumptions 1-3 hold. For $n\xrightarrow{}\infty$, the standardized sign-flip
score statistic is asymptotically second-moment null-invariant. The test that
rejects $H_{0}$ if (6) holds is an asymptotically $\alpha$ level test.

The following corollary enlarges the robustness property to the effective
sign-flip score test, as an immediate consequence from the proof of the
preceding theorem.

######  Corollary 1.

Assume that the variances are misspecified, that is, $V\neq\tilde{V}$, and
that Assumptions 1-3 hold. For $n\xrightarrow{}\infty$, the effective sign-
flip score statistic is asymptotically second-moment null-invariant. The test
that rejects $H_{0}$ if (4) holds is an asymptotically $\alpha$ level test.

##  6 Multivariate test

Until now we have considered hypotheses about a single parameter
$\beta\in\mathbb{R}$. We now generalize the test of the previous section to
$\beta\in\mathbb{R}^{d},\;d<n-q$. We consider a standard asymptotic setting
where $d$ is fixed while $n$ increases. The null hypothesis of interest is now
given by:

$H_{0}:\beta=\beta_{0}\in\mathbb{R}^{d}|\
(\gamma,\phi_{1},\dots,\phi_{n})\in\Gamma\times\Phi\times\dots\times\Phi,$

where $\Gamma\subseteq\mathbb{R}^{q}$ and $\Phi\subseteq(0,\infty)$, which
reduces to the null hypothesis (1) if $d=1$.

For this multivariate setting, we have to generalize the assumptions of
Section 2. Assumption 1 remains unchanged while Assumptions 2 and 3 are
replaced by their multivariate counterparts, respectively,

######  Assumption 4.

The $\lim_{n\xrightarrow{}\infty}n^{-1}\tilde{\mathcal{I}}_{\beta,\beta}$
converges to a positive definite matrix.

######  Assumption 5.

We require, for all $\epsilon>0$,

$\lim_{n\xrightarrow{}\infty}\frac{1}{n}\sum_{i=1}^{n}E[||\tilde{\nu}_{i,\beta}^{*}||^{2}\mathbf{1}_{\\{||\tilde{\nu}_{i,\beta}^{*}||/\sqrt{n}>\epsilon\\}}]\xrightarrow{}\mathbf{0},$

and $n^{-1}\sum_{i=1}^{n}\mbox{var}(\tilde{\nu}_{i,\beta}^{*})$ to converge to
a positive definite matrix, where $||\cdot||$ denotes the $\ell^{2}$ norm and
$\mathbf{0}$ is a $d$-dimensional zero vector.

Hemerik et al. (2020, Section 4) derived a generalization of the sign-flipped
effective score test statistic, as follows. Noting that the effective score is
now a $d$-dimensional vector $S(F)=(S^{1}(F),\dots,S^{d}(F))^{T}$, an
asymptotically exact $\alpha$-level test can be constructed by using the idea
of the nonparametric combination methodology (Pesarin, 2001). The test
statistic takes the form

$T(F)=\left\\{S(F)\right\\}^{T}M\left\\{S(F)\right\\}$

where $M$ is any non-zero matrix. Usually $M$ is chosen to be a symmetric
matrix, and in general this choice influences the distribution of the power
between the alternatives (see Hemerik et al. (2020) for details). The inverse
of an estimate of the effective Fisher information of $\beta$, if available,
is a common choice for $M$.

As in Section 4, we can improve control of type I error by standardizing the
score vector. Indeed, the same reasoning of Theorem 2 applies, showing that
the sign-flip standardized score vector is finite sample second-moment null-
invariant. The definition of the test is analogous by noting that now

$\mbox{var}\\{S(F)\\}=n^{-1}X^{T}W^{1/2}(I-H)F(I-H)F(I-H)W^{1/2}X+o_{p}(1)$

is a $d\times d$ matrix and hence the standardized score

$S^{*}(F)=S(F)/\mbox{var}\\{S(F)\\}^{1/2}$

is a $d$-dimensional vector. We can therefore define the test statistic as

$T^{*}(F)=\left\\{S^{*}(F)\right\\}^{T}M\left\\{S^{*}(F)\right\\}.$

Assume we observe values $T^{*}_{1},\dots,T^{*}_{g}$, where
$T^{*}_{1}=T^{*}(\mathbf{I})$ is the observed test statistic, while the sorted
values are $T^{*}_{(1)}\leq\dots\leq T^{*}_{(g)}$. Consider the test that
rejects the null hypothesis if

$T^{*}_{1}>T^{*}_{(\lceil(1-\alpha)g\rceil)}.$ (7)

The following Proposition states properties regarding null-invariance and
asymptotic exactness.

######  Proposition 4.

Assume that the variances are correctly specified, that is, $\tilde{V}=V$, and
that Assumptions 1, 4 and 5 hold. The $d$-dimensional standardized sign-flip
score vector is finite sample second-moment null-invariant. The test that
rejects $H_{0}$ if (7) holds is an asymptotically $\alpha$ level test.

Finally, the following theorem shows that the robustness property is inherited
by this multivariate extension.

######  Theorem 3.

Assume that the variances are misspecified, that is, $V\neq\tilde{V}$ and that
Assumptions 1,4,5 hold. For $n\xrightarrow{}\infty$, the $d$-dimensional
standardized sign-flip score vector is asymptotically second-moment null-
invariant. The test that rejects $H_{0}$ if (7) holds is an asymptotically
$\alpha$ level test.

##  7 Simulation study

Figure 1: Type I error control. (  , standardized;  , effective;  ,
parametric;  , sandwich). Figure 2: Power comparison (  , standardized;  ,
effective;  , parametric;  , sandwich).

We explore six different settings to compare empirically the type I error
control of the usual parametric approach (by considering the Wald test), the
effective and standardized flip-score tests and the Wald test based on the use
of the sandwich estimator of the variance to correct for variance
misspecification. A total of 5000 simulations have been carried out for each
setting. The covariates have been drown from a multivariate normal
distribution, with $X\in\mathbb{R}$ and $Z\in\mathbb{R}^{3}$. The three
nuisance covariates have correlation with the target variable equal to
$(0.5,0.1,0.1)$, while the true parameter is set to $\beta=0$. The null
hypothesis considered is $H_{0}:\beta=0$ against a two-sided alternative, with
a nominal significance level $\alpha=0.05$. Different sample sizes are
considered $(25,50,100,200,500,1000)$.

The top two plots of Figure 1 are settings with correctly specified models,
respectively a Poisson and a Logistic regression models. The middle two plots
represent normal models with neglected heteroscedasticity, which depends
either on a nuisance covariate or on the tested variable, i.e. respectively
$\mbox{var}(y_{i})=4z^{2}_{i}$ and $\mbox{var}(y_{i})=4x^{2}_{i}$. In the
bottom-left plot of Figure 1 a Poisson model was fitted when the true
distribution was negative binomial with dispersion parameter $\phi=1$, that
is, additional heteroscedasticity relative to the Poisson model that depends
on the mean. The bottom-right plot displays the results of a two-sample test
fitting a negative binomial. A common dispersion parameter was assumed for the
fitted negative binomial model, but the two groups are unbalanced (proportion
equal to 2/3 and 1/3), and generated from two distributions with different
dispersion parameter ($0.4$ and 1).

In the two top plots we observe that, in case of correctly specified models,
the standardized test is close to the nominal level, even with $n=25$,
improving the asymptotic convergence of the effective test. The parametric
test shows few rejections with small sample size for the logistic model, due
to a poor approximation of the likelihood. The sandwich shows slow convergence
in both cases, which is unsatisfactory since we are dealing with a well-
specified model. The middle-left plot has similar behavior.

The last three plots show the failure of the parametric test in presence of
some forms of variance misspecification, where the rejection fraction
converges to a level far from the nominal for increasing sample size. The
standardized test outperforms its competitors in all cases, being closest to
the nominal level. In particular, the improvement over the effective test and
further over the sandwich test is clear.

In the final two plots it is remarkable to see that the standardized test
seems to converge to the effective test faster than going to the nominal
level, leading to an initial worsening of the true level, with a recovery for
larger sample size.

Finally, Figure 2 contains an evaluation of the power of the tests with two
well-specified models, respectively a Poisson and a Normal model, with true
parameter equal respectively to 0.3 and 1. The results for the sandwich
estimator are given for completeness, although they are not comparable for
small sample size, since that method has no control of the type I error. We
see that the improvement of type I error control of the standardized test with
respect to the effective naturally costs some power. Analogously, we see some
power loss also with respect to the parametric model as expected, but this
difference is remarkably small.

## References

  * Agresti (2015) Alan Agresti.  _Foundations of Linear and Generalized Linear Models_.  Wiley, New York, 2015.
  * Azzalini (1996) Adelchi Azzalini.  _Statistical Inference based on the likelihood_.  Chapman and Hall, Boca Raton, FL, 1996.
  * Barndorff-Nielsen and Cox (1994) Ole E. Barndorff-Nielsen and David R. Cox.  _Inference and Asymptotics_.  Chapman and Hall, Boca Raton, FL, 1994.
  * Billingsley (1986) Patrick Billingsley.  _Probability and Measure_.  Wiley, New York, 1986.
  * Boos (1992) Dennis D. Boos.  On generalized score tests.  _The American Statistician_ , 46(4):327–333, 1992.
  * Breusch and Pagan (1979) Trevor S. Breusch and Adrian R. Pagan.  A simple test for heteroscedasticity and random coefficient variation.  _Econometrica_ , 47(5):1287–1294, 1979.
  * Chang and Pal (2008) Ching-Hui Chang and Nabendu Pal.  A revisit to the Behrens–Fisher problem: Comparison of five test methods.  _Communications in Statistics—Simulation and Computation_ , 37(6):1064–1085, 2008.
  * Commenges (2003) Daniel Commenges.  Transformations which preserve exchangeability and application to permutation tests.  _Nonparametric statistics_ , 15(2):171–185, 2003.
  * Cook and Weisberg (1983) R. Dennis Cook and Sanford Weisberg.  Diagnostics for heteroscedasticity in regression.  _Biometrika_ , 70(1):1–10, 1983.
  * Davidson and Flachaire (2008) Russell Davidson and Emmanuel Flachaire.  The wild bootstrap, tamed at last.  _Journal of Econometrics_ , 146(1):162–169, 2008.
  * Eicker (1967) Friedhelm Eicker.  Reducing tcb complexity for security-sensitive applications: three case studies.  _Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability_ , pages 59–82, 1967.
  * Fisher (1935) Ronald A. Fisher.  The fiducial argument in statistical inference.  _Annals of Eugenics_ , 6:391–398, 1935.
  * Fisher (1941) Ronald A. Fisher.  The asymptotic approach to Behrens’ integral with further tables for the d test of significance.  _Annals of Eugenics_ , 11:141–172, 1941.
  * Freedman (2006) David A. Freedman.  On the so-called “Huber sandwich estimator” and “robust standard errors".  _The American Statistician_ , 60(4):299–302, 2006.
  * Glejser (1969) Herbert Glejser.  A new test for heteroskedasticity.  _Journal of the American Statistical Association_ , 64(325):316–323, 1969.
  * Goldfeld and Quandt (1965) Stephen M. Goldfeld and Richard E. Quandt.  Some tests for homoscedasticity.  _Journal of the American Statistical Association_ , 60(310):539–547, 1965.
  * Hemerik and Goeman (2018) Jesse Hemerik and Jelle J. Goeman.  Exact testing with random permutations.  _TEST_ , 27:811–825, 2018.
  * Hemerik et al. (2020) Jesse Hemerik, Jelle. J. Goeman, and Livio Finos.  Robust testing in generalized linear models by sign flipping score contributions.  _Journal of The Royal Statistical Society Series B-statistical Methodology_ , 82(3):841–864, 2020.
  * Huber (1967) Peter J. Huber.  The behavior of maximum likelihood estimates under nonstandard conditions.  _Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability_ , pages 221–233, 1967.
  * Huber and Ronchetti (2009) Peter J. Huber and Elvezio M. Ronchetti.  _Robust Statistics_.  Wiley, New York, 2009.
  * Huh and Jhun (2001) Myung-Hoe Huh and Myoungshic Jhun.  Random permutation testing in multiple linear regression.  _Communications in Statistics - Theory and Methods_ , 30(10):2023–2032, 2001.
  * Jarque and Bera (1980) Carlos M. Jarque and Anil K. Bera.  Efficient tests for normality, homoscedasticity and serial independence of regression residuals.  _Economics Letters_ , 6(3):255–259, 1980.
  * Kauermann and Carroll (2000) Goeran Kauermann and Raymond J. Carroll.  The sandwich variance estimator: efficiency properties and coverage probability of confidence intervals.  _Discussion Paper 189_ , 2000.
  * Kherad-Pajouh and Renaud (2010) Sara Kherad-Pajouh and Olivier Renaud.  An exact permutation method for testing any effect in balanced and unbalanced fixed effect anova.  _Computational Statistics & Data Analysis_, 54(7):1881 – 1893, 2010.  ISSN 0167-9473.
  * Kim and Cohen (1998) Seock-Ho Kim and Allan S. Cohen.  On the Behrens–Fisher problem: a review.  _Journal of Educational and Behavioral Statistics_ , 23(4):356–377, 1998.
  * Long and Ervin (2000) J. Scott Long and Laurie H. Ervin.  Using heteroscedasticity consistent standard errors in the linear regression model.  _The American Statistician_ , 54(3):217–224, 2000.
  * Love et al. (2014) Michael I. Love, Wolfgang Huber, and Simon Anders.  Moderated estimation of fold change and dispersion for rna-seq data with deseq2.  _Genome biology_ , 15(12):1–21, 2014.
  * Maas and Hox (2004) Cora J. M. Maas and Joop J. Hox.  Robustness issues in multilevel regression analysis.  _Statistica Neerlandica_ , 58(2):127–137, 2004.
  * Magnus and Neudecker (2019) Jan Magnus and Heinz Neudecker.  _Matrix Differential Calculus with Applications in Statistics and Econometrics_.  Wiley, New York, 2019.
  * McCullagh and Nelder (1989) Peter McCullagh and John Nelder.  _Generalized Linear Models_.  Chapman and Hall, Boca Raton, FL, 1989.
  * Neyman and Scott (1948) Jerzy Neyman and Elizabeth L. Scott.  Consistent estimates based on partially consistent observations.  _Econometrica_ , 16:1–32, 1948.
  * Pesarin (2001) Fortunato Pesarin.  _Multivariate permutation tests: with applications in biostatistics._ Wiley, Chichester, 2001.
  * Rochon et al. (2012) Justine Rochon, Matthias Gondan, and Meinhard Kieser.  To test or not to test: Preliminary assessment of normality when comparing two independent samples.  _BMC Medical Research Methodology_ , 12(81), 2012.
  * Solari et al. (2014) Aldo Solari, Livio Finos, and Jelle J. Goeman.  Rotation-based multiple testing in the multivariate linear model.  _Biometrics_ , 70(4):954–961, 2014.
  * Van der Vaart (1998) Aad W. Van der Vaart.  _Asymptotic Statistics_.  Cambridge University Press, Cambridge, 1998.
  * White (1982) Halbert White.  Maximum likelihood estimation of misspecified models.  _Econometrica_ , 50(1):1–25, 1982.
  * Winkler et al. (2014) Anderson M. Winkler, Gerard R. Ridgway, Matthew A. Webster, Stephen M. Smith, and Thomas E. Nichols.  Permutation inference for the general linear model.  _Neuroimage_ , 92:381–397, 2014.

##  Appendix A Appendix: Proofs of the Theorems

We first introduce some Lemmas needed for the proofs of the Theorems.

The projection matrix for GLMs introduced in (2) is a proper projection matrix
for studentized units as shown in the following Lemma, which is mentioned for
instance in (Agresti, 2015, p. 136). We give a full proof for the sake of
completeness, since that is omitted in most textbooks.

######  Lemma 1.

The fitted and observed values in a GLM are connected through the relation

$V^{-1/2}(\hat{\mu}-\mu)=HV^{-1/2}(y-\mu)\\{1+o_{p}(1)\\}$

###### Proof.

The first-order approximation of the score function is

$\displaystyle s_{\hat{\beta}}$
$\displaystyle=s_{\beta}+\mathcal{I}_{\beta,\beta}(\hat{\beta}-\beta)+o_{p}(1)$
$\displaystyle 0$
$\displaystyle=X^{T}D{V}^{-1}(y-\mu)-X^{T}{W}X(\hat{\beta}-\beta)+o_{p}(1)$
$\displaystyle\hat{\beta}-\beta$
$\displaystyle=(X^{T}{W}X)^{-1}X^{T}D{V}^{-1}(y-\mu)+o_{p}(1)$

where $(o_{p}(1))$ is an error term asymptotically negligible. Then, by the
Delta method we have

$\displaystyle\hat{\mu}$
$\displaystyle=\mu+\frac{d\mu}{d\eta}\frac{d\eta}{d\beta}(\hat{\beta}-\beta)\\{1+o_{p}(1)\\}$
$\displaystyle\hat{\mu}-\mu$
$\displaystyle=DX^{T}(\hat{\beta}-\beta)\\{1+o_{p}(1)\\}$
$\displaystyle=DX^{T}(X^{T}{W}X)^{-1}X^{T}D{V}^{-1}(y-\mu)\\{1+o_{p}(1)\\}$
$\displaystyle=D{W}^{-1/2}{H}{W}^{1/2}D^{-1}(y-\mu)\\{1+o_{p}(1)\\}$
$\displaystyle={V}^{1/2}{H}{V}^{-1/2}(y-\mu)\\{1+o_{p}(1)\\}$
$\displaystyle{V}^{-1/2}(\hat{\mu}-\mu)$
$\displaystyle={H}{V}^{-1/2}(y-\mu)\\{1+o_{p}(1)\\}.$

where the asymptotically negligible error term has two sources, one related to
the second-order approximation of the likelihood and the other related to
possible non-linearity of the link function. ∎

######  Lemma 2.

Let $C$ be any $n$-dimensional matrix and $G$ be a nonsingular $n\times n$
matrix, then $C$ and $G^{-1}CG$ have the same set of eigenvalues (with the
same multiplicities).

###### Proof.

See (Magnus and Neudecker, 2019, p. 15). ∎

######  Lemma 3.

Let $C$ be any $n$-dimensional matrix and $\mathbb{F}$ be the set of all
$n$-dimensional flipping matrices, i.e. the set of all possible
$n$-dimensional diagonal matrices $F$ with elements $-1$ or $1$. Then

$\sum_{F\in\mathbb{F}}FCF=2^{n}\mbox{diag}(C),$

where $\mbox{diag}(C)$ is a diagonal matrix with the same diagonal elements of
$C$.

###### Proof.

First we note that the absolute value of each element of $C$ does not change
after the multiplication $FCF$. The sums for the off-diagonal elements contain
an equal number of terms with positive and negative sign, hence their sum over
all possible flips is zero, while the sign of each diagonal element is
positive for each term. By noting that the total number of flips is $2^{n}$ we
have the claim. ∎

The following result was adapted from Huber and Ronchetti (2009), extending
their theorem (Proposition 7.1, p. 156) for linear regression models to GLMs.

######  Lemma 4.

Assume that the regression coefficients of a generalized linear model are
consistently estimated, in the sense that for every $\varepsilon>0$, as
$n\rightarrow\infty$,

$\max_{1\leq i\leq n}\mbox{pr}(|\hat{\mu}_{i}-\mu_{i}|>\varepsilon)\rightarrow
0.$

Then

$\max_{1\leq i\leq n}h_{ii}\xrightarrow{}0,$

where $h_{ik}$ is the $ik$-th element of the matrix $H$, as defined in (2).

###### Proof.

Using Lemma 1 we have, for each $i$, $\hat{\mu}_{i}-\mu_{i}=$

$h_{ii}(y_{i}-\mu_{i})+\sum_{k\neq
i}v_{i}^{1/2}v_{k}^{-1/2}h_{ik}(y_{k}-\mu_{k})+o_{p}(1).$ (8)

Now we give a general probability result. Let $V_{1}$ and $V_{2}$ be two
independent random variables, for any $\varepsilon>0$ we have that

$\displaystyle\mbox{pr}(|V_{1}+V_{2}|$
$\displaystyle\geq\varepsilon)\geq\mbox{pr}(V_{1}\geq\varepsilon)\mbox{pr}(V_{2}\geq
0)+\mbox{pr}(V_{1}\leq-\varepsilon)\mbox{pr}(V_{2}\leq 0)$
$\displaystyle\geq\min\left\\{\mbox{pr}(V_{1}\geq\varepsilon),\mbox{pr}(V_{1}\leq-\varepsilon)\right\\}.$

Noting that $(y_{i}-\mu_{i})$ is independent from $(y_{k}-\mu_{k})$ for each
$k\neq i$ we can apply the result to expression (8) to obtain

$\mbox{pr}\left(|\hat{\mu}_{i}-\mu_{i}\geq\varepsilon|\right)\geq\min\left[\mbox{pr}\left\\{(y_{i}-\mu_{i})\geq\frac{\varepsilon}{h_{ii}}\right\\},\mbox{pr}\left\\{(y_{i}-\mu_{i})\leq-\frac{\varepsilon}{h_{ii}}\right\\}\right]=m_{i}^{n}.$

We have $\max_{1\leq i\leq
n}\mbox{pr}(|\hat{\mu}_{i}-\mu_{i}|>\varepsilon)\rightarrow 0$, so
$\max_{1\leq i\leq n}m_{i}^{n}\xrightarrow{}0$. Since $\varepsilon$ was
arbitrary, this implies that $\max_{1\leq i\leq n}h_{ii}\xrightarrow{}0$. ∎

######  Lemma 5.

The computational cost of the standardization constant in (5) is linear in
$n$.

###### Proof.

Define $W^{1/2}Z=U\Delta L^{T}$, that is, the singular value decomposition of
$W^{1/2}Z$, where $U$ is a semiorthogonal $n\times q$ matrix ($q$ equal to the
rank of $Z$), $\Delta$ a diagonal $q$ matrix and $L$ a $q\times q$ orthogonal
matrix. Therefore (2) can be written as

$H=W^{1/2}Z(Z^{T}WZ)^{-1}Z^{T}W^{1/2}=U\Delta L^{T}L\Delta U^{T}=UU^{T}.$

Now, let $a=(I-H)W^{1/2}X$ and $A=diag(a)$. Further, let $\mathbf{1}$ be an
$n$-dimensional vector of ones. The denominator of the standardized test
statistic becomes

$\displaystyle X^{T}W^{1/2}(I-H)F(I-H)F(I-H)W^{1/2}X=$ $\displaystyle=$
$\displaystyle a^{T}F(I-UU^{T})Fa=$ $\displaystyle=$ $\displaystyle
a^{T}FIFa-a^{T}FUU^{T}Fa$ $\displaystyle=$ $\displaystyle
a^{T}a-\mathbf{1}^{T}AFUU^{T}FA\mathbf{1}$ $\displaystyle=$ $\displaystyle
a^{T}a-\mathbf{1}^{T}FAUU^{T}AF\mathbf{1}$ $\displaystyle=$ $\displaystyle
a^{T}a-f^{T}CC^{T}f,$

where $f=F\mathbf{1}$ and $C=AU$.

Therefore, given that $a^{T}a$ is a constant and the computational cost of
$f^{T}CC^{T}f$ is linear with $n$ since we can write
$f^{T}CC^{T}f=\sum_{j=1}^{q}(\sum_{i=1}^{n}f_{i}C_{ij})^{2}$, the result of
the lemma follows. ∎

###  A.1 Proof of Theorem 1

######  Theorem.

Assume that the variances are correctly specified, that is, $\tilde{V}=V$, and
that Assumption 1-3 hold. For $n\xrightarrow{}\infty$, the test that rejects
$H_{0}$ if (4) holds is an asymptotically $\alpha$ level test.

###### Proof.

By the definition of the random sign-flipping transformations is trivial to
observe that the expected value of the test statistic $S(\mathcal{F})$ is
zero.

By Lemma 1 we can rewrite the effective score statistic (3) as

$S(\mathcal{F})=n^{-1/2}X^{T}W^{1/2}(I-H)\mathcal{F}(I-H)V^{-1/2}(y-\mu)+o_{p}(1).$

Let

$a=(I-H)W^{1/2}X.$ (9)

The variance conditional on $\mathcal{F}=F$ is

$\displaystyle\mbox{var}\left\\{S(F)\right\\}$
$\displaystyle=n^{-1}a^{T}F(I-H)V^{-1/2}E\\{(y-\mu)(y-\mu)^{T}\\}V^{-1/2}(I-H)Fa+o_{p}(1)$
$\displaystyle=n^{-1}a^{T}F(I-H)V^{-1/2}VV^{-1/2}(I-H)Fa+o_{p}(1)$
$\displaystyle=n^{-1}a^{T}F(I-H)Fa+o_{p}(1)$

and for $\mathcal{F}=\mathbf{I}$ we have

$\mbox{var}\left\\{S(\mathbf{I})\right\\}=n^{-1}a^{T}Ia+o_{p}(1).$

Taking the difference

$\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(F)\right\\}=n^{-1}a^{T}\left\\{I-F(I-H)F\right\\}a+o_{p}(1),$

we note that the first term is a quadratic form and we look at the matrix

$I-F(I-H)F=I-FF+FHF=FHF.$

Since $H$ is a projection matrix and $F^{-1}=F$, Lemma 2 implies that $FHF$ is
positive semidefinite, so that asymptotically
$\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(F)\right\\}\geq
0$.

Define

$h_{\sup}=\sup_{1\leq i\leq n}\left(h_{ii}\right).$

Let us make the randomness of the flips explicit. Let $\mathbb{F}$ denote the
set of all possible flipping matrices, and note that $|\mathbb{F}|=2^{n}$. By
Lemma 3 we have

$\displaystyle
E\left[\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(\mathcal{F})\mid
F\right\\}\right]$ $\displaystyle=$
$\displaystyle(2^{n})^{-1}\sum_{F\in\mathbb{F}}[\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(F)\right\\}]$
$\displaystyle=$ $\displaystyle n^{-1}a^{T}\mathrm{diag}(H)a+o_{p}(1)$
$\displaystyle\leq$ $\displaystyle n^{-1}\|a\|^{2}\cdot h_{\sup}+o_{p}(1),$

where, by Assumption 2 and Lemma 4, the limiting behavior is

$\lim_{n\xrightarrow{}\infty}n^{-1}\|a\|^{2}\cdot h_{\sup}=0.$

According to the law of total variance we have

$\mbox{var}\\{S(\mathcal{F})\\}=\mbox{var}\left[E\left\\{S(\mathcal{F})\mid
F\right\\}\right]+E\left[\mbox{var}\left\\{S(\mathcal{F})\mid
F\right\\}\right].$

We know that $E\\{S(\mathcal{F})\mid F\\}$ does not depend on $\mathcal{F}$,
which means that $\mbox{var}[E\\{S(\mathcal{F})\mid F\\}]=0$, so
$\mbox{var}\\{S(\mathcal{F})\\}=E[\mbox{var}\\{S(\mathcal{F})\mid F\\}].$ It
follows that, marginally over $F$,

$\lim_{n\xrightarrow{}\infty}\mbox{var}\\{S(I)\\}-\mbox{var}\\{S(\mathcal{F})\\}=0.$

Since the random sign-flipping transformations are all independent, the
corresponding test statistics are all uncorrelated, i.e, for all $1\leq
l<m\leq g$ we have $\mbox{cov}\\{S(\mathcal{F}_{l}),S(\mathcal{F}_{m})\\}=0$.
By the previous results and Assumption 3,
$(S(\mathbf{I}),\dots,S(\mathcal{F}_{g}))^{T}$ converges to a multivariate
normal distribution by the multivariate Lindberg-Feller central limit theorem
(Van der Vaart, 1998), with mean vector 0 and covariance matrix
$s^{2}\mathbf{I},$ where $s^{2}$ is the limiting variance of the flipped test
statistic. Finally, we use Lemma 1 of Hemerik et al. (2020) to conclude that
the test that rejects when (4) holds is an asymptotic $\alpha$ level test. ∎

###  A.2 Proof of Proposition 1

######  Proposition.

Consider a normal regression model with identity link. Assume that the
variances are correctly specified, that is, $\tilde{V}=V$, and that Assumption
1-3 hold. For finite sample size, the effective score statistic defined as in
(3) has $\mbox{var}\\{S(\mathbf{I})\\}>\mbox{var}\\{S(\mathcal{F})\\}$.

###### Proof.

By Lemma 1 we can rewrite the effective score statistic (3) as

$S(\mathcal{F})=n^{-1/2}X^{T}W^{1/2}(I-H)\mathcal{F}(I-H)V^{-1/2}(y-\mu)+o_{p}(1).$

The variance conditional on $\mathcal{F}=F$ is

$\displaystyle\mbox{var}\left\\{S(F)\right\\}$
$\displaystyle=n^{-1}a^{T}F(I-H)V^{-1/2}E\\{(y-\mu)(y-\mu)^{T}\\}V^{-1/2}(I-H)Fa$
$\displaystyle=n^{-1}a^{T}F(I-H)V^{-1/2}VV^{-1/2}(I-H)Fa$
$\displaystyle=n^{-1}a^{T}F(I-H)Fa$

and for $\mathcal{F}=\mathbf{I}$ we have

$\mbox{var}\left\\{S(\mathbf{I})\right\\}=n^{-1}a^{T}Ia.$

Taking the difference

$\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(F)\right\\}=n^{-1}a^{T}\left\\{I-F(I-H)F\right\\}a,$

we note that it is a quadratic form and we consider the matrix

$I-F(I-H)F=I-FF+FHF=FHF.$

Since $H$ is a projection matrix and $F^{-1}=F$, Lemma 2 implies that $FHF$ is
positive semidefinite, and therefore
$\mbox{var}\\{S(\mathbf{I})\\}-\mbox{var}\\{S(F)\\}\geq 0$.

We then prove the strict inequality for at least one flipping matrix $F$.
Taking any model with an intercept, by construction

$h_{\inf}=\inf_{1\leq i\leq n}\left(h_{ii}\right)>0.$

Note that $|\mathbb{F}|=2^{n}$. Since
$\mbox{var}\\{S(\mathbf{I})\\}-\mbox{var}\\{S(F)\\}\geq 0$ as proven above, it
suffices to show that

$\sum_{F\in\mathbb{F}}n^{-1}a^{T}FHFa>0,$

which means that we have a strictly inequality for some $F$. Using Lemma 3 we
have

$\sum_{F\in\mathbb{F}}FHF=2^{n}\,\mathrm{diag}(H).$

Therefore

$\displaystyle\sum_{F\in\mathbb{F}}n^{-1}a^{T}FHFa$ $\displaystyle=$
$\displaystyle 2^{n}n^{-1}a^{T}\textrm{diag}(H)a$ $\displaystyle\geq$
$\displaystyle 2^{n}n^{-1}\|a\|^{2}\cdot h_{\inf}>0.$

∎

###  A.3 Proof of Proposition 2

######  Proposition.

Assume that the variances are correctly specified, that is, $\tilde{V}=V$, and
that Assumptions 1-3 hold. The standardized sign-flip score statistic is
finite sample second-moment null-invariant. The test that rejects $H_{0}$ if
(6) holds is i) an exact $\alpha$ level test in case of normal response with
identity link and ii) asymptotically exact in all other cases.

###### Proof.

The proof is immediate. We simply observe that the expected value of the new
statistic (5) is left unchanged, while the standardization makes the variance
of each flipped new statistic equal to 1. The same argument given in the last
part of the proof of Theorem 1 applies to deduce the asymptotic exactness of
the test.

In the special case of the normal model with the identity link, we have that
$\mbox{var}(S(F))$ does not depend on any unknown nuisance parameters.
Moreover, in this model $S(\mathbf{I})$ and $S(\mathcal{F})$ are normally
distributed in finite samples, so the standardized sign-flip score statistic
achieves not just second-moment but full null-invariance, resulting in an
exact test by Theorem 2 of Hemerik and Goeman (2018). ∎

###  A.4 Proof of Proposition 3

######  Proposition.

Assume that Assumptions 1-3 hold. If the variances are misspecified by any
finite constant $c>0$, that is, $V=c\tilde{V}$, the standardized sign-flip
score statistic is second-moment null-invariant. The test that rejects $H_{0}$
if (6) holds is i) an exact $\alpha$ level test in case of normal response
with identity link and ii) asymptotically exact in all other cases.

###### Proof.

We observe that

$V\tilde{V}^{-1}=\phi\tilde{\phi}^{-1}I=cI$

and therefore

$\mbox{var}\\{S^{*}(F)\\}=c\quad\forall F\in\mathbb{F}.$

The result than follows from the proof of Proposition 2. ∎

###  A.5 Proof of Theorem 2

######  Theorem.

Assume that the variances are misspecified, that is, $V\neq\tilde{V}$ and that
Assumptions 1-3 hold. For $n\xrightarrow{}\infty$, the standardized sign-flip
score statistic is asymptotically second-moment null-invariant. The test that
rejects $H_{0}$ if (6) holds is an asymptotically $\alpha$ level test.

###### Proof.

It is trivial to see that the expected value of the test statistic is not
affected by this misspecification.

Let

$B=V\tilde{V}^{-1},$

note that it is a diagonal matrix and all its elements are finite and greater
than zero by Assumption 1. We compute again the variance of (3), but now we
consider the misspecification.

Let $\tilde{H}$ and $\tilde{a}$ be the quantities defined in (2) and (9) for
$W=\tilde{W}$. The variance can be written as

$\displaystyle\mbox{var}\left\\{S(F)\right\\}$
$\displaystyle=n^{-1}\tilde{a}^{T}F(I-\tilde{H})\tilde{V}^{-1/2}E\left\\{(y-\mu)(y-\mu)^{T}\right\\}\tilde{V}^{-1/2}(I-\tilde{H})F\tilde{a}+o_{p}(1)$
$\displaystyle=n^{-1}\tilde{a}^{T}F(I-\tilde{H})B(I-\tilde{H})F\tilde{a}+o_{p}(1).$

Take the difference

$\displaystyle\mbox{var}\left(S(\mathbf{I})\right)-\mbox{var}\left(S(F)\right)$
$\displaystyle=n^{-1}\tilde{a}^{T}\left\\{B-F(I-\tilde{H})B(I-\tilde{H})F\right\\}\tilde{a}+o_{p}(1)$
$\displaystyle=n^{-1}\tilde{a}^{T}\left[F\\{\tilde{H}B+(I-\tilde{H})B\tilde{H}\\}F\right]\tilde{a}+o_{p}(1).$

We notice that by Lemma 3

$\sum_{F\in\mathcal{F}}F(\tilde{H}B+B\tilde{H}-\tilde{H}B\tilde{H})F=2^{n}\,\mathrm{diag}(\tilde{H}B+B\tilde{H}-\tilde{H}B\tilde{H})$

where the $i$-th element of that diagonal matrix is

$2\tilde{h}_{ii}b_{i}-\sum_{k=1}^{n}\tilde{h}^{2}_{ik}b_{k}.$

Then we have

$\displaystyle
E\left[\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(\mathcal{F})\mid
F\right\\}\right]$
$\displaystyle=(2^{n})^{-1}\sum_{F\in\mathbb{F}}\left[\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(F)\right\\}\right]$
$\displaystyle=n^{-1}\tilde{a}^{T}\mathrm{diag}(\tilde{H}B+B\tilde{H}-\tilde{H}B\tilde{H})\tilde{a}+o_{p}(1)$
$\displaystyle\leq
n^{-1}\tilde{a}^{T}\mathrm{diag}(\tilde{H}B+B\tilde{H})\tilde{a}+o_{p}(1).$

By Assumption 1, note that there exists two finite positive constants
$c_{1},c_{2}$ such that for each $i$-th element

$2\tilde{h}_{ii}b_{i}\leq c_{1}b_{sup}\tilde{h}_{ii}\leq c_{2}\sup_{1\leq
i\leq n}\tilde{h}_{ii}=c_{2}\tilde{h}_{sup}$

where

$\tilde{h}_{\sup}=\sup_{1\leq i\leq n}(\tilde{h}_{ii})\quad
b_{\sup}=\sup_{1\leq i\leq n}(b_{i}).$

Therefore we can derive the upper bound

$E\left[\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(\mathcal{F})\mid
F\right\\}\right]\leq n^{-1}\|\tilde{a}\|^{2}c_{2}\cdot\tilde{h}_{\sup}.$

Using Assumption 2 and Lemma 4, the limiting behavior is

$\lim_{n\xrightarrow{}\infty}n^{-1}\|\tilde{a}\|^{2}c_{2}\cdot\tilde{h}_{\sup}=0.$

Meanwhile, for two positive constants $c_{3},c_{4}$ we have for each $i$-th
element

$-\sum_{k=1}^{n}\tilde{h}^{2}_{ik}b_{k}\geq-\sum_{k=1}^{n}\tilde{h}^{2}_{ik}c_{3}b_{\sup}=-c_{3}b_{\sup}\tilde{h}_{ii}\geq-
c_{4}\tilde{h}_{\sup}.$

Then using again Lemma 3 we can derive the lower bound

$\displaystyle
E\left[\mbox{var}\left\\{S(\mathbf{I})\right\\}-\mbox{var}\left\\{S(\mathcal{F})\mid
F\right\\}\right]$ $\displaystyle\geq
n^{-1}\tilde{a}^{T}\mathrm{diag}(-\tilde{H}B\tilde{H})\tilde{a}+o_{p}(1)$
$\displaystyle\geq-n^{-1}\|\tilde{a}\|^{2}c_{4}\cdot\tilde{h}_{\sup}+o_{p}(1)\,$

where, using Assumption 2 and Lemma 4, the limiting behavior is

$\lim_{n\xrightarrow{}\infty}-n^{-1}\|\tilde{a}\|^{2}c_{4}\cdot\tilde{h}_{\sup}=0.$

According to the law of total variance we have

$\mbox{var}\\{S(\mathcal{F})\\}=\mbox{var}\left[E\\{S(\mathcal{F})\mid
F\\}\right]+E\left[\mbox{var}\\{S(\mathcal{F})\mid F\\}\right].$

We know that $E\\{S(\mathcal{F})\mid F\\}$ does not depend on $\mathcal{F}$,
which means that $\mbox{var}[E\\{S(\mathcal{F})\mid F\\}]=0$, so
$\mbox{var}\\{S(\mathcal{F})\\}=E[\mbox{var}\\{S(\mathcal{F})\mid F\\}].$ It
follows that, marginally over $F$,

$\lim_{n\xrightarrow{}\infty}\mbox{var}\left\\{S(I)\right\\}-\mbox{var}\left\\{S(\mathcal{F})\right\\}=0.$

The same argument given in the last part of the proof of theorem 1 applies to
deduce the asymptotic exactness of the test considered. ∎

###  A.6 Proof of Proposition 4

######  Proposition.

Assume that the variances are correctly specified, that is, $\tilde{V}=V$, and
that Assumptions 1, 4 and 5 hold. The $d$-dimensional standardized sign-flip
score vector is finite sample second-moment null-invariant. The test that
rejects $H_{0}$ if (7) holds is an asymptotically $\alpha$ level test.

###### Proof.

The proof is analogous to Theorem 3 of Hemerik et al. (2020). The asymptotic
exactness in all the cases is a consequence of the use of the continuous
mapping theorem. ∎

###  A.7 Proof of Theorem 3

######  Theorem.

Assume that the variances are misspecified, that is, $V\neq\tilde{V}$ and that
Assumptions 1,4,5 hold. For $n\xrightarrow{}\infty$, the $d$-dimensional
standardized sign-flip score vector is asymptotically second-moment null-
invariant. The test that rejects $H_{0}$ if (7) holds is an asymptotically
$\alpha$ level test.

###### Proof.

It is trivial to observe that the expected value is a $d$-dimensional zero
vector. We focus on the variance, which now is a $d\times d$ matrix. We prove
the robustness elementwise.

The proof for the diagonal elements follows from Theorem 2. Then, we focus on
one covariance term of the covariance matrix of the flipped effective score
vector. Let $X=(X_{1},\dots,X_{d})^{T}$ and

$\displaystyle\tilde{a}_{1}$ $\displaystyle=(I-\tilde{H})\tilde{W}^{1/2}X_{1}$
$\displaystyle\tilde{a}_{2}$
$\displaystyle=(I-\tilde{H})\tilde{W}^{1/2}X_{2}.$

When the variances are misspecified

$\displaystyle\mbox{cov}\\{S^{1}(F),S^{2}(F)\\}=$
$\displaystyle\;=n^{-1}\tilde{a}_{1}^{T}F(I-\tilde{H})\tilde{V}^{-1/2}E\left\\{(y-\mu)(y-\mu)^{T}\right\\}\tilde{V}^{-1/2}(I-\tilde{H})F\tilde{a}_{2}+o_{p}(1)$
$\displaystyle\;=n^{-1}\tilde{a}_{1}^{T}F(I-\tilde{H})B(I-\tilde{H})F\tilde{a}_{2}+o_{p}(1).$

Take the difference

$\displaystyle\mbox{cov}\\{S^{1}(I),S^{2}(I)\\}-\mbox{cov}\\{S^{1}(F),S^{2}(F)\\}=$
$\displaystyle\;=n^{-1}\tilde{a}_{1}^{T}\left\\{B-F(I-\tilde{H})B(I-\tilde{H})F\right\\}\tilde{a}_{2}+o_{p}(1)$
$\displaystyle\;=n^{-1}\tilde{a}_{1}^{T}\left\\{F(\tilde{H}B+B\tilde{H}-\tilde{H}B\tilde{H})F\right\\}\tilde{a}_{2}+o_{p}(1).$

From this point the asymptotic second-moment null-invariance can be derived
directly from the proof of the Theorem 2, by applying the same reasoning to
each covariance term, replacing Assumptions 2-3 with Assumptions 4-5. The
asymptotic exactness of the test follows from the proof of Proposition 4. ∎


# Relational Data Selection for Data Augmentation of Speaker-dependent Multi-
band MelGAN Vocoder

###### Abstract

Nowadays, neural vocoders can generate very high-fidelity speech when a bunch
of training data is available. Although a speaker-dependent (SD) vocoder
usually outperforms a speaker-independent (SI) vocoder, it is impractical to
collect a large amount of data of a specific target speaker for most real-
world applications. To tackle the problem of limited target data, a data
augmentation method based on speaker representation and similarity measurement
of speaker verification is proposed in this paper. The proposed method selects
utterances that have similar speaker identity to the target speaker from an
external corpus, and then combines the selected utterances with the limited
target data for SD vocoder adaptation. The evaluation results show that,
compared with the vocoder adapted using only limited target data, the vocoder
adapted using augmented data improves both the quality and similarity of
synthesized speech.

Index Terms: neural vocoder, speaker similarity, data augmentation, multi-band
MelGAN, x-vector

##  1 Introduction

A vocoder [1] is a speech codec to analyze speech into acoustic features and
synthesize the acoustic features back to speech. Conventional vocoders [2, 3]
built following the source-filter model [4] have been widely used in speech
synthesis tasks. However, the quality of the synthetic speech is usually
degraded because the phase information and temporal details are discarded
during the analysis-synthesis process of the conventional vocoders. Many
neural speech generation models [5, 6, 7, 8, 9, 10, 11, 12, 13, 14] have been
proposed to directly model speech waveforms without many ad hoc assumptions of
speech generation. Using these models as vocoders [15, 16] to synthesize
speech based on the acoustic features extracted by conventional vocoders also
remarkably improve the naturalness of the synthetic speech.

Although a speaker-dependent (SD) model usually outperforms a speaker-
independent (SI) model [15, 16], collecting much data from a user is
impractical. An efficient way to develop an SD vocoder is to first train a
multi-speaker vocoder using a multi-speaker corpus, and then adapt the multi-
speaker vocoder to the SD vocoder using the few-shot target data. However,
this method still requires about five minutes of target data to develop a
stable SD vocoder [17, 18]. On the other hand, an SI vocoder trained with a
varied corpus may outperform an SD vocoder trained with a relatively small
corpus in terms of speech quality [19]. However, the difference of speaker
similarity has not been well investigated. In this paper, we first explore a
more challenging scenario, i.e., one-shot adaptation, where the available
target data is around 30s. Then, the speaker similarity difference between SI
and SD vocoders is investigated.

The performance of most neural models is highly correlated with the amount and
diversity of training data due to the data-driven nature. Therefore, the use
of generative models to generate augmented data is straightforward and has
been applied to many recognition tasks. For example, well-trained text-to-
speech (TTS) systems have been used to generate augmented data for training
automatic speech recognition (ASR) [20, 21, 22] and speech translation [23]
systems. A vocoder has been used to generate augmented data with a variety of
fundamental frequency ($F_{0}$) patterns to train an $F_{0}$ estimator [24],
and a voice conversion (VC) framework has been used to generate unseen speaker
data to train a speaker embedding extractor [25]. Even generation tasks, such
as TTS [26], can benefit from using augmented data generated by another TTS
system. However, there is still a performance gap between the models trained
with sufficient natural data and augmented synthetic data [21].

In this paper, different from using generative models to produce augmented
data, since leveraging natural data may avoid error propagation and keep the
entire framework simple, we propose a data selection method to select
augmented utterances from an external corpus. The speaker representation and
similarity measurement from speaker verification (SV) [27] are used to
formulate the selection criteria, and the selected utterances are presumed to
have similar speaker identities to the target speaker. An SI vocoder is first
trained with the multi-speaker external corpus, and then the one-shot target
data and selected augmented utterances are used together to adapt the SI
vocoder to the SD vocoder. Multi-band MelGAN [28] is adopted as the neural
vocoder because of its light architecture. According to our evaluation
results, the vocoder adapted with one-shot target data and augmented data
achieves higher quality and similarity of synthesized speech compared with the
vocoder adapted with only one-shot target data. To our best knowledge, our
method is the first approach to train SD neural vocoders by using SV
technology for data augmentation from an external corpus.

##  2 Baseline multi-band MelGAN vocoder

In this section, we introduce the baseline multi-band MelGAN vocoder, which is
a convolutional neural network (CNN)-based non-autoregressive (non-AR) raw
waveform generation model.

###  2.1 MelGAN with a multi-resolution STFT loss

A classic generative adversarial net (GAN) [29] architecture is adopted in
MelGAN to convert the input mel-spectrogram into speech samples. Instead of
one discriminator, MelGAN utilizes $K$ discriminators ($D_{k}$) running at
different temporal resolutions to capture the hierarchical structure of
speech. Given a generator ($G$), natural speech $\bm{y}$, and the input mel-
spectrogram $\bm{s}$, the discriminator adversarial loss ($L_{\mathrm{D}}$) is
formulated as

$\displaystyle L_{\mathrm{D}}(G,D)$
$\displaystyle=\displaystyle\sum\limits_{k=1}^{K}\frac{1}{K}(\mathbb{E}_{\bm{y}}\left[(1-D_{k}(\bm{y}))^{2}\right]+\mathbb{E}_{\bm{s}}\left[D_{k}(G(\bm{s}))^{2}\right]),$
(1)

where $k$ is the discriminator index.

Moreover, the generator loss ($L_{\mathrm{G}}$) includes an auxiliary loss in
addition to the original adversarial loss ($L_{\mathrm{adv}}$). Specifically,
to improve the stability of training, vanilla MelGAN adopts a feature matching
loss to regularize the discriminator intermediate feature maps of real and
fake data. However, since the multi-resolution short-time Fourier transform
(STFT) loss [12] is more meaningful and can make the training process converge
fast [28], the feature matching loss is replaced by the multi-resolution STFT
loss ($L_{\mathrm{sp}}$) when training our MelGAN vocoders. The generator loss
is formulated as

$\displaystyle
L_{\mathrm{G}}(G,D)=\lambda_{\mathrm{adv}}L_{\mathrm{adv}}(G,D)+L_{\mathrm{sp}}(G),$
(2)

where $\lambda_{\mathrm{adv}}$ is a balance weight, which is set to 2.5 in
this work. The $L_{\mathrm{adv}}$ loss is formulated as

$\displaystyle
L_{\mathrm{adv}}(G,D)=\mathbb{E}_{\bm{s}}\left[(1-D(G(\bm{s})))^{2}\right].$
(3)

The $L_{\mathrm{sp}}$ loss is calculated based on the STFT features extracted
using three setting groups, including the FFT size of [1024, 2048, 4096], the
frame shift of [120, 240, 480], and the window length of [600, 1200, 2400].

The generator and discriminators of MelGAN are fully convolutional networks.
The generator adopts several transposed CNN layers with residual networks and
dilated CNN layers [30] to gradually upsample the input mel-spectrogram to
match the temporal resolution of the output waveform samples. A LeakyReLU [31]
activation is adopted following each CNN layer except the last output layer,
which uses a tanh function to output waveforms. The multi-scale discriminators
have an identical network structure but different downsampling factors.
Downsampling is implemented using stride average pooling. More details can be
found in the open source repository111https://github.com/kan-
bayashi/ParallelWaveGAN.

###  2.2 Multi-band approach

Directly modeling speech samples with a high sampling frequency ($f_{s}$) is
challenging because of the speech signal has a high temporal resolution with a
very long-term dependency, which usually result in the consumption of time and
computing resources in the generation process. Decomposing the speech signal
into several subbands can significantly improve the generation efficiency,
because each subband signal is generated in parallel using a single network.
The multi-band approach has been successfully applied to many AR [32, 33, 34]
and non-AR [28] neural vocoders.

The $f_{s}$ of the speech signal processed in this paper is 44.1 KHz. The
analysis and synthesis filters in [34, 28] are used to decompose the speech
signal into five frequency bands, and the fullband signal is restored on the
basis of the subband signals. The generator is trained to generate the subband
signals in parallel. The inputs of the discriminators are the restored
fullband signal. To improve the stability of multi-band training, the multi-
resolution STFT loss is adopted for both fullband and subband signals. The
setting groups of subband STFT include the FFT size of [384, 683, 171], the
frame shift of [30, 60, 10], and the window length of [150, 300, 60]. More
details can be found in [34, 28] and the open source repository††footnotemark:
.

##  3 Relational data selection

To effectively develop the SD vocoder when very limited (one-shot) target data
is available, we propose a data augmentation framework leveraging an external
corpus (candidate pool) for adapting the SI vocoder to the target SD vocoder.
A hierarchical framework based on data relationships is used to select
suitable data for speaker adaptation. Three levels of relationships are
considered. First, the speaker similarity is measured by the inter-speaker
relationship between the target speaker and the candidate speaker. Second, the
inter-candidate-speaker relationships are established to verify the speaker-
wise confidence. Third, the reliability of each candidate utterance is
regularized by the relationship within the speaker.

###  3.1 Speaker similarity

Selecting external utterances with similar speaker identities to the target
speaker for speaker adaptation is straightforward. Identity modeling and
similarity measurement in SV technology can be used in this work. For identity
modeling, we use the state-of-the-art x-vector [27] speaker representation.
For similarity measurement, we use the probabilistic linear discriminant
analysis (PLDA) [35, 36]. X-vector is a speaker embedding extracted from the
intermediate layer of a speaker-identification neural network, and PLDA is a
distance measurement designed to be generalized for arbitrary unseen classes.
Therefore, the first selection criterion is formulated as

$\displaystyle\mathrm{PLDA}(\bm{x}_{n,i},\bm{x}_{Target}),$ (4)

where $\bm{x}$ denotes x-vector, $n$ is the candidate speaker index, $i$ is
the utterance index, and $\bm{x}_{Target}$ is the average x-vector of all
target utterances. The higher the PLDA score, the higher the speaker
similarity.

###  3.2 Speaker-wise divergence regularization

Since a robust speaker embedding should be almost independent of the phonetic
content, we assume that the utterance-wise x-vectors from the same speaker
should be similar. According to the assumption, if the distribution of the
x-vectors of a speaker is diverse, the reliability of these x-vectors may be
low. To model the speaker-wise confidence within the candidate pool, an SD
term is introduced to regularize the PLDA score. First, temperature sigmoid is
applied to make all PLDA scores have the same sign,

$\displaystyle\mathrm{PLDA}^{\prime}(\cdot)=\frac{1}{1+0.5\times
e^{-\mathrm{PLDA}(\cdot)}}.$ (5)

Then, the second selection criterion is formulated as

$\displaystyle\frac{\mathrm{PLDA}^{\prime}(\bm{x}_{n,i},\bm{x}_{Target})}{(\sigma_{n})^{\alpha}},$
(6)

where $\sigma_{n}$ is the square root of the average squared Euclidean
distance between each utterance x-vector and the mean x-vector of speaker $n$,
and the weight $\alpha$ is set to 0.1 in this paper. High x-vector diversity
of a speaker results in low speaker confidence.

###  3.3 Utterance-wise divergence regularization

Following the above assumption, the internal speaker relationship is
introduced to tackle the outlier utterances within each speaker. That is, if
the x-vector of an utterance is very different from the x-vectors of other
utterances of the same speaker, the utterance is considered an outlier, and
its x-vector is unreliable. Therefore, the denominator in Eq. (6) can be
combined with an utterance-wise regularizer to evaluate the utterance
reliability, and the third selection criterion is formulated as

$\displaystyle\frac{\mathrm{PLDA}^{\prime}(\bm{x}_{n,i},\bm{x}_{Target})}{(\sigma_{n}\lVert\bm{x}_{n,i}-\bm{u}_{n}\rVert_{2})^{\alpha}},$
(7)

where $\left\lVert\cdot\right\rVert_{2}$ denotes the Euclidean distance (L2
norm), and $\bm{u}_{n}$ is the mean x-vector of speaker $n$. The larger the
Euclidean distance, the lower the reliability.

In summary, the criteria in this section model different relationships among
the target speaker and the individual speakers and utterances in the candidate
pool. Each subsequent criterion is derived from the previous criterion.

##  4 Experiments

###  4.1 Corpus

The AIShell-3 and TST Mandarin corpora provided by the ICASSP2021-M2VoC
organizer [37] were used in the experiments. The training set of AIShell-3,
which includes 137 female and 37 male speakers, was used to train the baseline
SI vocoder and was used as the candidate pool. The female speakers have 50,117
utterances ($\sim$50 hours) in total, and the male speakers have 13,145
utterances ($\sim$13 hours) in total. One female and two male speakers in the
Track 1 subset of TST were used as the target speakers. Each target speaker
has 90 training utterances (6–10 mins) and 10 testing utterances. To simulate
the one-shot scenario, the first five training utterances ($\sim$30 s) of each
target speaker were used as the limited target data. The $f_{s}$ and bit-depth
of all utterances were set to 44.1 KHz and 16.

###  4.2 Experimental setting

The mel-spectrogram of 80 mel-filter banks was used as the input of the
vocoder. The hop size was 220 samples, and the FFT size was 2048. The pre-
trained models of the SITW (speakers in the wild) x-vector system222
https://kaldi-asr.org/models/m8was used for extraction of 512-dimensional
x-vectors and calculation of PLDA scores. All input audio files were
downsampled to 16 KHz to match the working $f_{s}$ of these speaker models.
The x-vector of each candidate speaker was the average x-vector of all the
utterances from that speaker, and the x-vector of each target speaker was also
the average x-vector of the available utterances corresponding to different
scenarios.

###  4.3 Model description

Six multi-band MelGAN vocoders were evaluated. That is, an SI (multi-speaker)
vocoder was first trained using the training set of AIShell-3, and then the
adapt5 and adapt90 SD vocoders were developed by adapting the SI vocoder using
five and 90 target utterances, respectively. The adapt90 vocoder was taken as
an upper bound in this section. For the proposed vocoders, to match the amount
of adaptation data of adapt90, 85 utterances were selected from the AIShell-3
training set using the proposed criteria 1–3, and then combined with the five
target utterances to form the adaptation sets DC1–DC3, where DC denotes the
data selection criterion. The proposed SD vocoders were adapted from the SI
vocoder using the DC1–DC3 sets, respectively.

The SI vocoder was trained for 1M iterations, and the discriminators were
jointly trained with the generator from the 200,000-th iteration. Adam
optimizer [38] was used, and the learning rate was set to 0.001 without decay.
During speaker adaptation, the discriminators and generator were updated. The
adaptation iteration number for the adapt5 vocoders was set to 1000, and the
iteration number for the adapt90 and DC1–DC3 vocoders was set to 9000. It was
difficult to find the optimal number of iterations for the adapt5 vocoders,
because the adaptation data was very limited. Therefore, the iteration number
of 1000 was a compromise choice for the three target speakers.

###  4.4 Objective evaluation

Table 1: Objective evaluation results. |  SI |  Adapt90 |  Adapt5 |  DC1 |  DC2 |  DC3
---|---|---|---|---|---|---
LSD (dB) |  1.07 |  1.00 |  1.08 |  1.09 |  1.08 |  1.04
MCD (dB) |  6.36 |  5.43 |  5.96 |  6.06 |  6.03 |  6.08
$F_{0}$ (Hz) |  70.8 |  67.3 |  69.5 |  68.2 |  67.9 |  70.0
$U/V$ (%) |  15.6 |  13.7 |  14.9 |  15.1 |  14.9 |  14.7
PLDA |  33.2 |  39.5 |  29.0 |  34.8 |  32.6 |  33.6
CosSim |  0.82 |  0.88 |  0.81 |  0.84 |  0.84 |  0.84

Objective evaluations based on spectral accuracy, source excitation accuracy,
and speaker similarity were conducted. The spectral accuracy was evaluated in
terms of log spectral distortion (LSD) and mel-cepstral distortion (MCD). For
source excitation accuracy, we measured the root mean square error (RMSE) of
$F_{0}$ and $U/V$ decision error. For speaker similarity measurement, the PLDA
score and cosine similarity (CosSim) of two x-vectors were used. Mel-cepstrum
($mcep$), $F_{0}$, and unvoiced/voicded ($U/V$) features were extracted using
WORLD [3]. The ground truth acoustic features and x-vectors were extracted
from the natural testing utterances.

The average evaluation results of the three target speakers are shown in Table
1. As expected, the adapt90 vocoders achieve the best performance in all
metrics, which shows the effectiveness of speaker adaptation when the
adaptation data is relatively sufficient. In contrast, the performance of the
adapt5 vocoders is much worse than that of the adapt90 vocoders. This is due
to the instability and quality degradation caused by much less adaptation
data.

For spectral accuracy evaluation, we can see that the adapt5 vocoders were
better than the SI vocoder in MCD, but worse in LSD. One possible reason is
that $mcep$ is dominated by spectral envelope components in low-frequency
bands, and due to improved SD component modeling, the adapt5 vocoders can
achieve higher formant modeling accuracy. When listening to the utterances
produced by adapt5, we could perceive a similar trend. That is, despite the
slight improvement in the similarity of timbre, the speech generated by adapt5
suffered from severe musical noise and oversmoothing. The musical noise and
oversmoothing effects may not be well modeled by $mcep$, but will be reflected
in the LSD measurement. Moreover, the spectral accuracy results show that the
the proposed methods are effective because the DC* vocoders achieve lower MCD
and similar/lower LSD than the SI vocoder. Compared with the adapt5 vocoders,
the DC3 vocoders achieve lower LSD and slightly higher MCD, which implies that
external data may slightly reduce the accuracy of SD component modeling.

Since both the voiced and unvoiced parts were involved in the RMSE calculation
of $F_{0}$, the error of $F_{0}$ is correlated with the $U/V$ error. Generally
speaking, prosodic characteristics are highly related to speaker identity. As
shown in Table 1, the SI vocoder yielded higher $F_{0}$ and $U/V$ errors than
the SD vocoders. The results confirm the assumption and indicate that prosodic
modeling can be improved by speaker adaptation. Moreover, the DC* vocoders
achieved similar or lower prosodic errors to the adapt5 vocoders, indicating
the potential of the proposed methods to further improve prosodic modeling.

For speaker similarity evaluation, the DC* vocoders achieved higher PLDA and
CosSim scores than the SI and adapt5 vocoders, which shows the effectiveness
of the proposed method in improving the speaker similarity of synthetic
speech. Again, one possible reason for the worst performance of the adapt5
vocoders is that the musical noise and oversmoothing can reduce the speaker
similarity. In addition, since DC2 and DC3 apply regularizations to the PLDA
score, it is reasonable that the DC1 vocoders are slightly better than the DC2
and DC3 vocoders in the PLDA score, but CosSim remains the same.

###  4.5 Subjective evaluation

Table 2: Subjective evaluation results (MOS values).

|  Natural |  SI |  Adapt90 |  Adapt5 |  DC3
---|---|---|---|---|---
Quality |  4.96$\pm$.03 |  3.47$\pm$.14 |  3.71$\pm$.13 |  2.65$\pm$.14 |  3.29$\pm$.13
Similarity |  - |  3.57$\pm$.16 |  4.23$\pm$.14 |  3.43$\pm$.16 |  3.71$\pm$.15

Table 3: Statistics of the selected utterances (female target).

|  Number of speakers |  Suspected utterances |  Utterance overlap (%) |  Speaker overlap (%)
---|---|---|---|---
DC1 |  18 |  9 |  reference |  reference
DC2 |  7 |  2 |  15.3 |  32.0
DC3 |  16 |  5 |  10.6 |  70.6

Two mean opinion score (MOS) evaluations were conducted for speech quality and
speaker similarity, respectively. In the quality evaluation, each utterance
was given a score in the range 1–5 by a listener to evaluate the naturalness.
In the similarity evaluation, the listener listened to a pair of natural and
synthetic utterances at a time, and gave a score between 1 and 5 to evaluate
the speaker similarity of the synthetic utterance to the natural utterance. In
both evaluations, the higher the score, the better the performance. Nine
native speakers participated in the evaluation using the same device. Five
types of speech were compared, including natural speech and the speech
produced by the SI, adapt5, adapt90, and DC3 vocoders. For each of the three
target speakers, there were 10 natural utterances and 40 synthetic utterances
produced by four vocoders. The 150 utterances were divided into two subsets,
and each subset was evaluated by at least five listeners. Demo samples can be
found on our website [39].

As shown in Table 3, the superior performance of the adapt90 vocoder in both
quality and similarity measurements proves the importance of speaker
adaptation to the vocoder. However, due to severe musical noise and
oversmoothing, the adapt5 vocoder achieved the worst naturalness and
similarity. It even gave worse performance than the SI vocoder. The results
show that one-shot speaker adaptation is challenging, and the vocoders adapted
with extremely limited data tend to be unstable and cannot be generalized.
However, the proposed DC3 vocoder significantly outperformed the adapt5
vocoder in both quality and similarity measurements. The results confirm the
effectiveness of the proposed data augmentation method derived from SV
technology for the speaker adaptation of SD vocoders.

Although the DC3 vocoder is superior to the SI vocoder in terms of speaker
similarity, the SI vocoder is robust to unseen speakers in terms of quality.
The result is reasonable because the SI vocoder is trained with a large amount
of data from many speakers. The great diversity of the training data allows
the SI vocoder to generate stable speech for unseen speakers, even if the
similarity is still insufficient. Therefore, we may conclude that the quality
of synthetic speech is highly related to the amount of training/adaptation
data, and the proposed data augmentation method can improve the speaker
similarity of the SI vocoder even when only 30s target data is available.
Moreover, according to the results, there is a significant quality gap between
natural and synthetic speech, which implies that the current multi-band MelGAN
vocoder may not be able to handle high $f_{s}$ speech generation. There is
still room for improving neural vocoders to generate signals with a high
$f_{s}$.

###  4.6 Discussion

Figure 1: PLDA score distributions of candidate utterances.

The score distributions for the female and the first male targets in Fig. 1
show the PLDA scores reflect the high correlation between gender and speaker
similarity as expected. The selection thresholds of different utterance
numbers imply the challenge to optimize the selection numbers because of the
rapid PLDA score decrease for the increased selected utterances.

Since the female pool is much larger, we report the statistics of the
selection utterances for the female target using DC1–3 in Table 3, and the
results of DC1 are taken as the reference. The speaker regularizer is adopted
to filter the utterances of low confidence speakers, and the speaker number of
the selected utterances does decrease in DC2. An utterance might be an outlier
if it is the only utterance from a speaker to be selected, and we call it a
suspected utterance. The utterance regularizer is utilized to remove the
outlier utterances, and the suspected utterance number does reduce in DC3.

The low speaker overlap rate between the DC1 and DC2 sets means the speaker
regularizer filters many speakers. The higher speaker but lower utterance
overlap rates between the DC1 and DC3 sets show that the utterance regularizer
makes DC3 select more representative utterances of each speaker, and the
influence of the speaker regularization is eased in DC3.

##  5 Conclusions

In this paper, we first investigated the influence of speaker adaptation and
adaptive data volume on the effectiveness of SD neural vocoder. Then, we
proposed selection criteria based on SV technology to leverage an external
corpus for improving speaker adaptation. The proposed criteria explored the
different relationships in the data space. The evaluation results show the
effectiveness of the proposed framework in one-shot speech synthesis. For
future work, the amount of selected data and the weight of regularizers should
be further optimized.

##  6 Acknowledgements

This work was supported by JSPS KAKENHI Grant Number 17H06101 and JST, CREST
Grant Number JPMJCR19A3.

## References

  * [1] H. Dudley, “Remaking speech,” _The Journal of the Acoustical Society of America_ , vol. 11, no. 2, pp. 169–177, 1939.
  * [2] H. Kawahara, I. Masuda-Katsuse, and A. De Cheveigne, “Restructuring speech representations using a pitch-adaptive time–frequency smoothing and an instantaneous-frequency-based F0 extraction: Possible role of a repetitive structure in sounds,” _Speech Communication_ , vol. 27, no. 3-4, pp. 187–207, 1999.
  * [3] M. Morise, F. Yokomori, and K. Ozawa, “WORLD: a vocoder-based high-quality speech synthesis system for real-time applications,” _IEICE Transactions on Information and Systems_ , vol. 99, no. 7, pp. 1877–1884, 2016\.
  * [4] R. McAulay and T. Quatieri, “Speech analysis/synthesis based on a sinusoidal representation,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 34, no. 4, pp. 744–754, 1986.
  * [5] A. van den Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, “WaveNet: A generative model for raw audio,” in _Proc. SSW9_ , Sept. 2016, p. 125.
  * [6] S. Mehri, K. Kumar, I. Gulrajani, R. Kumar, S. Jain, J. Sotelo, A. Courville, and Y. Bengio, “SampleRNN: An unconditional end-to-end neural audio generation model,” in _Proc. ICLR_ , Apr. 2017.
  * [7] N. Kalchbrenner, E. Elsen, K. Simonyan, S. Noury, N. Casagrande, E. Lockhart, F. Stimberg, A. van den Oord, S. Dieleman, and K. Kavukcuoglu, “Efficient neural audio synthesis,” in _Proc. ICML_ , July 2018, pp. 2415–2424.
  * [8] A. van den Oord, Y. Li, I. Babuschkin, K. Simonyan, O. Vinyals, K. Kavukcuoglu, G. van den Driessche, E. Lockhart, L. C. Cobo, F. Stimberg, N. Casagrande, D. Grewe, S. Noury, S. Dieleman, E. Elsen, N. Kalchbrenner, H. Zen, A. Graves, H. King, T. Walters, D. Belov, and D. Hassabis, “Parallel WaveNet: Fast high-fidelity speech synthesis,” in _Proc. ICML_ , July 2018, pp. 3915–3923.
  * [9] W. Ping, K. Peng, and J. Chen, “ClariNet: Parallel wave generation in end-to-end text-to-speech,” in _Proc. ICLR_ , May 2019.
  * [10] R. Prenger, R. Valle, and B. Catanzaro, “WaveGlow: A flow-based generative network for speech synthesis,” in _Proc. ICASSP_ , May 2019, pp. 3617–3621.
  * [11] S. Kim, S.-G. Lee, J. Song, J. Kim, and S. Yoon, “FloWaveNet : A generative flow for raw audio,” in _Proc. ICML_ , June 2019, pp. 3370–3378.
  * [12] R. Yamamoto, E. Song, and J.-M. Kim, “Parallel WaveGAN: A fast waveform generation model based on generative adversarial networks with multi-resolution spectrogram,” in _Proc. ICASSP_ , May 2020, pp. 6199–6203.
  * [13] K. Kumar, R. Kumar, T. de Boissiere, L. Gestin, W. Z. Teoh, J. Sotelo, A. de Brébisson, Y. Bengio, and A. C. Courville, “MelGAN: Generative adversarial networks for conditional waveform synthesis,” in _Proc. NeurIPS_ , Dec. 2019, pp. 14 910–14 921.
  * [14] M. Bińkowski, J. Donahue, S. Dieleman, A. Clark, E. Elsen, N. Casagrande, L. C. Cobo, and K. Simonyan, “High fidelity speech synthesis with adversarial networks,” in _Proc. ICLR_ , Apr. 2020.
  * [15] A. Tamamori, T. Hayashi, K. Kobayashi, K. Takeda, and T. Toda, “Speaker-dependent WaveNet vocoder,” in _Proc. INTERSPEECH_ , Aug. 2017, pp. 1118–1122.
  * [16] T. Hayashi, A. Tamamori, K. Kobayashi, K. Takeda, and T. Toda, “An investigation of multi-speaker training for WaveNet vocoder,” in _Proc. ASRU_ , Dec. 2017, pp. 712–718.
  * [17] Y.-C. Wu, P. L. Tobing, T. Hayashi, K. Kobayashi, and T. Toda, “The NU non-parallel voice conversion system for the Voice Conversion Challenge 2018,” in _Proc. Odyssey_ , June 2018, pp. 211–218.
  * [18] P. L. Tobing, Y.-C. Wu, T. Hayashi, K. Kobayashi, and T. Toda, “NU voice conversion system for the Voice Conversion Challenge 2018,” in _Proc. Odyssey_ , June 2018, pp. 219–226.
  * [19] J. Lorenzo-Trueba, T. Drugman, J. Latorre, T. Merritt, B. Putrycz, R. Barra-Chicote, A. Moinet, and V. Aggarwal, “Towards achieving robust universal neural vocoding,” in _Proc. Interspeech 2019_ , 2019, pp. 181–185. [Online]. Available: http://dx.doi.org/10.21437/Interspeech.2019-1424
  * [20] A. Tjandra, S. Sakti, and S. Nakamura, “Machine speech chain with one-shot speaker adaptation,” in _Proc. Interspeech 2018_ , 2018, pp. 887–891. [Online]. Available: http://dx.doi.org/10.21437/Interspeech.2018-1558
  * [21] A. Rosenberg, Y. Zhang, B. Ramabhadran, Y. Jia, P. Moreno, Y. Wu, and Z. Wu, “Speech recognition with augmented synthesized speech,” in _ASRU_. IEEE, 2019, pp. 996–1002.
  * [22] A. Laptev, R. Korostik, A. Svischev, A. Andrusenko, I. Medennikov, and S. Rybin, “You do not need more data: improving end-to-end speech recognition by text-to-speech data augmentation,” in _CISP-BMEI_. IEEE, 2020, pp. 439–444.
  * [23] Y. Jia, M. Johnson, W. Macherey, R. J. Weiss, Y. Cao, C.-C. Chiu, N. Ari, S. Laurenzo, and Y. Wu, “Leveraging weakly supervised data to improve end-to-end speech-to-text translation,” in _ICASSP_. IEEE, 2019, pp. 7180–7184.
  * [24] M. Airaksinen, L. Juvela, P. Alku, and O. Räsänen, “Data augmentation strategies for neural network f0 estimation,” in _ICASSP_. IEEE, 2019, pp. 6485–6489.
  * [25] H. Yamamoto, K. A. Lee, K. Okabe, and T. Koshinaka, “Speaker augmentation and bandwidth extension for deep speaker embedding,” in _Proc. Interspeech 2019_ , 2019, pp. 406–410. [Online]. Available: http://dx.doi.org/10.21437/Interspeech.2019-1508
  * [26] M.-J. Hwang, R. Yamamoto, E. Song, and J.-M. Kim, “TTS-by-TTS: TTS-driven data augmentation for fast and high-quality speech synthesis,” _arXiv preprint arXiv:2010.13421_ , 2020.
  * [27] D. Snyder, D. Garcia-Romero, G. Sell, D. Povey, and S. Khudanpur, “X-vectors: Robust dnn embeddings for speaker recognition,” in _ICASSP_. IEEE, 2018, pp. 5329–5333.
  * [28] G. Yang, S. Yang, K. Liu, P. Fang, W. Chen, and L. Xie, “Multi-band MelGAN: Faster waveform generation for high-quality text-to-speech,” in _Proc. SLT_ , Jan. 2021.
  * [29] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in _Proc. NIPS_ , Dec. 2014, pp. 2672–2680.
  * [30] F. Yu and K. Vladlen, “Multi-scale context aggregation by dilated convolutions,” in _Proc. ICLR_ , May 2016.
  * [31] A. L. Maas, A. Y. Hannun, and A. Y. Ng, “Rectifier nonlinearities improve neural network acoustic models,” in _Proc. ICML_ , June 2013, pp. 3–11.
  * [32] T. Okamoto, K. Tachibana, T. Toda, Y. Shiga, and H. Kawai, “An investigation of subband WaveNet vocoder covering entire audible frequency range with limited acoustic features,” in _ICASSP_. IEEE, 2018, pp. 5654–5658.
  * [33] T. Okamoto, T. Toda, Y. Shiga, and H. Kawai, “Improving FFTNet vocoder with noise shaping and subband approaches,” in _SLT_. IEEE, 2018, pp. 304–311.
  * [34] C. Yu, H. Lu, N. Hu, M. Yu, C. Weng, K. Xu, P. Liu, D. Tuo, S. Kang, G. Lei _et al._ , “DurIAN: Duration informed attention network for multimodal synthesis,” _arXiv preprint arXiv:1909.01700_ , 2019.
  * [35] S. Ioffe, “Probabilistic linear discriminant analysis,” in _ECCV_. Springer, 2006, pp. 531–542.
  * [36] S. J. Prince and J. H. Elder, “Probabilistic linear discriminant analysis for inferences about identity,” in _ICCV_. IEEE, 2007, pp. 1–8.
  * [37] Q. Xie, X. Tian, G. Liu, K. Song, L. Xie, Z. Wu, H. Li, S. Shi, H. Li, F. Hong, H. Bu, and X. Xu, “The multi-speaker multi-style voice cloning challenge 2021,” in _ICASSP_. IEEE, 2021.
  * [38] D. P. Kingma and J. L. Ba, “Adam: A method for stochastic optimization,” in _Proc. ICLR_ , May 2015.
  * [39] Y.-C. Wu, _RDS demo_ , Accessed: 2021. [Online]. Available: https://bigpon.github.io/RelationalDataSelection_demo/