Affine equivariant Tyler's M-estimator applied to tail parameter learning of elliptical distributions

Esa Ollila, Senior Member, IEEE, Daniel P. Palomar, Fellow, IEEE, and Frédéric Pascal Senior Member, IEEE

Esa Ollila is with the Department of Information and Communications Engineering, Aalto University, FI-00076 Aalto, Finland (e-mail: esa.ollila@aalto.fi). Daniel P. Palomar is with the Hong Kong University of Science and Technology, Hong Kong (e-mail: palomar@ust.hk). Frédéric Pascal is with the CNRS Université Paris-Saclay, CentraleSupélec, 91190 Gif-sur-Yvette, France (e-mail:frederic.pascal@centralesupelec.fr).

Affine equivariant Tyler's M-estimator
Ollila et al.

We propose estimating the scale parameter (mean of the eigenvalues) of the scatter matrix of an unspecified elliptically symmetric distribution using weights obtained by solving Tyler's M-estimator of the scatter matrix. The proposed Tyler's weights-based estimate (TWE)  of scale is then used to construct an affine equivariant Tyler's M-estimator as a weighted sample covariance matrix using normalized Tyler's weights. We then develop a unified framework for estimating the unknown  tail parameter of the elliptical distribution (such as the degrees of freedom (d.o.f.) $\nu$ of the multivariate $t$ (MVT) distribution). Using the proposed TWE of scale, a new robust estimate of the d.o.f. parameter of MVT distribution is proposed with excellent performance in heavy-tailed scenarios,  outperforming other competing methods. R-package is available that implements the proposed method.

Tyler's M-estimator, scatter matrix, covariance matrix, elliptical distributions

§ INTRODUCTION

We model the observed $p$-variate observations $\x_1,\ldots,\x_n$ as independent and identically distributed (i.i.d.) random samples from an unspecified centered  (i.e., symmetric around the origin) elliptically symmetric (ES) distribution [1, 2]. A continuous random vector $\x \in \mathbb{R}^p$ has centered ES distribution if
it possesses  a probability density function (pdf) of the form

\begin{equation*}
f(\x)=C_{p,g}    |\M|^{-1/2} g ( \x^\top \M^{-1} \x),
\end{equation*}

where $g:\mathbb{R}_{\geq 0} \to \mathbb{R}_{> 0}$ is called the density
generator,  $\M \succ 0$ is the positive definite symmetric matrix parameter, called the scatter matrix, and $C_{p,g}$ is a normalizing constant ensuring that $f(\x)$ integrates to 1. We let $\x \sim \mathcal E_p(\mathbf{0},\M,g)$ to denote this case.  For example, the centered multivariate normal (MVN) distribution  $\mathcal N_p(\mathbf{0},\M)$ is obtained when
$g(t) = \exp(-t/2)$ while the multivariate $t$ (MVT) distribution with $\nu>0$ degrees of freedom (d.o.f.)  is  obtained when

g(t)=(1+  t/ν)^-(p+ν)/2.
Parameter $\nu>0$ is a tail parameter of the density.
For $\nu \to \infty$, the MVT distribution reduces to the MVN distribution, while $\nu=1$ corresponds to the multivariate Cauchy distribution.  Also, many other subclasses of ES distributions are parametrized by a density generator indexed by an additional tail parameter that is unknown in practice and needs to be estimated. Learning this unknown parameter is also one of the goals of this paper.

We are interested in estimating the scale parameter of the scatter matrix, defined as the mean of its eigenvalues,

η= ()/p = 1/p ∑_i=1^p λ_i ,
where $\lambda_i >0$ denotes the $i$th eigenvalue of $\M$.
Formally, $\eta \equiv \eta(\M)$ is  a scale parameter if it verifies $\eta(\mathbf{I}) = 1$ and  $\eta(a\boldsymbol{\Sigma}) = a \eta(\boldsymbol{\Sigma})$ for all $a>0$ [3].

Tyler's M-estimator [4] is a popular robust M-estimator of the scatter matrix that has been extensively studied both in signal processing and statistics literature (e.g, [5, 6, 7, 8, 9, 10, 11, 12, 13]). Tyler's M-estimator is
defined as the solution to the fixed-point equation

=  1/n ∑_i=1^n  p _i^⊤^-1 _i     _i _i^⊤≜ℋ(; {_i}).
Note that me way also write the map $\mathcal H(\cdot; \cdot)$ in the form

\mathcal H(\hat \M ; \{\x_i\}) = \frac{1}{n} \sum_{i=1}^n  \hat w_i   \x_i \x_i^\top  \  \mbox{ with } \
\hat w_i = \cfrac{p}{ \x_i^\top \hat \M^{-1} \x_i} ,

where $\hat w_i$, $i=1,\ldots, n,$ are referred to as Tyler's weights.
Tyler's M-estimator  is  unique only up to a scaling factor, and therefore a common convention  is to consider a solution that verifies  $\tr(\hat{\M})=p$.  Thus, Tyler's M-estimator is actually an estimator of a shape matrix (normalized scatter matrix) $\boldsymbol{\Lambda}$, defined by
\boldsymbol{\Lambda} =  \M/\eta = p \M/\tr(\M),
and verifying $\tr(\boldsymbol{\Lambda})=p$.

In this paper, we propose an estimator of the scale $\eta$ based on Tyler's weights $\hat w_i$. The proposed scale estimate along with  Tyler's M-estimator $\hat \M$ are then jointly used for constructing affine equivariant robust  estimates of the scatter matrix $\M$ and the covariance matrix $\mathbf{R}=\cov(\x)$ (or their shrinkage versions). These developments are described in <ref>.  Then, in <ref>, we propose a unified framework allowing to estimate the tail parameter of the elliptical distribution using the proposed scale statistic $\hat \eta$. In the case of the MVT distribution, this leads to a new estimate of the d.o.f. parameter based on Tyler's weights. Finally, Section <ref> demonstrates the relevance of the proposed approach on simulated data, with concluding remarks in Section <ref>. In the R package  [14], the function  implements this method.

§ ESTIMATE OF SCALE, SCATTER, AND COVARIANCE MATRIX BASED ON TYLER'S WEIGHTS

Assuming that $\x \sim   \mathcal E_p(\mathbf{0},\boldsymbol{\Sigma},g)$  has finite 2nd-order moments, then its covariance matrix, $\mathbf{R} =  \E[ \x\x^\top]$ satisfies

𝐑 = θ·Σ      θ=  [ r^2] /p ,
where  $r^2 = \|\M^{-1/2} \x\|^2$ is the squared Mahalanobis distance of $\x$ w.r.t. $\M$, whose pdf is given by

f_r^2(t) = C t^p/2 -1 g(t) ,
where $C= \int_0^\infty t^{p/2 -1} g(t)  \mathrm{d} t$. Hence pdf of $r^2$ has a one-to-one correspondence with density generator $g$.
From (<ref>), we notice that the scatter matrix $\boldsymbol{\Sigma}$ is proportional to the covariance matrix $\mathbf{R}$ (assuming $\mathbf{R}$  exists).
In the MVN case,  $\theta=1$, while for the MVT distribution with density generator as in (<ref>) one obtains $\theta= \nu/(\nu-2)$ for all $\nu >2$.

§.§ Estimate of scale

As mentioned earlier, Tyler's M-estimator $\hat \M$ loses information of the scale $\eta$.  However, it is yet possible to construct an estimate of $\eta$   from  Tyler's weights.
Our Tyler's weights-based estimate (TWE) of scale is defined as  the harmonic mean of reciprocal of weights, $1/\hat w_i$'s, that is,

\begin{equation} \label{eq:hat_tau3}
\hat \eta_{\text{TWE}} =  \bigg( \frac{1}{n} \sum_{i=1}^n \hat w_i  \bigg)^{-1}   =  \bigg( \frac{p}{n} \sum_{i=1}^n  [\x_{i}^\top \hat{\M}^{-1} \x_{i}]^{-1} \bigg)^{-1} .
\end{equation}

One can esily verify that this is a proper scale estimate  in the sense that if $\hat \eta_{\text{TWE}}^*$ is computed on scaled observations, $\x_i^* = c \cdot\x_i$, $i=1,\ldots, n$, then $\hat \eta_{\text{TWE}}^* = c^2  \cdot \hat \eta_{\text{TWE}}$.  This follows because Tyler's M-estimator $\hat \M$ in (<ref>)  with trace constraint is invariant to scaling the data, so $\hat \M^* = \hat \M$.

The proposed estimate (<ref>) can also be motivated from the following result  derived in the high-dimensional random matrix theory (RMT) regime, where $p,n \to \infty$ with $n>p$ and  their ratio tending to constant: $
p/n \to c \in (0,1).$
Namely, let $ \hat{\M}$ be Tyler's M-estimator in (<ref>) verifying $\mathrm{tr}(\hat \M)=p$. Then, it was shown in [11, 13] that $\max_{\ell}  | \eta \hat w_{\ell} - 1 | \to 0$ almost surely.  The authors in [11] derived this result for the case that data is i.i.d. Gaussian $\mathcal N_p(\mathbf{0}, \mathbf{I})$ while [13] extended these results for more general distributions.  Thus since $1/\hat w_{i}$ concentrate on $\eta$,
(<ref>) is a natural robust estimator of scale.  Many other robust scale statistics could be constructed from $1/\hat w_i$'s, such as the median, trimmed mean, etc. In the next subsection, we illustrate why the proposed harmonic mean (<ref>) is the most natural.

§.§ Affine equivariant estimate of scatter matrix

Using the scale estimate $\hat \eta_{\text{TWE}}$ and Tyler's M-estimator $\hat \M$ (with $\tr(\hat \M)=p$), we can form an estimate of the scatter matrix

\begin{align} \label{eq:Tyler_Scatter}
\hat \M_{\text{TWE}}  &= \hat \eta_{\textup{TWE}}   \cdot \hat \M
\end{align}

referred  to as TWE of scatter matrix.
Thus $\hat \eta_{\textup{TWE}}$ is scale statistic derived from $\hat \M_{\text{TWE}}$ since $\hat \eta_{\textup{TWE}}= \tr(\hat \M_{\text{TWE}})/p$. Equivalently, the trace of $\hat \M_{\text{TWE}}$ can be easily computed as the harmonic mean of Tyler's quadratic form:

$$\tr(\hat \M_{\text{TWE}}) = \bigg( \frac{1}{n} \sum_{i=1}^n  [\x_{i}^\top \hat{\M}^{-1} \x_{i}]^{-1} \bigg)^{-1}.$$

Recalling (<ref>) we can write (<ref>) in the following more intuitive form:

\begin{align}
\hat \M_{\text{TWE}}   &=   \frac{1}{n}\sum_{i=1}^n   \hat v_i \x_i \x_i^\top,  \label{eq:TylerScatter2a}  \\
\hat v_{i} &=  \frac{\hat w_{i}}{\frac{1}{n}\sum_{\ell=1}^n \hat w_\ell } = \frac{  [\x_{i}^\top \hat \M^{-1} \x_{i}]^{-1} } { \frac{1}{n}\sum_{\ell=1}^n  [\x_\ell^\top \hat \M^{-1} \x_\ell]^{-1} } \label{eq:TylerScatter2b}
\end{align}

where $\hat v_1,\ldots,\hat v_n$ are normalized Tyler's weights that verify $\frac{1}{n}\sum_{i =1}^n \hat v_i = 1$.  Eq. (<ref>) and (<ref>)  illustrate that $\hat \M_{\text{TWE}}$ is a weighted sample covariance matrix (SCM) with weights $\hat v_{i}$.

Finally, we draw the parallel of Tyler's M-estimating equation and our estimator (<ref>). First, note that Tyler's M-estimating equation (<ref>) verifies broader invariance than just invariance with respect to scaling of the data matrix. Namely, denoting the unit-norm normalized data by $\tilde \x_i = \x_i/\| \x_i \|$, $i=1,\ldots,n$,  one can easily verify that the fixed-point equation in
(<ref>) can be rewritten as
$\hat \M  = \mathcal H(\hat \M ; \{\tilde \x_i\})$, so based on normalized data.
Furthermore, since $\tr(\hat \M) =p $, one has:

\begin{equation*}
\tr(\hat \M) = \frac{1}{n} \sum_{i=1}^n  \cfrac{p}{ \tilde\x_i^\top \hat{\M}^{-1} \tilde\x_i}  \,  \tr\left(\tilde\x_i \tilde\x_i^\top\right) = p,
\end{equation*}

or equivalently $\frac{1}{n} \sum_{i=1}^n  [\tilde\x_i^\top \hat{\M}^{-1} \tilde\x_i]^{-1} = 1$ since $\tr\left(\tilde\x_i \tilde\x_i^\top\right)=1.$ It follows that Tyler's M-estimator $\hat \M$ with $\tr(\hat \M)=p$ is the solution to the following fixed-point equation:

= 1/n∑_i=1^n w̃_i _i _i^⊤  w̃_i
=   p   [_i^⊤^-1 _i]^-1 / 1/n∑_ℓ=1^n  [_ℓ^⊤^-1 _ℓ]^-1 ,
where $\tilde{w}_{\ell}$'s are the normalized Tyler's weights computed on normalized (unit norm) observations $\tilde \x_i$'s. Thus, while Tyler's M-estimator $\hat \M$ can be interpreted as a weighted SCM based on normalized data $\{\tilde \x_i\}$ as shown in (<ref>),  TWE of scatter $\hat \M_{\textup{TWE}}$ in (<ref>) can be viewed as weighted SCM of actual (non-normalized) data $\{\x_i\}$.

It is worthwhile to point out that most robust estimators of scatter  are affine equivariant in the sense that an affine transformation on the data $\mathbf{x}_i  \mapsto \mathbf{A} \mathbf{x}_i$, $i = 1,\ldots,n$, induces following transformation on the estimate:

({ _i})= ({_i}) ^⊤,  ∀𝐀 ∈ℝ^p ×p
For example, robust Maronna's [15] M-estimators, S-estimators [16], or MM-estimators [17] are affine equivariant.
However, Tyler's scatter matrix $\hat \M$ is not affine equivariant since   (<ref>) only holds up to multiplicative scalar factor as shown in (<ref>).
Affine equivariance is desirable since if $\x \sim \mathcal E_p(\mathbf{0},\M, g) $, then $\A \x \sim \mathcal E_p(\mathbf{0},\A \M \A^\top, g)$.
Hence the transformed data shares the same elliptical distribution, but the scatter matrix parameter is mapped to $\A \M \A^\top$.
Thus a natural requirement to be imposed on any scatter matrix estimator is that it should verify this same equivariance  principle under transformations $\x_i \mapsto \A \x_i$.
This is shown next.

TWE of scatter matrix $\hat \M_{\textup{TWE}}$ is affine equivariant:  $\hat \M_{\textup{TWE}}(\{ \A \x_i \})= \A \hat \M(\{\x_i\}) \A^\top$,  $\forall$ invertible $\mathbf{A} \in \mathbb{R}^{p \times p}$.

It is straightforward to verify that Tyler's M-estimator (with $\tr(\hat \M)=p$) is equivariant in the sense that if $\hat \M^* = \hat{\M}(\{ \x_i^*\})$ denotes Tyler's M-estimator (verifying $\tr(\hat \M^*) =p$)  computed on data $\x_i^*=\A \x_i$, $i=1,\ldots,n$, then

^* =c  ^⊤,   c= p/ (^⊤)
for all invertible $\A \in \mathbb{R}^{p\times p}$.
Now let $\hat w_i^*$, $i=1,\ldots,n$   denote the corresponding Tyler's weights.
Then the scale estimate (<ref>) computed as the harmonic mean of  reciprocals of weights $1/\hat w_i^* = [ \x_i^*]^\top(\hat{\M}^*)^{-1} \x_i^*$ is

\begin{align}
\hat \eta^*_{\textup{TWE}} &= \Big(\frac 1 n \sum_{i=1}^n \hat w_i^* \Big)^{-1}
=  \bigg( \frac{p}{n} \sum_{i=1}^n  [\x_{i}^\top \A^\top (\hat{\M}^*)^{-1} \A \x_{i}]^{-1} \bigg)^{-1}  \notag \\
&= c^{-1} \bigg( \frac{p}{n} \sum_{i=1}^n  [\x_{i}^\top \hat{\M}^{-1}  \x_{i}]^{-1} \bigg)^{-1}  = c^{-1} \hat{\eta}_{\textup{TWE}}
\end{align}

where in the 2nd to last identity we simply utilized (<ref>). Thus it follows that

\begin{align*}
\hat{\M}_{\textup{TWE}}^* &= \hat \eta_{\textup{TWE}}^* \cdot \hat \M^* = c^{-1} \hat \eta_{\textup{TWE}} \cdot  c  \A \hat{\M}\A^{\top} \\ &= \A ( \hat \eta_{\textup{TWE}}\hat \M)\A^\top
= \A \hat\M_{\textup{TWE}} \A^{\top}.
\end{align*}

In the case of large dimensional data, one can also consider a shrinkage Tyler's M-estimator of the scatter matrix as

_,β = β_  + (1-β) η̂_ I,
where the data adaptive shrinkage parameter $\beta \in [0,1]$ is computed as described in <cit.>. However, unlike the estimator
in [18], the shrinkage TWE in  (<ref>) provides  an estimator  of scatter instead of shape matrix.

§.§ An estimator of covariance matrix

If the  density generator $g$ (and hence the underlying ES distribution) is specified, then the value of  $\theta$ in (<ref>) can be determined,  and  we can use relationship (<ref>) to obtain a covariance matrix estimator as
\hat{\mathbf{R}}_{\text{TWE}}  = \theta \cdot \hat \M_{\text{TWE}} .
For example, if the data has an MVN distribution,  then
$\theta=1$ while $\theta = \nu/(\nu-2)$ in the case of an MVT distribution with $\nu$ d.o.f.
However, often the underlying parametric family is known, but the underlying tail parameter, say $\nu$, indexing the density generator is unknown.
As is shown in <ref>, we can form an estimate of $\nu$, denoted $\hat \nu_{\text{TWE}}$, using Tyler's weights.
Since $\theta = h(\nu)$ (cf. Eq. (<ref>) below), a TWE of covariance matrix can be computed as
$\hat{\mathbf{R}}_{\textup{TWE}}  =  \hat \theta_{\textup{TWE}} \cdot \hat \M_{\textup{TWE}}$,
where $ \hat \theta_{\text{TWE}} = h( \hat \nu_{\text{TWE}})$.

§ ESTIMATING THE TAIL PARAMETER OF ES DISTRIBUTION

From (<ref>) we can induce the following relationship between the scale parameter $\eta _{\text{cov}} = p^{-1}\tr(\mathbf{R})$ of the covariance matrix and scale $\eta = p^{-1}\tr(\M)$ of the scatter matrix:

η_cov =  θη⇔θ= η_cov/η.
Note that  a natural estimate of $\eta_{\text{cov}}$ is $p^{-1}\tr(\mathbf{S})$, where  $\mathbf{S} = \frac 1 n \sum_{i=1}^n \x_i \x_i^\top$ denotes the sample covariance matrix (SCM). On the other hand,  if density generator $g$ is specified up to unknown tail parameter $\nu$, thus indexed by $g_{\nu}(\cdot)$, then $\theta$ in (<ref>) is a following function of the tail parameter $\nu$:

θ= ∫_0^∞ t f_r^2(t; ν) d t

where the pdf $f_{r^2}(\cdot)$, defined in (<ref>), is one-to-one with $g_\nu(\cdot)$. We do not need numerical integration in most practical cases as often closed-form expression for $h(\nu)$ can be derived.  Then, after solving the inverse mapping, $\nu = h^{-1}(\theta) = h^{-1}(\eta_{\text{cov}}/\eta)$, Algorithm <ref> offers a unified approach for estimating the tail parameter of an ES distribution:

Distribution tail parameter learning
Input: Data  $\{\x_i\}_{i=1}^n$

Output: Estimated tail parameter $\nu$

1. Compute Tyler's M-estimator  $\hat \M$ and weights $\hat w_i$'s in (<ref>);

2. Compute $\hat \eta_{\text{TWE}}$ in (<ref>), and set  $\hat \theta_{\text{TWE}} = p^{-1}\tr(\S)/\hat \eta_{\text{TWE}}$;

3. Using (<ref>), estimate $\nu$ as $ \hat \nu_{\text{TWE}} = h^{-1}(\hat \theta_{\text{TWE}})$.

As an example, if $\x$ follows an MVT distribution with $\nu>2$ d.o.f., one has that  $\theta = h(\nu)=  \frac{\nu}{\nu -2}$, which
unfolds the relation:

ν= h^-1(θ) = 2 θ/ θ-1   θ>1.
Note that $\nu >2$ is required for the covariance matrix $\mathbf{R}$ to exist. The obtained estimator $ \hat \nu_{\text{TWE}}$  is closely related to estimator in <cit.>, referred to as OPP estimator for short.  OPP is an iterative approach that iteratively (re-)computes the maximum likelihood estimator (MLE)  $\hat \M$ of the MVT distribution with $\nu$ given by the current estimate of d.o.f. parameter $\nu^{(k)}$. It then computes $ \hat \theta = \tr(\S)/\tr(\hat \M)$ which provides an update $\nu^{(k+1)}=h^{-1}(\hat \theta)$ via (<ref>). The algorithm iterates for $k=0,1,2 \ldots$ until convergence,
starting from an initial start $\nu^{(0)} = \hat \nu_{\text{kurt}}$, where  $\hat \nu_{\text{kurt}}$  is an estimate of $\nu$ based on elliptical kurtosis,  proposed in  [18], and referred to as kurtosis estimator. We also proposed an improved version of OPP estimator in [19], which, however, is impractical for large $n$ and $p$. ML estimation of $\nu$ via the Expectation-Maximization (EM) approach  is considered in  [20]. This method is unfortunately rather unstable [21].

§ SIMULATION STUDIES

Boxplots d.o.f. estimates $\hat \nu$ as a function of $n$ when $\nu=5$ (top) and $\nu=3$ (bottom);  $p=100$, $\rho=0.6$.

MSE of d.o.f estimates $\hat \nu$ as a function of $n$ when $\nu=5$ (left panel) and $\nu=3$ (right panel);  $p=100$, $\rho=0.6$.

We first consider the case where scatter matrix
$\M$ has an autoregressive model (AR(1)) structure,
$(\boldsymbol{\Sigma})_{ij} = \eta \varrho^{|i-j|}$,
where $\eta = \tr(\M)/p$ is the scale parameter  and $\varrho$ is the correlation parameter, $\varrho \in (-1,1)$.   Since Tyler's M-estimator $\hat \M$ is invariant to the data scaling, we can set $\eta=1$ without favoring any estimator over the other.  The number of  Monte-Carlo runs is 5,000, and samples are generated from an MVT distribution with different choices of d.o.f. parameter $\nu$.

First, we investigate how the TWE of d.o.f. parameter $\nu$ compares against OPP and kurtosis estimators. <ref> displays the boxplots in the case that $p=100$ and varying sample lengths when
$\nu=5$ or $\nu=3$.
As can be noted, the proposed TWE attains the best accuracy as well as the smallest variability. Moreover, for $\nu=5$, its median values are right on  the spot. The kurtosis estimator obviously performs poorly when $\nu=3$ since the 4th-order moment does not exist in this case. <ref> shows the average mean squared error (MSE),  $(\hat  \nu - \nu)^2$, which further illustrates the benefits and high accuracy of the proposed TWE against its competitors.

Estimation errors  of scale estimates $\hat \eta$ as a function of $n$ when $\nu=5$ (top) and $\nu=3$ (bottom);  $p=100$, $\rho=0.6$.

<ref> displays the boxplots of different estimates of scale $\eta$. Here we compare $\hat \eta_{\text{TWE}}$ to OPP estimate of scale, defined as
$\hat \eta_{\text{OPP}} = \tr(\hat \M)/p $, where $\hat \M$ is the MLE of scatter based on $\nu = \hat \nu_{\text{OPP}}$.   We also compare with the  scale estimate provided by the SCM, defined as $\hat \eta_{\cov} = \tr(\S)/p$, but multiplied by $\theta^{-1}=(\nu-2)/\nu$ to obtain an estimate of $\eta$; recall (<ref>).
We can notice from <ref>  that TWE slightly underestimates the true scale while  OPP is overestimating. We also notice that the SCM estimator is clearly unbiased for $\nu=5$, but has huge variability. <ref> displays the median values of $\hat \nu$
for range of $\nu$ values when $n=150$ (and $p=100$ and $\varrho=0.6$ as earlier). The proposed TWE estimator
significantly outperforms the other estimators for all d.o.f. $\nu \in [3,8]$.

Estimated values of $\nu$; $p=100$, $\rho=0.6$, $n=150$.  The black line indicates the true value.

Estimated values of $\nu$ for synthetic stock market data; $p=100, \nu=4$.

We now consider an example based on stock market data. We generate synthetic data ($p=100$ assets) with heavy tails following MVT distribution with d.o.f. $\nu=4$ and covariance matrix as measured from stocks of S&P 500 index. Figure <ref> compares the estimated value of $\nu$ versus the number of observations for the following methods: kurtosis estimator, OPP estimator [18], and the proposed estimator, with the latter being clearly superior, illustrating its promising performance for real-world financial data.

§ CONCLUDING REMARKS

We  proposed a new robust estimator of scale parameter of an elliptical distribution based on the weights from Tyler's M-estimator, which was further used to construct an affine equivariant Tyler's M-estimator.
We then proposed a unified framework to estimate the tail parameter of an elliptical distribution.
Finally, it should be noted that this method generalizes to complex-valued data in a straightforward manner.

[1]
K.-T. Fang, S. Kotz, and K.-W. Ng, Symmetric Multivariate and Related
Distributions.1em plus 0.5em minus 0.4emLondon: Chapman and
hall, 1990.

[2]
E. Ollila, D. E. Tyler, V. Koivunen, and H. V. Poor, “Complex elliptically
symmetric distributions: survey, new results and applications,” IEEE
Transactions on Signal Processing, vol. 60, no. 11, pp. 5597–5625, 2012.

[3]
D. Paindaveine, “A canonical definition of shape,” Statistics &
probability letters, vol. 78, no. 14, pp. 2240–2247, 2008.

[4]
D. E. Tyler, “A distribution-free M-estimator of multivariate scatter,”
Annals of Statistics, vol. 15, no. 1, pp. 234–251, 1987.

[5]
Y. Sun, P. Babu, and D. P. Palomar, “Regularized Tyler's scatter estimator:
Existence, uniqueness, and algorithms,” IEEE Transactions on Signal
Processing, vol. 62, no. 19, pp. 5143–5156, 2014.

[6]
F. Pascal, Y. Chitour, and Y. Quek, “Generalized robust shrinkage estimator
and its application to STAP detection problem,” IEEE Transactions
on Signal Processing, vol. 62, no. 21, pp. 5640–5651, 2014.

[7]
E. Ollila and D. E. Tyler, “Regularized $M$-estimators of scatter matrix,”
IEEE Transactions on Signal Processing, vol. 62, no. 22, pp.
6059–6070, 2014.

[8]
Y. Chen, A. Wiesel, and A. O. Hero, “Robust shrinkage estimation of
high-dimensional covariance matrices,” IEEE Transactions on Signal
Processing, vol. 59, no. 9, pp. 4097–4107, 2011.

[9]
R. Couillet and M. McKay, “Large dimensional analysis and optimization of
robust shrinkage covariance matrix estimators,” Journal of
Multivariate Analysis, vol. 131, pp. 99–120, 2014.

[10]
I. Soloveychik and A. Wiesel, “Performance analysis of Tyler's covariance
estimator,” IEEE Transactions on Signal Processing, vol. 63, no. 2,
pp. 418–426, 2015.

[11]
T. Zhang, X. Cheng, and A. Singer, “Marčenko–pastur law for Tyler's
M-estimator,” Journal of Multivariate Analysis, vol. 149, pp.
114–123, 2016.

[12]
J. Goes, G. Lerman, and B. Nadler, “Robust sparse covariance estimation by
thresholding Tyler's M-estimator,” The Annals of Statistics,
vol. 48, no. 1, pp. 86–110, 2020.

[13]
E. Romanov, G. Kur, and B. Nadler, “Tyler's and Maronna's M-estimators:
Non-asymptotic concentration results,” Journal of Multivariate
Analysis, vol. 196, p. 105184, 2023.

[14]
D. P. Palomar, R. Zhou, X. Wang, F. Pascal, and E. Ollila, fitHeavyTail:
Mean and Covariance Matrix Estimation under Heavy Tails, 2023, r package
version 0.2.0. [Online]. Available:

[15]
R. A. Maronna, “Robust M-estimators of multivariate location and scatter,”
Annals of Statistics, vol. 5, no. 1, pp. 51–67, 1976.

[16]
P. L. Davies, “Asymptotic behaviour of S-estimates of multivariate location
parameters and dispersion matrices,” Annals of Statistics, pp.
1269–1292, 1987.

[17]
K. S. Tatsuoka and D. E. Tyler, “On the uniqueness of S-functionals and
M-functionals under nonelliptical distributions,” Annals of
Statistics, pp. 1219–1243, 2000.

[18]
E. Ollila, D. P. Palomar, and F. Pascal, “Shrinking the eigenvalues of
M-estimators of covariance matrix,” IEEE Transactions on Signal
Processing, vol. 69, pp. 256–269, 2021.

[19]
F. Pascal, E. Ollila, and D. P. Palomar, “Improved estimation of the degree of
freedom parameter of multivariate $ t $-distribution,” in 2021 29th
European Signal Processing Conference (EUSIPCO).1em plus 0.5em minus
0.4emIEEE, 2021, pp. 860–864.

[20]
C. Liu and D. B. Rubin, “ML estimation of the t-distribution using EM and
its extensions, ECM and ECME,” Statistica Sinica, vol. 5, no. 1,
pp. 19–39, 1995.

[21]
C. Fernandez and M. F. J. Steel, “Multivariate Student-t regression models:
Pitfalls and inference,” Biometrika, vol. 86, no. 1, pp. 153–167,


Diffusion models for Handwriting Generation

Troy Luhman***Equal contribution Eric Luhman11footnotemark: 1
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Abstract

  * In this paper, we propose a diffusion probabilistic model for handwriting generation. Diffusion models are a class of generative models where samples start from Gaussian noise and are gradually denoised to produce output. Our method of handwriting generation does not require using any text-recognition based, writer-style based, or adversarial loss functions, nor does it require training of auxiliary networks. Our model is able to incorporate writer stylistic features directly from image data, eliminating the need for user interaction during sampling. Experiments reveal that our model is able to generate realistic , high quality images of handwritten text in a similar style to a given writer. Our implementation can be found at https://github.com/tcl9876/Diffusion-Handwriting-Generation.

##  1 Introduction

Deep generative models have been able to produce realistic handwritten text.
Handwriting data can be stored in online or offline format. Online data is
rendered as a sequence of pen strokes, and offline data is stored directly as
an image. Online temporal data is easier to work with due to its lower
dimensionality, but is harder to collect since each stroke must be recorded
during writing. Each sample of handwriting data is associated with a text
sequence label describing the content of the handwritten text. Each
handwritten sample can also be associated with another example of handwriting
from the same writer, which provides information about the style of the
writer.

The task of handwriting generation was first done by [Graves, 2013], who used
a RNN to synthesize online handwriting. Since then, most work on online
handwriting generation has used an RNN based architecture, while most work on
offline handwriting generation has used a GAN. Both methods have been able to
generate realistic, diverse samples in a given writer’s handwriting style. RNN
based models [Chung et al., 2015, Aksan et al., 2018, Kotani et al., 2020] are
simple to train and sample from, but one challenge they face is that they
require the writer style data to be in online format. This poses a challenge
for incorporating stylistic information, as the user’s pen strokes must be
recorded during sampling. GAN based methods for handwriting generation [Alonso
et al., 2019, Kang et al., 2020, Fogel et al., 2020] have been used as an
alternative to online autoregressive models. GAN based methods have the
advantage of being able to incorporate offline stylistic information, as they
work with images directly. The downside to GAN-based methods is that they can
be difficult to train and may suffer from lack of diversity due to mode
collapse. Additionally, previous work on GAN-based handwriting generation
requires training at least two auxiliary networks: the discriminator and a
text recognition network to control the content of the generated text.

Diffusion probabilistic models [Sohl-Dickstein et al., 2015, Ho et al., 2020]
use a Markov chain to convert a known distribution (e.g. Gaussian) into a more
complex data distribution. A diffusion process converts the data distribution
into a simple distribution by iteratively adding Gaussian noise to the data,
and the generative model learns to reverse this diffusion process. Diffusion
models can be trained to optimize a weighted variational lower bound of the
data likelihood. This objective is similar to that of noise conditional score
networks [Song and Ermon, 2019], which estimates the gradients of the data
distribution. Both noise conditional score networks and diffusion models
generate samples by starting from Gaussian noise and gradually removing this
noise.

In this paper, we propose a diffusion probabilistic model for online
handwriting generation. Our proposed model has advantages over both
autoregressive and GAN based methods of handwriting generation. Our model,
though it generates samples in online format, is able to incorporate writer
stylistic features from offline data, eliminating the need for user
interaction during sampling. Our method of handwriting generation does not
require using any text based, style based, or adversarial loss functions, nor
does it require training of auxiliary networks. As a result, our model has a
very simple training and sampling procedure. Our model is able to generate
realistic samples of handwritten text in a similar style to the original
writer.

##  2 Diffusion Models for Handwriting Generation

###  2.1 Diffusion Probabilistic Models

Let $q(\text{y}_{0})$ be the data distribution, and let $\text{y}_{1}$, …,
$\text{y}_{T}$ be a series of $T$ latent variables with the same
dimensionality as $\text{y}_{0}$. A diffusion model consists of two processes:
a noise adding diffusion process, and a reverse process [Sohl-Dickstein et
al., 2015]. The posterior $q\left(\text{y}_{1:T}\ \mid\text{y}_{0}\right)$, or
the diffusion process, is defined as a fixed Markov chain where Gaussian noise
is added at each iteration based on a fixed noise schedule
$\beta_{1},\ldots,\beta_{T}$:

$\displaystyle q\left(\text{y}_{1:T}\ \mid\text{y}_{0}\right)=\
\prod_{t=1}^{T}q\left(\text{y}_{t\ }\mid\text{y}_{t-1}\right),\hskip
20.0ptq\left(\text{y}_{t}\ \mid\text{y}_{t-1}\right)=\
\mathcal{N}\left(\text{y}_{t};\sqrt{1-\beta_{t}}\text{y}_{t-1},\beta_{t}\bm{I}\right)$
(1)

The reverse process is defined as a Markov chain parameterized by $\theta$:

$\displaystyle p(\text{y}_{T})=\mathcal{N}(\text{y}_{T};0,\bm{I}),\hskip
20.0ptp_{\theta}\left(\text{y}_{0:T}\ \right)=p\left(\text{y}_{T}\right)\
\prod_{t=1}^{T}p_{\theta}\left(\text{y}_{t\ }\mid\text{y}_{t-1}\right)$ (2)

where $p_{\theta}\left(\text{y}_{t-1}\ \mid\text{y}_{t}\right)$ intends to
reverse the effect of the noise adding process $q\left(\text{y}_{t}\
\mid\text{y}_{t-1}\right)$:

$\displaystyle p_{\theta}\left(\text{y}_{t-1}\ \mid\text{y}_{t}\right)=\
\mathcal{N}\left(\text{y}_{t-1};\mu_{\theta}\left(\text{y}_{t},t\right),\
\Sigma_{\theta}\left(\text{y}_{t},t\right)\right)$ (3)

where $\Sigma_{\theta}\left(\text{y}_{t},t\right)=\sigma_{t}^{2}\bm{I}$ and
$\sigma_{t}^{2}$ is a constant related to $\beta_{t}$. The forward process
posterior is defined as:

$\displaystyle
q(\text{y}_{t-1}\mid\text{y}_{t},\text{y}_{0})=\mathcal{N}(\text{y}_{t-1};\tilde{\mu}\left(\text{y}_{t},\text{y}_{0}\right),\sigma_{t}^{2}\bm{I})$
(4)

[Ho et al., 2020] showed that the ELBO can be calculated in closed form by
expanding it into a series of KL Divergences between Gaussian distributions:

$\displaystyle\textbf{ELBO}=-\mathbb{E}_{q}\left(D_{KL}\left(q(\text{y}_{T}\mid\text{y}_{0})\
\middle\|\
p(\text{y}_{T})\right)+\sum_{t=2}^{T}D_{KL}\left(q(\text{y}_{t-1}\mid\text{y}_{t},\text{y}_{0})\
\middle\|\ p_{\theta}(\text{y}_{t-1}\mid\text{y}_{t})\right)-\log
p_{\theta}\left(\text{y}_{0}\mid\text{y}_{1}\right)\right)$ (5)

By defining some constants, $\text{y}_{t}$ can be calculated in closed form
for any step t:

$\displaystyle\alpha_{t}=1-\beta_{t},\
\bar{\alpha_{t}}=\prod_{s=1}^{t}\alpha_{s},\
\epsilon\thicksim\mathcal{N}(0,\bm{I}),\
\text{y}_{t}=\sqrt{\bar{\alpha_{t}}}\text{y}_{0}+\sqrt{1-\bar{\alpha_{t}}}\epsilon$
(6)

The loss function then becomes the following for some timestep $t-1$, where
$\mu_{\theta}$ is a model that predicts the forward process posterior mean
$\tilde{\mu}$:

$\displaystyle
L_{t-1}=E_{\text{y}_{0},\epsilon}\left(\frac{1}{2\sigma^{2}}\left\|\frac{1}{\sqrt{\alpha_{t}}}\left(\text{y}_{t}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon\right)-\mu_{\theta}\left(\text{y}_{t},t\right)\right\|^{2}\right)$
(7)

[Ho et al., 2020] observed that diffusion probabilistic models can be
reparameterized to resemble score-based generative models [Song and Ermon,
2019]. Score-based generative models estimate the gradient of the logarithmic
data density $\nabla_{x}\log\left(p\left(x\right)\right)$ with a denoising
objective. Under this reparameterization, the diffusion model,
$\epsilon_{\theta}$, predicts $\epsilon$ similar to the objective of score-
based models. The loss function for predicting $\epsilon$ becomes:

$\displaystyle
L_{t-1}=\mathbb{E}_{t,\epsilon}\left[C_{t}\|\epsilon-\epsilon_{\theta}(\text{y}_{t},t)\|_{2}^{2}\right],C_{t}=\frac{\beta^{2}}{2\sigma_{t}^{2}\alpha_{t}\left(1-\bar{\alpha_{t}}\right)}$
(8)

During sampling, diffusion probablilistic models iteratively remove the noise
added in the diffusion process, by sampling $\text{y}_{t-1}$ for
$t=T,\ldots,1$:

$\displaystyle\text{y}_{t-1}=\frac{1}{\sqrt{a_{t}}}\left(\text{y}_{t}-\frac{\beta_{t}}{\sqrt{1-\bar{\alpha_{t}}}}\epsilon_{\theta}\left(\text{y}_{t},t\right)\right)+\sigma_{t}z$
(9)

where $z\thicksim\mathcal{N}(0,\bm{I})$ and $\sigma_{t}$ is a constant related
to $\beta_{t}$. For our experiments, we used $\sigma_{t}^{2}=\beta_{t}$.

In our experiments, we found it beneficial to make a modification to the
original sampling procedure. First, as our model predicts $\epsilon$, we can
use the estimate of $\epsilon$ to approximate $\text{y}_{0}$ at any timestep
$t$:

$\displaystyle\text{y}_{0}\approx\hat{\text{y}}_{0}=\frac{1}{\sqrt{\bar{\alpha_{t}}}}\left(\text{y}_{t}-\sqrt{1-\bar{\alpha_{t}}}\epsilon_{\theta}(\text{y}_{t},t)\right)$
(10)

We now use our approximation for $\text{y}_{0}$ to find $\text{y}_{t-1}$:

$\displaystyle\text{y}_{t-1}=\sqrt{\bar{\alpha}_{t-1}}\hat{\text{y}}_{0}+\sqrt{1-\bar{\alpha}_{t-1}}z,\
z\thicksim\mathcal{N}(0,\bm{I})$ (11)

We first estimate $\hat{\text{y}}_{0}$ from $\text{y}_{t}\text{and}\
\epsilon_{\theta}$ using Equation 10, then find $\text{y}_{t-1}$ according to
Equation 11. Combining these two, we get:

$\displaystyle\text{y}_{t-1}=\frac{1}{\sqrt{\alpha_{t}}}\left(\text{y}_{t}-\sqrt{1-\bar{\alpha_{t}}}\epsilon_{\theta}(\text{y}_{t},t)\right)+\sqrt{1-\bar{\alpha}_{t-1}}z,\
z\thicksim\mathcal{N}(0,\bm{I})$ (12)

The training procedure is performed by minimizing Equation 8 with the $C_{t}$
term removed. We compared sampling $\text{y}_{t-1}$ according to Equation 12
as opposed to Equation 9 and found that the Equation 12 led to more realistic
samples (see ablation study in Table 1), at the cost of a slight decrease in
diversity.

###  2.2 Conditional Handwriting Generation

Figure 1: Samples from the IAM Database.

We will now address the task of conditional online handwriting generation.
Each data point $\text{x}_{0}$ is composed of a sequence of $N$ vectors
$x_{1}\ldots x_{N}$. Each individual vector in the sequence
$x_{n}\in\mathbb{R}^{2}\times\\{0,1\\}$ is composed of a real valued pair
which represents the pen offset from the previous stroke in the x and y
direction, and a binary entry that has a value of 0 if the pen was down when
writing the stroke and 1 otherwise. Each handwritten sequence is associated
with a discrete character sequence c describing what was written. Each
sequence is also associated with an offline image containing writer’s style
information, denoted by s.

There is one technical issue that needs to be addressed. The reverse process
in Equation 3 is parameterized by a Gaussian distribution. We cannot
parameterize the binary variable representing whether the stroke was drawn by
a Gaussian distribution as we did for the real valued pen strokes. However, we
can instead parameterize it with a Bernoulli distribution, which can also be
optimized in closed form.

We therefore split each data point $\text{x}_{0}$ into two sequences
$\text{y}_{0}$ and $\text{d}_{0}$ of equal length, with $\text{y}_{0}$
representing the real valued pen strokes, and $\text{d}_{0}$ representing
whether the stroke was drawn. At each step t, our model
$\text{d}_{\theta}(\text{y}_{t},\text{c},\text{s},\sqrt{\bar{\alpha}})$
returns an estimate $\hat{\text{d}_{0}}$ of whether the pen was down.
$\text{d}_{\theta}$ shares all parameters with $\epsilon_{\theta}$. Our stroke
loss and pen-draw loss are shown in Equations 13 and 14:

$\displaystyle
L_{\text{stroke}}(\theta)=\|\epsilon-\epsilon_{\theta}\left(\text{y}_{t},\text{c},\text{s},\sqrt{\bar{\alpha}}\right)\|_{2}^{2}$
(13) $\displaystyle
L_{\text{drawn}}(\theta)=-\text{d}_{0}\log\left(\hat{\text{d}_{0}}\right)-\left(1-\text{d}_{0}\right)\log\left(1-\hat{\text{d}_{0}}\right)$
(14)

We found it beneficial to weight the pen-draw loss according to the noise
level, since it is more difficult to predict the pen-draws at higher noise
levels. We weight Equation 14 by $\bar{\alpha}$ during training (see Algorithm
1).

[Song and Ermon, 2019, 2020] noted that the choice of noise schedule is
crucial to generating high quality samples. [Chen et al., 2020] proposed to
condition the model on the continuous noise level $\sqrt{\bar{\alpha}}$ as
opposed to the discrete index $t$. This allows for the use of different noise
schedules during sampling without retraining the model. In order to condition
the network on the continuous noise level, we first define noise schedule l
where $l_{0}=1,l_{t}=\sqrt{\bar{\alpha_{t}}}$ During training, we condition
the model on $\sqrt{\bar{\alpha}}\thicksim\text{Uniform}(l_{t-1},l_{t})$,
where $t\thicksim\text{Uniform}(\\{1,\ldots,T\\})$.

One desirable property of handwriting synthesis networks is the ability to
generate handwriting in the style of a given writer. Previous handwriting
synthesis methods are able to control the style of the generated handwriting
by conditioning on a sample written by the writer. Although our model
generates online output samples, our model accepts offline images as input for
writer stylistic features. Previous RNN-based models require online data
during sampling, which is difficult to collect. To incorporate style
information, we extract features from the image using Mobilenet [Sandler et
al., 2018] and use these features as input to our model.

We describe our training procedure in Algorithm 1 and our sampling procedure
in Algorithm 2. Note that our training procedure does not require training any
auxiliary networks and does not use text recognition or style-based losses.

1

2while not converged do

3 $\text{y}_{0}\thicksim q(\text{y}_{0})$

4 $t\thicksim\text{Uniform}(\\{1,\ldots,T\\})$

5 $\sqrt{\bar{\alpha}}\thicksim\text{Uniform}(l_{t-1},l_{t})$

6 $\epsilon\thicksim\mathcal{N}(0,\bm{I})$

7 $\text{y}_{t}=\sqrt{\bar{\alpha}}\text{y}_{0}+\sqrt{1-\bar{\alpha}}\epsilon$

8 Take gradient descent step on:

9
$\nabla_{\theta}(L_{\text{stroke}}(\theta)+\bar{\alpha}L_{\text{drawn}}(\theta))$

10 end while

11

Algorithm 1 Training

1 $\text{y}_{T}\thicksim\mathcal{N}(0,\bm{I})$

2 for $t=T,\ldots,1$ do

3 $z\thicksim\mathcal{N}(0,\bm{I})$

4
$\text{y}_{t-1}=\frac{1}{\sqrt{\alpha_{t}}}\left(\text{y}_{t}-\sqrt{1-\bar{\alpha_{t}}}\epsilon_{\theta}(\text{y}_{t},\text{c},\text{s},\sqrt{\bar{\alpha_{t}}})\right)$

5 $\text{y}_{t-1}=\text{y}_{t-1}+\sqrt{1-\bar{\alpha}_{t-1}}z$ if $t>1$

6
$\text{d}_{0}=\text{d}_{\theta}(\text{y}_{t},\text{c},\text{s},\sqrt{\bar{\alpha_{t}}})$

7 end for

8return $\text{y}_{0},\text{d}_{0}$

Algorithm 2 Sampling

##  3 Model Architecture

Figure 2: Overview of our full model.

Our full model is shown in Figure 2. It consists of two parts: a text and
style encoder to represent the desired text and stylistic features, and a
diffusion probabilistic model to predict $\epsilon_{\theta}$. The noise level
$\sqrt{\bar{\alpha}}$ is passed through a feedforward network consisting of
two fully connected layers. Both the encoder and the diffusion model are
conditioned on the noise level.

Conditioning on the noise level: To condition our model on the noise level, we
use affine transformations along the channel axis. Each affine
transformation’s scales and biases are parameterized by the output of a fully
connected layer. Our affine transformation is detailed in Figure 6 in Appendix
6.

Text and Style Conditioning: For writer conditioning, we first extract local
features from an image of their handwriting using a MobileNetV2 [Sandler et
al., 2018] pretrained on Imagenet. Character-level embeddings are used to
represent the text sequence. We then compute attention between the text
sequence and the extracted features, allowing different text characters to
attend to different portions of the given writer sample. This output is then
added to the text sequence representation, before being passed through a
feedforward network. Full architectural details on the encoder can be seen in
Figure 10 in Appendix 10.

Diffusion model: Our diffusion model consists of downsampling blocks, followed
by upsampling blocks, and uses long range convolutional skip connections. We
use two main types of blocks, convolutional blocks and attentional blocks.
Full information on our diffusion model architecture is detailed in Figure 9
in Appendix 9.

Our convolutional blocks consist of 3 convolutional layers and a convolutional
skip connection. We apply conditional affine transformations to the output of
every convolutional layer. Further architectural information on our
convolutional block is shown in Figure 7 in Appendix 7.

Our attentional blocks consist of 2 multi-head attention layers, and a feed
forward network. The first attention layer performs attention between the
stroke sequences x and the output of the text-style encoder, while the second
performs self-attention. Sinusoidal positional encodings [Vaswani et al.,
2017] are added to the queries and keys at every attentional layer. For
attentional layers at higher stroke resolutions, the positions of the text
sequence are multiplied before the positional encoding to allow for easier
alignment between the text sequence and the longer stroke sequence. We use
layer normalization [Ba et al., 2016] followed by conditional affine
transformations after every attentional layer and feed forward network.
Further architectural information on our attentional block is shown in Figure
8 in appendix 8.

##  4 Related Work

Our work builds off previous work on diffusion probabilistic models [Sohl-
Dickstein et al., 2015, Ho et al., 2020] and Noise Conditional Score Networks
or NSCNs [Vincent, 2011, Song and Ermon, 2019], which are closely related. [Ho
et al., 2020] addressed the task of image generation using a diffusion model,
and [Kong et al., 2020, Chen et al., 2020] used diffusion models for speech
generation. [Song and Ermon, 2019, 2020] used NSCNs for image generation, and
[Cai et al., 2020] used NCSNs for shape generation. [Jolicoeur-Martineau et
al., 2020] combined the NCSN objective with an adversarial objective to
address image generation.

Handwriting synthesis was first explored in [Graves, 2013], which used a
recurrent neural network (RNN) to predict each stroke one at a time. Since
then, most work on online handwriting generation has made use of a type of RNN
knows as a Variational RNN [Chung et al., 2015]. Variational RNNs (VRNNs)
incorporate a Variational Autoencoder [Kingma and Welling, 2014], or VAE , at
each timestep to generate sequences. [Aksan et al., 2018] uses a conditional
VRNN that is able to separate content and style elements, allowing for editing
of the generated samples. [Kotani et al., 2020] incorporates character and
writer level style information for further model flexibity. RNN based
architectures exhibit high quality, diverse samples and have many
capabilities. These drawback of autoregressive models is that they accept
online writing examples in order to incorporate the writer’s stylistic
information, requiring the writer to interact during sampling. [Mayr et al.,
2020] alleviates this problem by approximating an online representation of an
offline image with the downside being that it is highly dependent on the
performance of the approximation model.

As an alternative to autoregressive RNN models, generative adversarial
networks (GANs) [Goodfellow et al., 2014] have been proposed for handwriting
synthesize. [Alonso et al., 2019] first used a GAN to generate images of
handwritten words. [Kang et al., 2020, Fogel et al., 2020] both use the GAN
framework and improve sample quality over [Alonso et al., 2019], in addition
to being able to condition on writer stylistic features. GAN based methods of
handwriting generation are not autoregressive, and are able to condition on
offline style features, which is beneficial since no user interaction is
required during sampling. The downside of GANs is that adversarial training
procedure of GANs can be unstable and result in mode collapse, which reduces
the diversity of samples produced. Previous works on handwritten text
generation with GANs also have used a text recognition objective for the
generated samples which can complicate the training procedure.

##  5 Experiments

Figure 3: Four of the above samples are real, drawn from the test dataset. The
remaining 8 were randomly generated by our model. We encourage the reader to
try to determine which ones are fake, and which are real. Answer in Appendix B

###  5.1 Experiment Details

We use the IAM Online Database [Liwicki and Bunke, 2005], a dataset of
approximately 12000 handwritten lines that are each associated with a
character string label. Each sample is also associated with an sample of
handwritten text by the writer, which we store as an offline image. We use the
same splits as [Graves, 2013]. To preprocess the data, we first divide each
example by its standard deviation. We discard samples containing strokes of
lengths 15 standard deviations above the mean. These very long strokes usually
are a result of large spacing between words. We combine strokes pointing in
approximately the same direction as their neighbor, which reduces the
dimensionality without visual differences.

We used $T=60$ diffusion steps for our models. For our choice of noise
schedule $\beta_{1},\ldots,\beta_{T}$, we used $\beta_{t}$ = 0.02 +
Exponential($1\times 10^{-5},0.4$). Exponential($1\times 10^{-5},0.4$) denotes
a geometric sequence from $1\times 10^{-5}$ to 0.4 . We trained our model for
60000 steps with a batch size of 96 on a single Nvidia V-100. We use the Adam
optimizer [Kingma and Ba, 2017] with $\beta_{1}=0.9$ , $\beta_{2}=0.98$, and
we clip the norm of the gradients to 100. For our learning rate, we use the
inverse square root schedule in [Vaswani et al., 2017] with 10000 warmup steps
and the $d_{model}$ argument at 256. All our models have the same size of 10.0
million parameters. These values were not sweeped over.

We evaluate the performance of our models with Frechet Inception Distance
(FID) [Heusel et al., 2017] and Geometry Score (GS) [Khrulkov and Oseledets,
2018]. FID uses the InceptionV3 network [Szegedy et al., 2016] to compare the
similarity of real images to generated images. Since our model outputs online
sequences as output, we first plot the sequences and then convert them to
offline images. The Geometry score compares the geometric properties of the
real and generated data manifolds to measure quality and diversity.

To conduct the experiment, we draw two sets of 4500 examples from the training
dataset. We compare the first set of examples with our model’s predictions of
the second set of examples given the corresponding text and style information.
For the FID evaluation, we resize all the plotted images to size M $\times$ N.
For the Geometry Score evaluation, we keep the generated data in online form,
and use the same values of $L_{0},\gamma,i_{max}$ as [Khrulkov and Oseledets,
2018]. We run both experiments once.

###  5.2 Results and Discussion

We report our scores for the objective evaluation metrics in Table 1.

Table 1: Objective evaluation metrics (lower is better) | FID | GS
---|---|---
Ground Truth | 2.91 | $5.4\times 10^{-4}$
Our Model | 7.10 | $3.3\times 10^{-3}$
Ablated Model | 8.05 | $2.7\times 10^{-3}$

Due to differences in the dataset and format of the generated output, we do
not directly compare our results with previous works. Our images for FID are
first plotted on a graph, then the graph is converted to an image. This is in
contrast to previous methods that used FID, which generated offline images
directly. Our GS is computed from online data, which cannot be compared with
scores from offline data. We instead compare our results with online samples
from the dataset itself. The ground truth score is calculated with the real
examples in the second set as opposed to the model’s predictions of examples
in the second set.

Ablation Study We conduct an ablation study to compare the quality of samples
generated with our sampling procedure (in Algorithm 2) to those generated with
the original sampling procedure, which uses Equation 9. We compare FID and GS
for the generated samples in Table 1. We used the same trained model for both
models, and only changed the sampling procedure. We found that in the case of
online handwriting generation, our sampling procedure yielded relatively
higher FID scores but relatively lower Geometry Score. For qualitative
comparison between the two, see Figures 15 and 16 in the Appendix.

Attention Weights Learning to properly align the text sequence to the stroke
sequence is critical for generating longer sequences of realistic text without
missing, misplaced or confused letters. The diagonal line in Figure 4
represents the predicted alignment between the strokes and the text sequence.
Empirically, when the attention weights diverged from the diagonal line, this
led to spelling mistakes. We observed that the text-to-stroke attention
weights remain well-aligned throughout the reverse process, even in the
earliest stages where the text itself is not recognizable.

Figure 4: Text - to - stroke attention weights at $t=5,t=30,t=55$. The
diagonal line shows the alignment predicted by the model between the strokes
and the texts.

Style Interpolation Given two samples of handwriting from two different
writers, we would like to be able to interpolate a sample that shares
stylistic aspects with both writers. In order to incorporate style
information, our model accepts a feature vector s that is obtained by extract
features from the given image using MobileNet. To perform style interpolation,
we first obtain feature vectors $\text{s}_{0}^{\prime},\
\text{s}_{1}^{\prime}$ from $\text{s}_{0},\ \text{s}_{1}$ and condition our
model on $\hat{s}$, where
$\hat{s}=\lambda\text{s}_{0}^{\prime}+(1-\lambda)\text{s}_{1}^{\prime}$.
$\lambda\in[0,1]$ controls the relative importance of $\text{s}_{0}^{\prime}$.

Figure 5: Style Interpolation. The top and bottom images $\text{s}_{0},\
\text{s}_{1}$ contain writer style information and are drawn from the test
dataset. To perform style interpolation, we keep the generated text constant,
and gradually decrease $\lambda$. The generated samples share calligraphic
aspects with the bottom picture $\text{s}_{1}$ more and less with the top
picture $\text{s}_{0}$ less as $\lambda$ decreases.

##  6 Conclusion

In this paper, we have proposed a diffusion probabilistic model for online
handwriting generation. Our work builds off of previous work in both diffusion
probabilistic models [Ho et al., 2020], and closely related score based
generative models [Song and Ermon, 2019]. Our training procedure is more
simple than most previous handwriting generation methods, as ours does not
require adversarial, text-recognition based or writer-recognition based
objectives. We also demonstrate the ability for our model to incorporate
calligraphic features from offline data, which is easier to collect during
sampling time. Empirically, the samples from our model are realistic and
diverse.

## References

  * Aksan et al. [2018] E. Aksan, F. Pece, and O. Hilliges.  Deepwriting: Making digital ink editable via deep generative modeling, 2018.
  * Alonso et al. [2019] E. Alonso, B. Moysset, and R. Messina.  Adversarial generation of handwritten text images conditioned on sequences, 2019.
  * Ba et al. [2016] J. Ba, J. Kiros, and G. Hinton.  Layer normalization.  07 2016.
  * Block et al. [2020] A. Block, Y. Mroueh, A. Rakhlin, and J. Ross.  Fast mixing of multi-scale langevin dynamics under the manifold hypothesis, 2020.
  * Cai et al. [2020] R. Cai, G. Yang, H. Averbuch-Elor, Z. Hao, S. Belongie, N. Snavely, and B. Hariharan.  Learning gradient fields for shape generation, 2020.
  * Chen et al. [2020] N. Chen, Y. Zhang, H. Zen, R. J. Weiss, M. Norouzi, and W. Chan.  Wavegrad: Estimating gradients for waveform generation, 2020.
  * Chung et al. [2015] J. Chung, K. Kastner, L. Dinh, K. Goel, A. Courville, and Y. Bengio.  A recurrent latent variable model for sequential data, 2015.
  * Fogel et al. [2020] S. Fogel, H. Averbuch-Elor, S. Cohen, S. Mazor, and R. Litman.  Scrabblegan: Semi-supervised varying length handwritten text generation, 2020.
  * Goodfellow et al. [2014] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial nets.  In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 27_ , pages 2672–2680. Curran Associates, Inc., 2014.  URL http://papers.nips.cc/paper/5423-generative-adversarial-nets.pdf.
  * Graves [2013] A. Graves.  Generating sequences with recurrent neural networks, 2013.
  * Heusel et al. [2017] M. Heusel, H. Ramsauer, T. Unterthiner, B. Nessler, and S. Hochreiter.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 6626–6637. Curran Associates, Inc., 2017.  URL http://papers.nips.cc/paper/7240-gans-trained-by-a-two-time-scale-update-rule-converge-to-a-local-nash-equilibrium.pdf.
  * Ho et al. [2020] J. Ho, A. Jain, and P. Abbeel.  Denoising diffusion probabilistic models, 2020.
  * Jolicoeur-Martineau et al. [2020] A. Jolicoeur-Martineau, R. Piché-Taillefer, R. T. des Combes, and I. Mitliagkas.  Adversarial score matching and improved sampling for image generation, 2020.
  * Kang et al. [2020] L. Kang, P. Riba, Y. Wang, M. Rusiñol, A. Fornés, and M. Villegas.  Ganwriting: Content-conditioned generation of styled handwritten word images, 2020.
  * Khrulkov and Oseledets [2018] V. Khrulkov and I. Oseledets.  Geometry score: A method for comparing generative adversarial networks, 2018.
  * Kingma and Ba [2017] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization, 2017.
  * Kingma and Welling [2014] D. P. Kingma and M. Welling.  Auto-encoding variational bayes, 2014.
  * Kong et al. [2020] Z. Kong, W. Ping, J. Huang, K. Zhao, and B. Catanzaro.  Diffwave: A versatile diffusion model for audio synthesis, 2020.
  * Kotani et al. [2020] A. Kotani, S. Tellex, and J. Tompkin.  Generating handwriting via decoupled style descriptors, 2020.
  * Liu et al. [2016] Q. Liu, J. D. Lee, and M. I. Jordan.  A kernelized stein discrepancy for goodness-of-fit tests and model evaluation, 2016.
  * Liwicki and Bunke [2005] M. Liwicki and H. Bunke.  Iam-ondb - an on-line english sentence database acquired from handwritten text on a whiteboard.  _Eighth International Conference on Document Analysis and Recognition (ICDAR’05)_ , pages 956–961 Vol. 2, 2005.
  * Mayr et al. [2020] M. Mayr, M. Stumpf, A. Nikolaou, M. Seuret, A. Maier, and V. Christlein.  Spatio-temporal handwriting imitation, 2020.
  * Pang et al. [2020] T. Pang, K. Xu, C. Li, Y. Song, S. Ermon, and J. Zhu.  Efficient learning of generative models via finite-difference score matching, 2020.
  * Sandler et al. [2018] M. Sandler, A. Howard, M. Zhu, A. Zhmoginov, and L.-C. Chen.  Mobilenetv2: Inverted residuals and linear bottlenecks.  pages 4510–4520, 06 2018.  doi: 10.1109/CVPR.2018.00474.
  * Sohl-Dickstein et al. [2015] J. Sohl-Dickstein, E. A. Weiss, N. Maheswaranathan, and S. Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics, 2015.
  * Song and Ermon [2019] Y. Song and S. Ermon.  Generative modeling by estimating gradients of the data distribution, 2019\.
  * Song and Ermon [2020] Y. Song and S. Ermon.  Improved techniques for training score-based generative models, 2020.
  * Szegedy et al. [2016] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna.  Rethinking the inception architecture for computer vision.  06 2016.  doi: 10.1109/CVPR.2016.308.
  * Vaswani et al. [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin.  Attention is all you need.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 5998–6008. Curran Associates, Inc., 2017.  URL http://papers.nips.cc/paper/7181-attention-is-all-you-need.pdf.
  * Vincent [2011] P. Vincent.  A connection between score matching and denoising autoencoders.  _Neural Computation_ , 23:1661–1674, 2011.

##  Appendix A Architecture Details

Figure 6: Our Conditional Affine Transformation. Figure 7: Our Convolutional
block.

Figure 8: Our Attentional block. Figure 9: Architecture of the diffusion
model. Figure 10: Architecture of the Encoder.

##  Appendix B Samples

The phrases ”as I ever saw”, ”that the attractive forces”, ”the traverser had
eaten”, and ”you will find a number of” in Figure 3 are real.

Figure 11: Samples generated from our model. On the left: ”the quick brown fox
jumps”. On the right: ”Samples from a diffusion model” Figure 12: Samples
generated from our model. Figure 13: Samples generated from our model. Figure
14: Samples generated from our model. Both of the texts, ”night bring only a
disturbed tossing” and ”here in the same words in” are from the test dataset.
Figure 15: In each column, the top line is a sample of a writer from the test
dataset’s style. The bottom two lines are produced by the generative model
using the style information s in the first line.  Figure 16: The left column
shows samples from the ablated model, which uses Equation 9 to compute
$\text{y}_{t-1}$. The right column shows samples from our final model.  Figure
17: The left column shows samples from the ablated model, which uses Equation
9 to compute $\text{y}_{t-1}$. The right column shows samples from our final
model.