Fast, Asymptotically Efficient, Recursive Estimation in a Riemannian Manifold

Abstract

Stochastic optimisation in Riemannian manifolds, especially the Riemannian stochastic gradient method, has attracted much recent attention. The present work applies stochastic optimisation to the task of recursive estimation of a statistical parameter which belongs to a Riemannian manifold. Roughly, this task amounts to stochastic minimisation of a statistical divergence function. The following problem is considered: how to obtain fast, asymptotically efficient, recursive estimates, using a Riemannian stochastic optimisation algorithm with decreasing step sizes. In solving this problem, several original results are introduced. First, without any convexity assumptions on the divergence function, we proved that, with an adequate choice of step sizes, the algorithm computes recursive estimates which achieve a fast non-asymptotic rate of convergence. Second, the asymptotic normality of these recursive estimates is proved by employing a novel linearisation technique. Third, it is proved that, when the Fisher information metric is used to guide the algorithm, these recursive estimates achieve an optimal asymptotic rate of convergence, in the sense that they become asymptotically efficient. These results, while relatively familiar in the Euclidean context, are here formulated and proved for the first time in the Riemannian context. In addition, they are illustrated with a numerical application to the recursive estimation of elliptically contoured distributions.

1. Introduction

Over the last five years, the data science community has devoted significant attention to stochastic optimisation in Riemannian manifolds. This was impulsed by Bonnabel, who proved the convergence of the Riemannian stochastic gradient method [1]. Later on [2], the rate of convergence of this method was studied in detail and under various convexity assumptions on the cost function. More recently, asymptotic efficiency of the averaged Riemannian stochastic gradient method was proved in [3]. Previously, for the specific problem of computing Riemannian means, several results on the convergence and asymptotic normality of Riemannian stochastic optimisation methods had been obtained [4,5]. The framework of stochastic optimisation in Riemannian manifolds is far-reaching, and encompasses applications to principal component analysis, dictionary learning, and tensor decomposition, to give only a few examples [6,7,8].

The present work moves in a different direction, focusing on recursive estimation in Riemannian manifolds. While recursive estimation is a special case of stochastic optimisation, it has its own geometric structure, given by the Fisher information metric. Here, several original results will be introduced, which show how this geometric structure can be exploited, to design Riemannian stochastic optimisation algorithms which compute fast, asymptotically efficient, recursive estimates, of a statistical parameter which belongs to a Riemannian manifold. For the first time in the literature, these results extend, from the Euclidean context to the Riemannian context, the classical results of [9,10].

The mathematical problem, considered in the present work, is formulated in Section 2. This involves a parameterised statistical model P of probability distributions $P_{\theta }$, where the statistical parameter $\theta$ belongs to a Riemannian manifold $\mathsf{\Theta }$. Given independent observations, with distribution $P_{{\theta }^{*}}$ for some ${\theta }^{*}\in \mathsf{\Theta }$, the aim is to estimate the unknown parameter ${\theta }^{*}$. In principle, this is done by minimising a statistical divergence function $D\left(\theta \right)$, which measures the dissimilarity between $P_{\theta }$ and $P_{{\theta }^{*}}$. Taking advantage of the observations, there are two approaches to minimising $D\left(\theta \right)$: stochastic minimisation, which leads to recursive estimation, and empirical minimisation, which leads to classical techniques, such as maximum-likelihood estimation [11,12].

The original results, obtained in the present work, are stated in Section 3. In particular, these are Propositions 2, 4, and 5. Overall, these propositions show that recursive estimation, which requires less computational resources than maximum-likelihood estimation, can still achieve the same optimal performance, characterised by asymptotic efficiency [13,14].

To summarise these propositions, consider a sequence of recursive estimates ${\theta }_n$, computed using a Riemannian stochastic optimisation algorithm with decreasing step sizes (n is the number of observations already processed by the algorithm). Informally, under assumptions which guarantee that ${\theta }^{*}$ is an attractive local minimum of $D\left(\theta \right)$, and that the algorithm is neither too noisy, nor too unstable, in the neighborhood of ${\theta }^{*}$,

• Proposition 2 states that, with an adequate choice of step sizes, the ${\theta }_n$ achieve a fast non-asymptotic rate of convergence to ${\theta }^{*}$. Precisely, the expectation of the squared Riemannian distance between ${\theta }_n$ and ${\theta }^{*}$ is $O\left(n^{-1}\right)$. This is called a fast rate, because it is the best achievable, for any step sizes which are proportional to $n^{-q}$ with $q\in (1/2,1]$ [9,15]. Here, this rate is obtained without any convexity assumptions, for twice differentiable $D\left(\theta \right)$. It would still hold for non-differentiable, but strongly convex, $D\left(\theta \right)$ [2].

• Proposition 4 states that the distribution of the ${\theta }_n$ becomes asymptotically normal, centred at ${\theta }^{*}$, when n grows increasingly large, and also characterises the corresponding asymptotic covariance matrix. This proposition is proved using a novel linearisation technique, which also plays a central role in [3].

• Proposition 5 states that, if the Riemannian manifold $\mathsf{\Theta }$ is equipped with the Fisher information metric of the statistical model P, then Riemannian gradient descent with respect to this information metric, when used to minimise $D\left(\theta \right)$, computes recursive estimates ${\theta }_n$ which are asymptotically efficient, achieving the optimal asymptotic rate of convergence, given by the Cramér-Rao lower bound. This is illustrated, with a numerical application to the recursive estimation of elliptically contoured distributions, in Section 4.

The end result of Proposition 5 is asymptotic efficiency, achieved using the Fisher information metric. In [3], an alternative route to asymptotic efficiency is proposed, using the averaged Riemannian stochastic gradient method. This method does not require any prior knowledge of the Fisher information metric, but has an additional computational cost, which comes from computing on-line Riemannian averages.

The proofs of Propositions 2, 4, and 5, are detailed in Section 6, and Appendix A and Appendix B. Necessary background, about the Fisher information metric (in short, this will be called the information metric), is recalled in Appendix C. Before going on, the reader should note that the summation convention of differential geometry is used throughout the following, when working in local coordinates.

2. Problem Statement

Let $P=(P,\mathsf{\Theta },X)$ be a statistical model, with parameter space $\mathsf{\Theta }$ and sample space X. To each $\theta \in \mathsf{\Theta }$, the model P associates a probability distribution $P_{\theta }$ on X. Here, $\mathsf{\Theta }$ is a $C^r$ Riemannian manifold with $r>3$, and X is any measurable space. The Riemannian metric of $\mathsf{\Theta }$ will be denoted $⟨·,·⟩$, with its Riemannian distance $d(·,·)$. In general, the metric $⟨·,·⟩$ is not the information metric of the model P.

Let $(\Omega ,\mathcal{F},\mathbb{P})$ be a complete probability space, and $(x_n;n=1,2,\dots )$ be i.i.d. random variables on $\Omega$, with values in X. While the distribution of $x_n$ is unknown, it is assumed to belong to the model P. That is, $\mathbb{P}\circ x_n^{-1}=P_{{\theta }^{*}}$ for some ${\theta }^{*}\in \mathsf{\Theta }$, to be called the true parameter.

Consider the following problem: how to obtain fast, asymptotically efficient, recursive estimates ${\theta }_n$ of the true parameter ${\theta }^{*}$, based on observations of the random variables $x_n$? The present work proposes to solve this problem through a detailed study of the decreasing-step-size algorithm, which computes, similar to [1]

(1a)${\theta }_{n+1}={\mathrm{Exp}}_{{\theta }_n}\left({\gamma }_{n+1}u({\theta }_n,x_{n+1})\right)n=0,1,\dots$

starting from an initial guess ${\theta }_0$.

This algorithm has three ingredients. First, $\mathrm{Exp}$ denotes the Riemannian exponential map of the metric $⟨·,·⟩$ of $\mathsf{\Theta }$ [16]. Second, the step sizes ${\gamma }_n$ are strictly positive, decreasing, and verify the usual conditions for stochastic approximation [10,17]

(1b)$\sum {\gamma }_n=\infty \sum {\gamma }_n^2<\infty$

Third, $u(\theta ,x)$ is a continuous vector field on $\mathsf{\Theta }$ for each $x\in X$, which generalises the classical concept of score statistic [13,18]. It will become clear, from the results given in Section 3, that the solution of the above-stated problem depends on the choice of each one of these three ingredients.

A priori knowledge about the model P is injected into Algorithm (1a) using a divergence function $D\left(\theta \right)=D(P_{{\theta }^{*}},P_{\theta })$ (note that ${\theta }^{*}$ is unknown, though). As defined in [19], this is a positive function, equal to zero if and only if $P_{\theta }=P_{{\theta }^{*}}$, and with positive definite Hessian at $\theta ={\theta }^{*}$. Since one expects that minimising $D\left(\theta \right)$ will lead to estimating ${\theta }^{*}$, it is natural to require that

(1c)$E_{{\theta }^{*}}u(\theta ,x)=-\nabla D\left(\theta \right)$

In other words, that $u(\theta ,x)$ is an unbiased estimator of minus the Riemannian gradient of $D\left(\theta \right)$. With $u(\theta ,x)$ given by (1c), Algorithm (1a) is a Riemannian stochastic gradient descent, of the form considered in [1,2,3]. However, as explained in Remark 2, (1c) may be replaced by the weaker condition (9), which states that $D\left(\theta \right)$ is a Lyapunov function of Algorithm (1a), without affecting the results in Section 3. In this sense, Algorithm (1a) is more general than Riemannian stochastic gradient descent.

In practice, a suitable choice of $D\left(\theta \right)$ is often the Kullback-Leibler divergence [20],

(2a)$D\left(\theta \right)=-E_{{\theta }^{*}}logL\left(\theta \right)L\left(\theta \right)=\frac{d{P}_{\theta }}{d{P}_{{\theta }^{*}}}$

where $P_{\theta }$ is absolutely continuous with respect to $P_{{\theta }^{*}}$ with Radon-Nikodym derivative $L\left(\theta \right)$ (the likelihood function). Indeed, if $D\left(\theta \right)$ is chosen to be the Kullback-Leibler divergence, then (1c) is satisfied by

(2b)$u(\theta ,x)=\nabla logL\left(\theta \right)$

which, in many practical situations, can be evaluated directly, without any knowledge of ${\theta }^{*}$.

Before stating the main results, in the following Section 3, it may be helpful to recall some general background on recursive estimation [10]. For simplicity, let $D\left(\theta \right)$ be the Kullback-Leibler divergence (2a). The problem of estimating the true parameter ${\theta }^{*}$ is equivalent to the problem of finding a global minimum of $D\left(\theta \right)$. Of course, this problem cannot be tackled directly, since $D\left(\theta \right)$ cannot be computed without knowledge of ${\theta }^{*}$. There exist two routes around this difficulty.

The first route is empirical minimisation, which replaces the expectation in (2a) with an empirical mean over observed data. Given the first n observations, instead of minimising $D\left(\theta \right)$, one minimises the empirical divergence $D_n\left(\theta \right)$,

(3)$D_n\left(\theta \right)=-\frac{1}{n}\sum _{m=1}^{n}logL(\theta ,x_m)$

where L is the likelihood function of (2a). Now, given the minus sign ahead of the sum in (3), it is clear that minimising $D_n\left(\theta \right)$ amounts to maximising the sum of log-likelihoods. Thus, one is lead to the method of maximum-likelihood estimation.

It is well-known that maximum-likelihood estimation under general regularity conditions is asymptotically efficient [13]. Roughly, this means the maximum-likelihood estimator has the least possible asymptotic variance, equal to the inverse of the Fisher information. On the other hand, as the number n of observations grows, it can be especially difficult to deal with the empirical divergence $D_n\left(\theta \right)$ of Equation (3). In the process of searching for the minimum of $D_n\left(\theta \right)$, each evaluation of this function, or of its derivatives, will involve a massive number of operations, ultimately becoming unpractical.

Aiming to avoid this difficulty, the second route recursive estimation is based on observation-driven updates, following the general scheme of algorithm (1a). In this scheme, each new recursive estimate ${\theta }_{n+1}$ is computed using only the new observation $x_{n+1}$. Therefore, as the number of observations grows very large, the overall computational effort required by recursive estimation remains the same.

The main results in the following section show that recursive estimation can achieve the same asymptotic performance as maximum-likelihood estimation as the number n of observations grows large. As a word of caution, it should be mentioned that recursive estimation does not, in general, fare better than maximum-likelihood estimation for moderate values of the number n of observations, and models with a small number of parameters. However, it is a very desirable substitute to maximum-likelihood estimation for models with a large number of parameters, which typically require a very large number of observations in order to be estimated correctly.

3. Main Results

The motivation of the following Propositions 1 to 5 is to provide general conditions, which guarantee that Algorithm (1a) computes fast, asymptotically efficient, recursive estimates ${\theta }_n$ of the true parameter ${\theta }^{*}$. In the statement of these propositions, it is implicitly assumed that conditions (1b) and (1c) are verified. Moreover, the following assumptions are considered.

• (d1) 

  the divergence function $D\left(\theta \right)$ has an isolated stationary point at $\theta ={\theta }^{*}$, and Lipschitz gradient in a neighborhood of this point.

• (d2) 

  this stationary point is moreover attractive: $D\left(\theta \right)$ is twice differentiable at $\theta ={\theta }^{*}$, with positive definite Hessian at this point.

• (u1) 

  in a neighborhood of $\theta ={\theta }^{*}$, the function $V\left(\theta \right)=E_{{\theta }^{*}}{\parallel u(\theta ,x)\parallel }^2$ is uniformly bounded.

• (u2) 

  in a neighborhood of $\theta ={\theta }^{*}$, the function $R\left(\theta \right)=E_{{\theta }^{*}}{\parallel u(\theta ,x)\parallel }^4$ is uniformly bounded.

For Assumption (d1), the definition of a Lipschitz vector field on a Riemannian manifold may be found in [21]. For Assumptions (u1) and (u2), $\parallel ·\parallel$ denotes the Riemannian norm.

Assumptions (u1) and (u2) are so-called moment control assumptions. They imply that the noise in Algorithm (1a) does not cause the iterates ${\theta }_n$ to diverge, and are also crucial to proving the asymptotic normality of these iterates.

Let ${\mathsf{\Theta }}^{*}$ be a neighborhood of ${\theta }^{*}$ which verifies (d1), (u1), and (u2). Without loss of generality, it is assumed that ${\mathsf{\Theta }}^{*}$ is compact and convex (see the definition of convexity in [16,22]). Then, ${\mathsf{\Theta }}^{*}$ admits a system of normal coordinates $({\theta }^{\alpha };\alpha =1,\dots ,d)$ with origin at ${\theta }^{*}$. With respect to these coordinates, denote the components of $u({\theta }^{*},x)$ by $u^{\alpha }\left({\theta }^{*}\right)$ and let ${\sum }^{*}=\left({\sum }_{\alpha \beta }^{*}\right)$,

(4a)${\sum }_{\alpha \beta }^{*}=E_{{\theta }^{*}}\left[u^{\alpha }\left({\theta }^{*}\right)u^{\beta }\left({\theta }^{*}\right)\right]$

When (d2) is verified, denote the components of the Hessian of $D\left(\theta \right)$ at $\theta ={\theta }^{*}$ by $H=\left(H_{\alpha \beta }\right)$,

(4b)$H_{\alpha \beta }={\left.\frac{{\partial }^{2}D}{\phantom{(}\partial {\theta }^{\alpha }\partial {\theta }^{\beta }}\right|}_{{\theta }^{\alpha }=0}$

Then, the matrix $H=\left(H_{\alpha \beta }\right)$ is positive definite [23]. Denote by $\lambda >0$ its smallest eigenvalue.

Propositions 1 to 5 require the condition that the recursive estimates ${\theta }_n$ are stable, which means that all the ${\theta }_n$ lie in ${\mathsf{\Theta }}^{*}$, almost surely. The need for this condition is discussed in Remark 3. Note that, if ${\theta }_n$ lies in ${\mathsf{\Theta }}^{*}$, then ${\theta }_n$ is determined by its normal coordinates ${\theta }_n^{\alpha }$.

Proposition 1

(consistency). assume (d1) and (u1) are verified, and the recursive estimates ${\theta }_n$ are stable. Then, $lim{\theta }_n={\theta }^{*}$ almost surely.

Proposition 2

(mean-square rate). assume (d1), (d2) and (u1) are verified, the recursive estimates ${\theta }_n$ are stable, and ${\gamma }_n=\frac{a}{n}$ where $2\lambda a>1$. Then

(5)$\mathbb{E}{d}^2({\theta }_n,{\theta }^{*})=O\left(n^{-1}\right)$

Proposition 3

(almost-sure rate). assume the conditions of Proposition 2 are verified. Then,

(6)$d^2({\theta }_n,{\theta }^{*})=o\left(n^{-p}\right)forp\in (0,1)almostsurely$

Proposition 4

(asymptotic normality). assume the conditions of Proposition 2, as well as (u2), are verified. Then, the distribution of the re-scaled coordinates $\left(n^{1/2}{\theta }_n^{\alpha }\right)$ converges to a centred d-variate normal distribution, with covariance matrix Σ given by Lyapunov’s equation

(7)$A\sum +\sum A=-a^2{\sum }^{*}$

where $A=\left(A_{\alpha \beta }\right)$ with $A_{\alpha \beta }=\frac{1}{2}{\delta }_{\alpha \beta }-a{H}_{\alpha \beta }$ (here, δ denotes Kronecker’s delta).

Proposition 5

(asymptotic efficiency). assume the Riemannian metric $⟨·,·⟩$ of Θ coincides with the information metric of the model P, and let $D\left(\theta \right)$ be the Kullback-Leibler divergence (2a). Further, assume (d1), (d2), (u1) and (u2) are verified, the recursive estimates ${\theta }_n$ are stable, and ${\gamma }_n=\frac{a}{n}$ where $2a>1$. Then,

(i) the rates of convergence (5) and (6) hold true.

(ii) if $a=1$, the distribution of the re-scaled coordinates $\left(n^{1/2}{\theta }_n^{\alpha }\right)$ converges to a centred d-variate normal distribution, with covariance matrix ${\sum }^{*}$.

(iii) if $a=1$, and $u(\theta ,x)$ is given by (2b), then ${\sum }^{*}$ is the identity matrix, and the recursive estimates ${\theta }_n$ are asymptotically efficient.

(iv) the following rates of convergence also hold

(8a) $\begin{array}{c}\hfill \mathbb{E}D\left({\theta }_n\right)=O\left(n^{-1}\right)\end{array}$

(8b) $\begin{array}{c}\hfill D\left({\theta }_n\right)=o\left(n^{-p}\right)forp\in (0,1)almostsurely\end{array}$

The following remarks are concerned with the scope of Assumptions (d1), (d2), (u1), and (u2), and with the applicability of Propositions 1 to 5.

Remark 1. 

(d2), (u1) and (u2) do not depend on the Riemannian metric $⟨·,·⟩$ of Θ. Precisely, if they are verified for one Riemannian metric on Θ, then they are verified for any Riemannian metric on Θ. Moreover, if the function $D\left(\theta \right)$ is $C^2$, then the same is true for (d1). In this case, Propositions 1 to 5 apply for any Riemannian metric on Θ, so that the choice of the metric $⟨·,·⟩$ is a purely practical matter, to be decided according to applications.

Remark 2.

the conclusion of Proposition 1 continues to hold, if (1c) is replaced by

(9) $E_{{\theta }^{*}}⟨u(\theta ,x),\nabla D\left(\theta \right)⟩<0for\theta \ne {\theta }^{*}$

Then, it is even possible to preserve Propositions 2, 3, and 4, provided (d2) is replaced by the assumption that the mean vector field, $X\left(\theta \right)=E_{{\theta }^{*}}u(\theta ,x)$, has an attractive stationary point at $\theta ={\theta }^{*}$. This generalisation of Propositions 1 to 4 can be achieved following essentially the same approach as laid out in Section 6. However, in the present work, it will not be carried out in detail.

Remark 3.

the condition that the recursive estimates ${\theta }_n$ are stable is standard in all prior work on stochastic optimisation in manifolds [1,2,3]. In practice, this condition can be enforced through replacing Algorithm (1a) by a so-called projected or truncated algorithm. This is identical to (1a), except that ${\theta }_n$ is projected back onto the neighborhood ${\mathsf{\Theta }}^{*}$ of ${\theta }^{*}$, whenever it falls outside of this neighborhood [10,17]. On the other hand, if the ${\theta }_n$ are not required to be stable, but (d1) and (u1) are replaced by global assumptions,

• (d1’) 

  $D\left(\theta \right)$ has compact level sets and globally Lipschitz gradient.

• (u1’) 

  $V\left(\theta \right)\le C(1+D(\theta \left)\right)$ for some constant C and for all $\theta \in \mathsf{\Theta }$.

then, applying the same arguments as in the proof of Proposition 1, it follows that the ${\theta }_n$ converge to the set of stationary points of $D\left(\theta \right)$, almost surely.

Remark 4.

from (ii) and (iii) of Proposition 5, it follows that the distribution of $n{d}^2({\theta }_n,{\theta }^{*})$ converges to a ${\chi }^2$-distribution with d degrees of freedom. This provides a practical means of confirming the asymptotic efficiency of the recursive estimates ${\theta }_n$.

4. Application: Estimation of ECD

Here, the conclusion of Proposition 5 is illustrated, by applying Algorithm (1a) to the estimation of elliptically contoured distributions (ECD) [24,25]. Precisely, in the notation of Section 2, let $\mathsf{\Theta }={\mathcal{P}}_m$ the space of $m\times m$ positive definite matrices, and $X={\mathbb{R}}^m$. Moreover, let each $P_{\theta }$ have probability density function

(10)$p\left(x\right|\theta )\propto exp\left[h\left(x^{†}{\theta }^{-1}x\right)-\frac{1}{2}logdet\left(\theta \right)\right]\theta \in {\mathcal{P}}_m,x\in {\mathbb{R}}^m$

where $h:\mathbb{R}\to \mathbb{R}$ is fixed, has negative values, and is decreasing, and ${}^{†}$ denotes the transpose. Then, $P_{\theta }$ is called an ECD with scatter matrix $\theta$. To begin, let $(x_n;n=1,2,\dots )$ be i.i.d. random vectors in ${\mathbb{R}}^m$, with distribution $P_{{\theta }^{*}}$ given by (10), and consider the problem of estimating the true scatter matrix ${\theta }^{*}$. The standard approach to this problem is based on maximum-likelihood estimation [25,26]. An original approach, based on recursive estimation, is now introduced using Algorithm (1a).

As in Proposition 5, the parameter space ${\mathcal{P}}_m$ will be equipped with the information metric of the statistical model P just described. In [27], it is proved that this information metric is an affine-invariant metric on ${\mathcal{P}}_m$. In other words, it is of the general form [28]

(11a)${⟨u,u⟩}_{\theta }=I_1\mathrm{tr}{\left({\theta }^{-1}u\right)}^2+I_2{\mathrm{tr}}^2\left({\theta }^{-1}u\right)u\in T_{\theta }{\mathcal{P}}_m$

parameterised by constants $I_1>0$ and $I_2\ge 0$, where $\mathrm{tr}$ denotes the trace and ${\mathrm{tr}}^2$ the squared trace, and where $T_{\theta }{\mathcal{P}}_m$ denotes the tangent space at $\theta$ to the manifold ${\mathcal{P}}_m$. Precisely [27], for the information metric of the model P,

(11b)$I_1=\frac{\phi }{2m^2(m+2)}{I}_2=\frac{\phi }{m^2}-\frac{1}{4}$

where $\phi$ is a further constant, given by the expectation

(11c)$\phi =E_e{\left[h^{\prime }\left(x^{†}x\right)\left(x^{†}x\right)\right]}^2$

with $e\in {\mathcal{P}}_m$ the identity matrix, and $h^{\prime }$ the derivative of h. This expression of the information metric can now be used to specify Algorithm (1a).

First, since the information metric is affine-invariant, it is enough to recall that all affine-invariant metrics on ${\mathcal{P}}_m$ have the same Riemannian exponential map [25,29],

(12a)${\mathrm{Exp}}_{\theta }\left(u\right)=\theta exp\left({\theta }^{-1}u\right)$

where exp denotes the matrix exponential. Second, as in (ii) of Proposition 5, choose the sequence of step sizes

(12b)${\gamma }_n=\frac{1}{n}$

Third, as in (iii) of Proposition 5, let $u(\theta ,x)$ be the vector field on ${\mathcal{P}}_m$ given by (2b),

(12c)$u(\theta ,x)={\nabla }^{\left(inf\right)}logL\left(\theta \right)={\nabla }^{\left(inf\right)}logp\left(x\right|\theta )$

where ${\nabla }^{\left(inf\right)}$ denotes the gradient with respect to the information metric, and $L\left(\theta \right)$ is the likelihood ratio, equal to $p\left(x\right|\theta )$ divided by $p\left(x\right|{\theta }^{*})$. Now, replacing (12) into (1a) defines an original algorithm for recursive estimation of the true scatter matrix ${\theta }^{*}$.

To apply this algorithm in practice, one may evaluate $u(\theta ,x)$ via the following steps. Denote $g(\theta ,x)$ the gradient of $logp\left(x\right|\theta )$ with respect to the affine-invariant metric of [29], which corresponds to $I_1=1$ and $I_2=0$. By direct calculation from (10), this is given by

(13a)$g(\theta ,x)=-\frac{1}{2}\theta -h^{\prime }\left(x^{†}{\theta }^{-1}x\right)x{x}^{†}$

Moreover, introduce the constants $J_1=I_1$ and $J_2=I_1+m{I}_2$. Then, $u(\theta ,x)$ can be evaluated,

(13b)$u(\theta ,x)=J_1^{-1}{\left(g(\theta ,x)\right)}^{\perp }+J_2^{-1}{\left(g(\theta ,x)\right)}^{\Vert }$

from the orthogonal decomposition of $g=g(\theta ,x)$,

(13c)$g^{\Vert }=\mathrm{tr}\left({\theta }^{-1}g\right)\frac{\theta }{m}{g}^{\perp }=g-g^{\Vert }$

Figure 1 and Figure 2 below display numerical results from an application to Kotz-type distributions, which correspond to $h\left(t\right)=-\frac{t^s}{\phantom{(}2}$ in (10) and $\phi =s^2\frac{m}{2s}\left(\frac{m}{2s}+1\right)$ in (11c) [24,27]. These figures were generated from ${10}^3$ Monte Carlo runs of the algorithm defined by (1a) and (12), with random initialisation, for the specific values $s=4$ and $m=7$. Essentially the same numerical results could be observed for any $s\le 9$ and $m\le 50$.

Figure 1 confirms the fast non-asymptotic rate of convergence (5), stated in (i) of Proposition 5. On a log-log scale, it shows the empirical mean ${\mathbb{E}}_{\mathrm{MC}}{d}^2({\theta }_n,{\theta }^{*})$ over Monte Carlo runs, as a function of n. This decreases with a constant negative slope equal to $-1$, starting roughly at $logn=4$. Here, the Riemannian distance $d({\theta }_n,{\theta }^{*})$ induced by the information metric (11) is given by [28]

(14)$d^2(\theta ,{\theta }^{*})=I_1\mathrm{tr}{\left[log\left({\theta }^{-1}{\theta }^{*}\right)\right]}^2+I_2{\mathrm{tr}}^2\left[log\left({\theta }^{-1}{\theta }^{*}\right)\right]\theta ,{\theta }^{*}\in \mathsf{\Theta }$

where log denotes the symmetric matrix logarithm [30]. Figure 2 confirms the asymptotic efficiency of the recursive estimates ${\theta }_n$, stated in (iii) of Proposition 5, using Remark 4. It shows a kernel density estimate of $n{d}^2({\theta }_n,{\theta }^{*})$ where $n={10}^5$ (solid blue curve). This agrees with a ${\chi }^2$-distribution with 28 degrees of freedom (dotted red curve), where $d=28$ is indeed the dimension of the parameter space ${\mathcal{P}}_m$ for $m=7$.

5. Conclusions

Recursive estimation is a subject that is over fifty years old [10], with its foundation in the general theory of stochastic optimisation [9,15]. Its applications are very wide-ranging, as they cover areas from control theory to machine learning [17].

With the increasing role of Riemannian manifolds in statistical inference and machine learning, it was natural to generalise the techniques of stochastic optimisation, from Euclidean space to Riemannian manifolds. Indeed, this started with the work of Bonnabel [1], which impulsed a series of successive works, such as [2,3].

These works have mostly sought to directly adapt classical results, known in Euclidean space, which concern optimal rates of convergence to a unique attractive minimum of a cost function. The present work also belongs to this line of thinking. It shows that when dealing with a recursive estimation problem, where the unknown statistical parameter belongs to a differentiable manifold, equipping this manifold with the information metric of the underlying statistical model, leads to optimal algorithm performance, which is moreover automatic (does not involve parameter tuning).

The results obtained in the present work (as well as in [2,3]) suffer from inherent limitations. Indeed, being only focused on convergence to a unique attractive minimum, it does not tackle the following important open problems:

• stability of stochastic optimisation algorithms finding verifiable conditions which ensure a stochastic optimisation algorithm remains within a compact set. A more general form of this problem is computing the probability of a stochastic optimisation algorithm exiting a certain neighborhood of a stationary point (whether attractive or not) within a finite number of iterations.

• non-asymptotic performance of stochastic optimisation algorithms: this involves computing explicitly the outcome which the algorithm is able to achieve, after a given finite number of iterations. This provides a much stronger theoretical guarantee, to the user, than standard results which compute a rate of convergence.

These problems have attracted much attention and generated well-known results when considered in the Euclidean case [31,32], but remain open in the context of Riemannian manifolds. They involve much richer interaction between Riemannian geometry and stochastic optimisation, due to their global nature.

6. Proofs of Main Results

Proof of Proposition 1. 

The proof is a generalisation of the original proof in [1], itself modeled on the proof for the Euclidean case in [33]. Throughout the following, let ${\mathcal{X}}_n$ be the $\sigma$-field generated by $x_1,\dots ,x_n$ [20]. Recall that $(x_n;n=1,2,\dots )$ are i.i.d. with distribution $P_{{\theta }^{*}}$. Therefore, by (1a), ${\theta }_n$ is ${\mathcal{X}}_n$-measurable and $x_{n+1}$ is independent from ${\mathcal{X}}_n$. Thus, using elementary properties of conditional expectation [20],

(15a)$\begin{array}{c}\hfill \mathbb{E}\left[u({\theta }_n,x_{n+1})|{\mathcal{X}}_n\right]=-D\left({\theta }_n\right)\end{array}$

(15b)$\begin{array}{c}\hfill \mathbb{E}\left[\parallel u({\theta }_n,x_{n+1}){\parallel }^2|{\mathcal{X}}_n\right]=V\left({\theta }_n\right)\end{array}$

where (15a) follows from (1c), and (15b) from (u1). Let L be a Lipschitz constant for $\nabla D\left(\theta \right)$, and C be an upper bound on $V\left(\theta \right)$, for $\theta \in {\mathsf{\Theta }}^{*}$. The following inequality is now proved, for any positive integer n,

(16)$\mathbb{E}\left[D\left({\theta }_{n+1}\right)-D\left({\theta }_n\right)|{\mathcal{X}}_n\right]\le {\gamma }_{n+1}^{2}LC-{\gamma }_{n+1}{\parallel \nabla D\left({\theta }_n\right)\parallel }^2$

once this is done, Proposition 1 is obtained by applying the Robbins-Siegmund theorem [9].

Proof of (16): let $c\left(t\right)$ be the geodesic connecting ${\theta }_n$ to ${\theta }_{n+1}$ with equation

(17a)$c\left(t\right)={\mathrm{Exp}}_{{\theta }_n}\left(t{\gamma }_{n+1}u({\theta }_n,x_{n+1})\right)$

From the fundamental theorem of calculus,

(17b)$D\left({\theta }_{n+1}\right)-D\left({\theta }_n\right)={\gamma }_{n+1}⟨u({\theta }_n,x_{n+1}),\nabla D\left({\theta }_n\right)⟩+{\gamma }_{n+1}{\int }_0^1\left[{⟨\dot{c},\nabla D⟩}_{c\left(t\right)}-{⟨\dot{c},\nabla D⟩}_{c\left(0\right)}\right]dt$

Since the recursive estimates ${\theta }_n$ are stable, ${\theta }_n$ and ${\theta }_{n+1}$ both lie in ${\mathsf{\Theta }}^{*}$. Since ${\mathsf{\Theta }}^{*}$ is convex, the whole geodesic $c\left(t\right)$ lies in ${\mathsf{\Theta }}^{*}$. Then, since $\nabla D\left(\theta \right)$ is Lipschitz on ${\mathsf{\Theta }}^{*}$, it follows from (17b),

(17c)$D\left({\theta }_{n+1}\right)-D\left({\theta }_n\right)\le {\gamma }_{n+1}⟨u({\theta }_n,x_{n+1}),\nabla D\left({\theta }_n\right)⟩+{\gamma }_{n+1}^{2}L{\parallel u({\theta }_n,x_{n+1})\parallel }^2$

Taking conditional expectations in this inequality, and using (15a) and (15b),

(17d)$\mathbb{E}\left[D\left({\theta }_{n+1}\right)-D\left({\theta }_n\right)|{\mathcal{X}}_n\right]\le -{\gamma }_{n+1}{\parallel \nabla D\left({\theta }_n\right)\parallel }^2+{\gamma }_{n+1}^{2}LV\left({\theta }_n\right)$

so (16) follows since (u1) guarantees $V\left({\theta }_n\right)\le C$. □ Conclusion: by the Robbins-Siegmund theorem, inequality (16) implies that, almost surely,

(18a)$limD\left({\theta }_n\right)=D_{\infty }<\infty and\sum _{n=1}^{\infty }{\gamma }_{n+1}{\parallel \nabla D\left({\theta }_n\right)\parallel }^2<\infty$

In particular, from the first condition in (1b), convergence of the sum in (18a) implies

(18b)$lim\parallel \nabla D\left({\theta }_n\right)\parallel =0\mathrm{almost}\mathrm{surely}$

Now, since the sequence of recursive estimates ${\theta }_n$ lies in the compact set ${\mathsf{\Theta }}^{*}$, it has at least one point of accumulation in this set, say ${\theta }_{*}$. If ${\theta }_{n\left(k\right)}$ is a subsequence of ${\theta }_n$, converging to ${\theta }_{*}$,

$\parallel \nabla D\left({\theta }_{*}\right)\parallel =lim\parallel \nabla D\left({\theta }_{n\left(k\right)}\right)\parallel =lim\parallel \nabla D\left({\theta }_n\right)\parallel =0\mathrm{almost}\mathrm{surely}$

where the third equality follows from (18b). This means that ${\theta }_{*}$ is a stationary point of $D\left(\theta \right)$ in ${\mathsf{\Theta }}^{*}$. Thus, (d1) implies ${\theta }_{*}={\theta }^{*}$ is the unique point of accumulation of ${\theta }_n$. In other words, $lim{\theta }_n={\theta }^{*}$ almost surely. □

Proof of Proposition 2. 

The proof is modeled on the proofs for the Euclidean case, given in [10,15]. It relies on the following geometric Lemmas 1 and 2. Lemma 1 will be proved in Appendix A. On the other hand, Lemma 2 is the same as the trigonometric distance bound of [2]. For Lemma 1, recall that $\lambda >0$ denotes the smallest eigenvalue of the matrix H defined in (4b).

Lemma 1. 

for any $\mu <\lambda$, there exists a neighborhood ${\overline{\mathsf{\Theta }}}^{*}$ of ${\theta }^{*}$, contained in ${\mathsf{\Theta }}^{*}$, with

(19a) $⟨{\mathrm{Exp}}_{\theta }^{-1}\left({\theta }^{*}\right),\nabla D\left(\theta \right)⟩\le -\mu d^2(\theta ,{\theta }^{*})for\theta \in {\overline{\mathsf{\Theta }}}^{*}$

Lemma 2. 

let $-{\kappa }^2$ be a lower bound on the sectional curvature of Θ in ${\mathsf{\Theta }}^{*}$, and $C_{\kappa }=R\kappa coth\left(R\kappa \right)$ where R is the diameter of ${\mathsf{\Theta }}^{*}$. For $\tau ,\theta \in {\mathsf{\Theta }}^{*}$, where $\tau ={\mathrm{Exp}}_{\theta }\left(u\right)$,

(19b) $d^2(\tau ,{\theta }^{*})\le d^2(\theta ,{\theta }^{*})-2⟨{\mathrm{Exp}}_{\theta }^{-1}\left({\theta }^{*}\right),u⟩+C_{\kappa }{\parallel u\parallel }^2$

Proof of (5): let ${\gamma }_n=\frac{a}{n}$ with $2\lambda a>2\mu a>1$ for some $\mu <\lambda$, and let ${\overline{\mathsf{\Theta }}}^{*}$ be the neighborhood corresponding to $\mu$ in Lemma 1. By Proposition 1, the ${\theta }_n$ converge to ${\theta }^{*}$ almost surely. Without loss of generality, it can be assumed that all the ${\theta }_n$ lie in ${\overline{\mathsf{\Theta }}}^{*}$, almost surely. Then, (1a) and Lemma 2 imply, for any positive integer n,

(20a)$d^2({\theta }_{n+1},{\theta }^{*})\le d^2({\theta }_n,{\theta }^{*})-2{\gamma }_{n+1}⟨{\mathrm{Exp}}_{{\theta }_n}^{-1}\left({\theta }^{*}\right),u({\theta }_n,x_{n+1})⟩+{\gamma }_{n+1}^2{C}_{\kappa }{\parallel u({\theta }_n,x_{n+1})\parallel }^2$

Indeed, this follows by replacing $\tau ={\theta }_{n+1}$ and $\theta ={\theta }_n$ in (19b). Taking conditional expectations in (20a), and using (15a) and (15b),

$\mathbb{E}\left[d^2({\theta }_{n+1},{\theta }^{*})|{\mathcal{X}}_n\right]\le d^2({\theta }_n,{\theta }^{*})+2{\gamma }_{n+1}⟨{\mathrm{Exp}}_{{\theta }_n}^{-1}\left({\theta }^{*}\right),\nabla D\left({\theta }_n\right)⟩+{\gamma }_{n+1}^2{C}_{\kappa }V\left({\theta }_n\right)$

Then, by (u1) and (19a) of Lemma 1,

(20b)$\mathbb{E}\left[d^2({\theta }_{n+1},{\theta }^{*})|{\mathcal{X}}_n\right]\le d^2({\theta }_n,{\theta }^{*})(1-2{\gamma }_{n+1}\mu )+{\gamma }_{n+1}^2{C}_{\kappa }C$

where C is an upper bound on $V\left(\theta \right)$, for $\theta \in {\mathsf{\Theta }}^{*}$. By further taking expectations

(20c)$\mathbb{E}{d}^2({\theta }_{n+1},{\theta }^{*})\le \mathbb{E}{d}^2({\theta }_n,{\theta }^{*})(1-2{\gamma }_{n+1}\mu )+{\gamma }_{n+1}^2{C}_{\kappa }C$

Using (20c), the proof reduces to an elementary reasoning by recurrence. Indeed, replacing ${\gamma }_n=\frac{a}{n}$ into (20c), it follows that

(21a)$\mathbb{E}{d}^2({\theta }_{n+1},{\theta }^{*})\le \mathbb{E}{d}^2({\theta }_n,{\theta }^{*})\left(1-\frac{2\mu a}{n+1}\right)+\frac{a^2{C}_{\kappa }C}{{(n+1)}^2}$

On the other hand, if $b\left(n\right)=\frac{b}{n}$ where $b>a^2{C}_{\kappa }C{(2\mu a-1)}^{-1}$, then

(21b)$b(n+1)\ge b\left(n\right)\left(1-\frac{2\mu a}{n+1}\right)+\frac{a^2{C}_{\kappa }C}{{(n+1)}^2}$

Let b be sufficiently large, so (21b) is verified and $\mathbb{E}{d}^2({\theta }_{n_o},{\theta }^{*})\le b\left(n_o\right)$ for some $n_o$. Then, by recurrence, using (21a) and (21b), one also has that $\mathbb{E}{d}^2({\theta }_n,{\theta }^{*})\le b\left(n\right)$ for all $n\ge n_o$. In other words, (5) holds true. □

Proof of Proposition 3. 

the proof is modeled on the proof for the Euclidean case in [10]. To begin, let $W_n$ be the stochastic process given by

(22a)$W_n=n^p{d}^2({\theta }_n,{\theta }^{*})+n^{-q}\mathrm{where}q\in (0,1-p)$

The idea is to show that this process is a positive supermartingale, for sufficiently large n. By the supermartingale convergence theorem [20], it then follows that $W_n$ converges to a finite limit, almost surely. In particular, this implies

(22b)$lim{n}^p{d}^2({\theta }_n,{\theta }^{*})={\ell }_p<\infty \mathrm{almost}\mathrm{surely}$

Then, ${\ell }_p$ must be equal to zero, since p is arbitrary in the interval $(0,1)$. Precisely, for any $\epsilon \in (0,1-p)$,

${\ell }_p=lim{n}^p{d}^2({\theta }_n,{\theta }^{*})=lim{n}^{-\epsilon }{n}^{p+\epsilon }{d}^2({\theta }_n,{\theta }^{*})=\left(lim{n}^{-\epsilon }\right){\ell }_{p+\epsilon }=0$

It remains to show that $W_n$ is a supermartingale, for sufficiently large n. To do so, note that by (20b) from the proof of Proposition 2,

$\mathbb{E}\left[W_{n+1}-W_n|{\mathcal{X}}_n\right]\le d^2({\theta }_n,{\theta }^{*})\frac{p-2\mu a}{{(n+1)}^{1-p}}+\frac{a^2{C}_{\kappa }C}{{(n+1)}^{2-p}}-\frac{q}{{(n+1)}^{q+1}}$

Here, the first term on the right-hand side is negative, since $2\mu a>1>p$. Moreover, the third term dominates the second one for sufficiently large n, since $q<1-p$. Thus, for sufficiently large n, the right-hand side is negative, and $W_n$ is a supermartingale. □

Proof of Proposition 4. 

the proof relies on the following geometric Lemmas 3 and 4, which are used to linearise Algorithm (1a), in terms of the normal coordinates ${\theta }^{\alpha }$. This idea of linearisation in terms of local coordinates also plays a central role in [3].

Lemma 3. 

let ${\theta }_n,{\theta }_{n+1}$ be given by (1a) with ${\gamma }_n=\frac{a}{n}$. Then, in a system of normal coordinates with origin at ${\theta }^{*}$,

(23a) ${\theta }_{n+1}^{\alpha }={\theta }_n^{\alpha }+{\gamma }_{n+1}^{\phantom{2}}{u}_{n+1}^{\alpha }+{\gamma }_{n+1}^2{\pi }_{n+1}^{\alpha }\mathbb{E}\left|{\pi }_{n+1}^{\alpha }\right|=O\left(n^{-1/2}\right)$

where $u_{n+1}^{\alpha }$ are the components of $u({\theta }_n,x_{n+1})$.

Lemma 4. 

let $v_n=\nabla D\left({\theta }_n\right)$. Then, in a system of normal coordinates with origin at ${\theta }^{*}$,

(23b) $v_n^{\alpha }=H_{\alpha \beta }^{\phantom{2}}{\theta }_n^{\beta }+{\rho }_n^{\alpha }{\rho }_n^{\alpha }=o\left(d({\theta }_n,{\theta }^{*})\right)$

where $v_n^{\alpha }$ are the components of $v_n$ and the $H_{\alpha \beta }$ were defined in (4b).

Linearisation of (1a): let $u({\theta }_n,x_{n+1})=-v_n+w_{n+1}$. Then, it follows from (23a) and (23b),

(24a)${\theta }_{n+1}^{\alpha }={\theta }_n^{\alpha }-{\gamma }_{n+1}^{\phantom{2}}{H}_{\alpha \beta }^{\phantom{2}}{\theta }_n^{\beta }-{\gamma }_{n+1}^{\phantom{2}}{\rho }_n^{\alpha }+{\gamma }_{n+1}^{\phantom{2}}{w}_{n+1}^{\alpha }+{\gamma }_{n+1}^2{\pi }_{n+1}^{\alpha }$

Denote the re-scaled coordinates $n^{1/2}{\theta }_n^{\alpha }$ by ${\eta }_n^{\alpha }$, and recall ${\gamma }_n=\frac{a}{n}$. Then, using the estimate ${(n+1)}^{1/2}=n^{1/2}(1+{\left(2n\right)}^{-1}+O\left(n^{-2}\right))$, it follows from (24a) that

(24b)${\eta }_{n+1}^{\alpha }={\eta }_n^{\alpha }+\frac{A_{\alpha \beta }}{n+1}{\eta }_n^{\beta }+\frac{a}{{(n+1)}^{1/2}}\left[B_{\alpha \beta }^{\phantom{2}}{\theta }_n^{\beta }-{\rho }_n^{\alpha }+w_{n+1}^{\alpha }+\frac{a{\pi }_{n+1}^{\alpha }}{n+1}\right]$

where $A_{\alpha \beta }=\frac{1}{2}{\delta }_{\alpha \beta }-a{H}_{\alpha \beta }$ and $B_{\alpha \beta }=O\left(n^{-1}\right)$. Equation (24b) is a first-order, inhomogeneous, linear difference equation, for the “vector” ${\eta }_n^{\phantom{2}}$ of components ${\eta }_n^{\alpha }$. □

Study of equation (24b): switching to vector-matrix notation, equation (24b) is of the general form

(25a)${\eta }_{n+1}^{\phantom{2}}=\left(I+\frac{A}{n+1}\right){\eta }_n^{\phantom{2}}+\frac{a{\xi }_{n+1}^{\phantom{2}}}{{(n+1)}^{1/2}}$

where I denotes the identity matrix, A has matrix elements $A_{\alpha \beta }$, and $\left({\xi }_n\right)$ is a sequence of inputs. The general solution of this equation is [10,34]

(25b)${\eta }_n^{\phantom{2}}=A_{n,m}^{\phantom{2}}{\eta }_m^{\phantom{2}}+\sum _{k=m+1}^n{A}_{n,k}^{\phantom{2}}\frac{a{\xi }_k^{\phantom{2}}}{k^{1/2}}\mathrm{for}n\ge m$

where the transition matrix $A_{n,k}$ is given by

(25c)$A_{n,k}^{\phantom{2}}=\prod _{j=k+1}^n\left(I+\frac{A}{j}\right)A_{n,n}^{\phantom{2}}=I$

Since $2\lambda a>1$, the matrix A is stable. This can be used to show that [10,34]

(25d)$q>\frac{1}{2}\mathrm{and}\mathbb{E}\left|{\xi }_n^{\phantom{2}}\right|=O\left(n^{-q}\right)⟹lim{\eta }_n^{\phantom{2}}=0\mathrm{in}\mathrm{probability}$

where $|{\xi }_n|$ denotes the Euclidean vector norm. Then, it follows from (25d) that ${\eta }_n^{\phantom{2}}$ converges to zero in probability, in each one of the three cases

${\xi }_{n+1}^{\alpha }=B_{\alpha \beta }^{\phantom{2}}{\theta }_n^{\beta };{\xi }_{n+1}^{\alpha }={\rho }_n^{\alpha };{\xi }_{n+1}^{\alpha }=\frac{{\pi }_{n+1}^{\alpha }}{n+1}$

Indeed, in the first two cases, the condition required in (25d) can be verified using (5), whereas in the third case, it follows immediately from the estimate of $\mathbb{E}|{\pi }_{n+1}^{\alpha }|$ in (23a). □

Conclusion: by linearity of (24b), it is enough to consider the case ${\xi }_{n+1}^{\alpha }=w_{n+1}^{\alpha }$ in (25a). Then, according to (25b), ${\eta }_n^{\phantom{2}}$ has the same limit distribution as the sums

(26)${\tilde{\eta }}_n^{\phantom{2}}=\sum _{k=1}^n{A}_{n,k}^{\phantom{2}}\frac{a{w}_k}{k^{1/2}}$

By (15), $\left(w_k\right)$ is a sequence of square-integrable martingale differences. Therefore, to conclude that the limit distribution of ${\tilde{\eta }}_n^{\phantom{2}}$ is a centred d-variate normal distribution, with covariance matrix $\sum$ given by (7), it is enough to verify the conditions of the martingale central limit theorem [35],

(27a)$lim\underset{k\le n}{max}\left|A_{n,k}^{\phantom{2}}\frac{a{w}_k}{k^{1/2}}\right|=0\mathrm{in}\mathrm{probability}$

(27b)$sup\mathbb{E}{\left|{\tilde{\eta }}_n^{\phantom{2}}\right|}^2<\infty$

(27c)$lim\sum _{k=1}^n\frac{a^2}{k}{A}_{n,k}^{\phantom{2}}{\sum }_k^{\phantom{2}}{A}_{n,k}^{\phantom{2}}=\sum \mathrm{in}\mathrm{probability}$

where ${\sum }_k^{\phantom{2}}$ is the conditional covariance matrix

(28)${\sum }_k^{\phantom{2}}=\mathbb{E}\left[w_k^{\phantom{†}}{w}_k^{†}|{\mathcal{X}}_{k-1}\right]$

Conditions (27) are verified in Appendix B, which completes the proof. □

Proof of Proposition 5. 

Denote ${\partial }_{\alpha }=\frac{\partial }{\phantom{(}\partial {\theta }^{\alpha }}$ the coordinate vector fields of the normal coordinates ${\theta }^{\alpha }$. Since $⟨·,·⟩$ coincides with the information metric of the model P, it follows from (4b) and (A10),

(29a)$H_{\alpha \beta }={⟨{\partial }_{\alpha },{\partial }_{\beta }⟩}_{{\theta }^{*}}$

However, by the definition of normal coordinates [16], the ${\partial }_{\alpha }$ are orthonormal at ${\theta }^{*}$. Therefore,

(29b)$H_{\alpha \beta }={\delta }_{\alpha \beta }$

Thus, the matrix H is equal to the identity matrix, and its smallest eigenvalue is $\lambda =1$.

Proof of (i): this follows directly from Propositions 2 and 3. Indeed, since $\lambda =1$, the conditions of these propositions are verified, as soon as $2a>1$. Therefore, (5) and (6) hold true. □

Proof of (ii): this follows from Proposition 4. The conditions of this proposition are verified, as soon as $2a>1$. Therefore, the distribution of the re-scaled coordinates $\left(n^{1/2}{\theta }_n^{\alpha }\right)$ converges to a centred d-variate normal distribution, with covariance matrix $\sum$ given by Lyapunov’s equation (7). If $a=1$, then (29b) implies $A_{\alpha \beta }=-\frac{1}{2}{\delta }_{\alpha \beta }$, so that Lyapunov’s equation (7) reads $\sum ={\sum }^{*}$, as required. □

For the following proof of (iii), the reader may wish to recall that summation convention is used throughout the present work. That is [16], summation is implicitly understood over any repeated subscript or superscript from the Greek alphabet, taking the values $1,\dots ,d$.

Proof of (iii): let $\ell \left(\theta \right)=logL\left(\theta \right)$ and assume $u(\theta ,x)$ is given by (2b). Then, by the definition of normal coordinates [16], the following expression holds

(30a)$u^{\alpha }\left({\theta }^{*}\right)={\left.\frac{\partial \ell }{\partial {\theta }^{\alpha }}\right|}_{{\theta }^{\alpha }=0}$

Replacing this into (4a) gives

(30b)${\sum }_{\alpha \beta }^{*}=E_{{\theta }^{*}}{\left[\frac{\partial \ell }{\partial {\theta }^{\alpha }}\frac{\partial \ell }{\partial {\theta }^{\beta }}\right]}_{{\theta }^{\alpha }=0}=-E_{{\theta }^{*}}{\left.\frac{{\partial }^2\ell }{\phantom{(}\partial {\theta }^{\alpha }\partial {\theta }^{\beta }}\right|}_{{\theta }^{\alpha }=0}={\left.\frac{{\partial }^{2}D}{\phantom{(}\partial {\theta }^{\alpha }\partial {\theta }^{\beta }}\right|}_{{\theta }^{\alpha }=0}$

where the second equality is the so-called Fisher’s identity (see [19], Page 28), and the third equality follows from (2a) by differentiating under the expectation. Now, by (4b) and (29b), ${\sum }^{*}$ is the identity matrix.

To show that the recursive estimates ${\theta }_n$ are asymptotically efficient, let $({\tau }^{\alpha };\alpha =1,\dots ,d)$ be any local coordinates with origin at ${\theta }^{*}$ and let ${\tau }_n^{\alpha }={\tau }^{\alpha }\left({\theta }_n\right)$. From the second-order Taylor expansion of each coordinate function ${\tau }^{\alpha }$, it is straightforward to show that

(31a)$n^{1/2}{\tau }_n^{\alpha }={\left(\frac{\partial {\tau }^{\alpha }}{\partial {\theta }^{\gamma }}\right)}_{{\theta }^{*}}\left(n^{1/2}{\theta }_n^{\gamma }\right)+{\sigma }^{\alpha }\left({\theta }_n\right)\left(n^{1/2}{d}^2({\theta }_n,{\theta }^{*})\right)$

where the subscript ${\theta }^{*}$ indicates the derivative is evaluated at ${\theta }^{*}$, and where ${\sigma }^{\alpha }$ is a continuous function in the neighborhood of ${\theta }^{*}$. By (6), the second term in (31a) converges to zero almost surely. Therefore, the limit distribution of the re-scaled coordinates $\left(n^{1/2}{\tau }_n^{\alpha }\right)$ is the same as that of the first term in (31a). By (ii), this is a centred d-variate normal distribution with covariance matrix ${\sum }^{\tau }$ given by

(31b)${\sum }_{\alpha \beta }^{\tau }={\left(\frac{\partial {\tau }^{\alpha }}{\partial {\theta }^{\gamma }}\right)}_{{\theta }^{*}}{\sum }_{\gamma \kappa }^{*}{\left(\frac{\partial {\tau }^{\beta }}{\partial {\theta }^{\kappa }}\right)}_{{\theta }^{*}}={\left(\frac{\partial {\tau }^{\alpha }}{\partial {\theta }^{\gamma }}\right)}_{{\theta }^{*}}{\left(\frac{\partial {\tau }^{\beta }}{\partial {\theta }^{\gamma }}\right)}_{{\theta }^{*}}$

where the second equality follows because ${\sum }_{\gamma \kappa }^{*}={\delta }_{\gamma \kappa }$ since ${\sum }^{*}$ is the identity matrix.

It remains to show that ${\sum }^{\tau }$ is the inverse of the information matrix $I^{\tau }$ as in (A12). According to (A10), this is given by

(31c)$I_{\alpha \beta }^{\tau }={\left.\frac{{\partial }^{2}D}{\phantom{(}\partial {\tau }^{\alpha }\partial {\tau }^{\beta }}\right|}_{{\tau }^{\alpha }=0}=-E_{{\theta }^{*}}{\left.\frac{{\partial }^2\ell }{\phantom{(}\partial {\tau }^{\alpha }\partial {\tau }^{\beta }}\right|}_{{\tau }^{\alpha }=0}=E_{{\theta }^{*}}{\left[\frac{\partial \ell }{\partial {\tau }^{\alpha }}\frac{\partial \ell }{\partial {\tau }^{\beta }}\right]}_{{\tau }^{\alpha }=0}$

where the second equality follows from (2a), and the third equality from Fisher’s identity (see [19], Page 28). Now, a direct application of the chain rule yields the following

$I_{\alpha \beta }^{\tau }=E_{{\theta }^{*}}{\left[\frac{\partial \ell }{\partial {\tau }^{\alpha }}\frac{\partial \ell }{\partial {\tau }^{\beta }}\right]}_{{\tau }^{\alpha }=0}={\left(\frac{\partial {\theta }^{\gamma }}{\partial {\tau }^{\alpha }}\right)}_{{\theta }^{*}}{E}_{{\theta }^{*}}{\left[\frac{\partial \ell }{\partial {\theta }^{\gamma }}\frac{\partial \ell }{\partial {\theta }^{\kappa }}\right]}_{{\theta }^{\gamma }=0}{\left(\frac{\partial {\theta }^{\kappa }}{\partial {\tau }^{\beta }}\right)}_{{\theta }^{*}}$

By the first equality in (30b), this is equal to

(31d)$I_{\alpha \beta }^{\tau }={\left(\frac{\partial {\theta }^{\gamma }}{\partial {\tau }^{\alpha }}\right)}_{{\theta }^{*}}{\sum }_{\gamma \kappa }^{*}{\left(\frac{\partial {\theta }^{\kappa }}{\partial {\tau }^{\beta }}\right)}_{{\theta }^{*}}={\left(\frac{\partial {\theta }^{\gamma }}{\partial {\tau }^{\alpha }}\right)}_{{\theta }^{*}}{\left(\frac{\partial {\theta }^{\gamma }}{\partial {\tau }^{\beta }}\right)}_{{\theta }^{*}}$

because ${\sum }_{\gamma \kappa }^{*}={\delta }_{\gamma \kappa }$ is the identity matrix. Comparing (31b) to (31d>), it is clear that ${\sum }^{\tau }$ is the inverse of the information matrix $I^{\tau }$ as in (A12).

Proof of(iv): (8a) and (8b) follow from (5) and (6), respectively, by using (A11). Precisely, it is possible to write (A11) in the form

(32a)$D\left({\theta }_n\right)=\frac{1}{2}{d}^2({\theta }_n,{\theta }^{*})+\omega \left({\theta }_n\right)d^2({\theta }_n,{\theta }^{*})$

where $\omega$ is a continuous function in the neighborhood of ${\theta }^{*}$, equal to zero at $\theta ={\theta }^{*}$. To obtain (8a), it is enough to take expectations in (32a) and note that $\omega$ is bounded above in the neighborhood of ${\theta }^{*}$. Then, (8a) follows directly from (5).

To obtain (8b), it is enough to multiply (32a) by $n^p$ where $p\in (0,1)$. This gives the following expression

(32b)$n^{p}D\left({\theta }_n\right)=\frac{1}{2}{n}^p{d}^2({\theta }_n,{\theta }^{*})\left(1+\omega \left({\theta }_n\right)\right)$

From (6), $n^p{d}^2({\theta }_n,{\theta }^{*})$ converges to zero almost surely. Moreover, by continuity of $\omega$, it follows that $\omega \left({\theta }_n\right)$ converges to $\omega \left({\theta }^{*}\right)=0$ almost surely. Therefore, by taking limits in (32b), it is readily seen that

(32c)$lim{n}^{p}D\left({\theta }_n\right)=\frac{1}{2}\left(lim{n}^p{d}^2({\theta }_n,{\theta }^{*})\right)\left(1+lim\omega \left({\theta }_n\right)\right)=0$

almost surely. However, this is equivalent to the statement that $D\left({\theta }_n\right)=o\left(n^{-p}\right)$ for $p\in (0,1)$, almost surely. Thus, (8b) is proved. □

Author Contributions

Data curation, J.Z.; Software, J.Z.; Writing—original draft, S.S.; Writing—review & editing, S.S.

Funding

Agence Nationale de la Recherche: ANR-17-ASTR-0015.

Acknowledgments

At the end, we thank the support from the ANR MARGARITA (Modern Adaptive Radar: Great Advances in Robust and Inference Techniques and Application) under Grant ANR-17-ASTR-0015 for our works.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Proofs of Geometric Lemmas

Appendix B. Conditions of the Martingale CLT

For the verification of Conditions (27), the following inequality (A5) will be useful. Let $\nu =a\lambda -\frac{1}{2}$, so $-\nu$ is the largest eigenvalue of the matrix A in (25a). There exists a constant $C_A$ such that the transition matrices $A_{n,k}$ in (25c) satisfy [10,34]

(A5)

${\left|A_{n,k}\right|}_{\mathrm{Op}}\le C_A{\left(\frac{k}{n}\right)}^{\nu }$

where ${\left|A_{n,k}\right|}_{\mathrm{Op}}$ denotes the Euclidean operator norm, equal to the largest singular value of the matrix $A_{n,k}$.

Condition (27a): to verify this condition, note that for arbitrary $\epsilon >0$,

(A6a)

$\mathbb{P}\left(\underset{k\le n}{max}\left|A_{n,k}^{\phantom{2}}\frac{a{w}_k}{k^{1/2}}\right|>\epsilon \right)\le \sum _{k=1}^n\mathbb{P}\left(\left|A_{n,k}^{\phantom{2}}\frac{a{w}_k}{k^{1/2}}\right|>\epsilon \right)\le \sum _{k=1}^n\mathbb{P}\left(C_A{\left(\frac{k}{n}\right)}^{\nu }\left|\frac{a{w}_k}{k^{1/2}}\right|>\epsilon \right)$

where the second inequality follows from (A5). However, it follows from (u2) that there exists a uniform upper bound $M_w$ on the fourth-order moments of $|w_k|$. Therefore, by Chebyshev’s inequality [20]

(A6b)

$\sum _{k=1}^n\mathbb{P}\left(C_A{\left(\frac{k}{n}\right)}^{\nu }\left|\frac{a{w}_k}{k^{1/2}}\right|>\epsilon \right)\le {\left(\frac{a{C}_A}{\epsilon }\right)}^4\frac{M_w}{n^{4\nu }}\sum _{k=1}^n{k}^{4\nu -2}$

Since $\nu >0$, the right-hand side of (A6b) has limit equal to 0 as $n\to \infty$, by the Euler–Maclaurin formula [37]. Replacing this limit from (A6b) into (A6a) immediately yields Condition (27a). □

Condition (27b): to verify this condition, recall that $\left(w_k\right)$ is a sequence of square-integrable martingale differences. Therefore, from (26)

(A7a)

$\mathbb{E}{\left|{\tilde{\eta }}_n^{\phantom{2}}\right|}^2=\sum _{k=1}^n\frac{a^2}{k}\mathbb{E}\mathrm{tr}\left(A_{n,k}^2{\sum }_k^{\phantom{2}}\right)$

where ${\sum }_k$ is the conditional covariance matrix in (28). Applying (A5) to each term under the sum in (A7a), it follows that

(A7b)

$\mathbb{E}{\left|{\tilde{\eta }}_n^{\phantom{2}}\right|}^2\le d^{\frac{1}{\phantom{(}2}}\sum _{k=1}^n\frac{a^2}{k}\mathbb{E}{\left|A_{n,k}\right|}_{\mathrm{Op}}^2{\left|{\sum }_k^{\phantom{2}}\right|}_{\mathrm{F}}\le \left(d^{\frac{1}{\phantom{(}2}}{a}^2{C}_A^2\right)\frac{1}{n^{2\nu }}\sum _{k=1}^n{k}^{2\nu -1}\mathbb{E}{\left|{\sum }_k^{\phantom{2}}\right|}_{\mathrm{F}}$

where d is the dimension of the parameter space $\mathsf{\Theta }$, and ${\left|{\sum }_k^{\phantom{2}}\right|}_{\mathrm{F}}$ denotes the Frobenius matrix norm. However, it follows from (u1) that there exists a uniform upper bound S on ${\left|{\sum }_k^{\phantom{2}}\right|}_{\mathrm{F}}$. Therefore, by (A7b)

(A7c)

$\mathbb{E}{\left|{\tilde{\eta }}_n^{\phantom{2}}\right|}^2\le \left(d^{\frac{1}{\phantom{(}2}}{a}^2{C}_A^2\right)\frac{S}{n^{2\nu }}\sum _{k=1}^n{k}^{2\nu -1}$

Since $\nu >0$, the right-hand side of (A7c) remains bounded as $n\to \infty$, by the Euler-Maclaurin formula [37]. This immediately yields Condition (27b). □

Condition (27c): to verify this condition, it is first admitted that the following limit is known to hold

(A8a)

$lim\mathbb{E}\left({\sum }_k\right)={\sum }^{*}$

where ${\sum }^{*}$ was defined in (4a). Then, let the sum in (27c) be written

(A8b)

$\sum _{k=1}^n\frac{a^2}{k}{A}_{n,k}^{\phantom{2}}{\sum }_k^{\phantom{2}}{A}_{n,k}^{\phantom{2}}=\sum _{k=1}^n\frac{a^2}{k}{A}_{n,k}^{\phantom{2}}{\sum }^{*}{A}_{n,k}^{\phantom{2}}+\sum _{k=1}^n\frac{a^2}{k}{A}_{n,k}^{\phantom{2}}\left[{\sum }_k^{\phantom{2}}-{\sum }^{*}\right]A_{n,k}^{\phantom{2}}$

Due to the equivalence $A_{n,k}\sim exp(ln(n/k)A)$ (see [10], Page 125), the first term in (A8b) is a Riemann sum for the integral [10,34]

$a^2{\int }_0^1{e}^{-ln\left(s\right)A}{\sum }^{*}{e}^{-ln\left(s\right)A}dln\left(s\right)=a^2{\int }_0^{\infty }{e}^{-tA}{\sum }^{*}{e}^{-tA}dt$

which is known to be the solution $\sum$ of Lyapunov’s equation (7). The second term in (A8b) can be shown to converge to zero in probability, using inequality (A5) and the limit (A8a), by a similar argument to the ones in the verification of Conditions (27a) and (27b). Then, Condition (27c) follows immediately. □

Proof of (A8a): recall that $w_k=u_k+v_{k-1}$ where $u_k=u({\theta }_{k-1},x_k)$ and $v_{k-1}=\nabla D\left({\theta }_{k-1}\right)$. Since $\left(w_k\right)$ is a sequence of square-integrable martingale differences, it is possible to write, in the notation of (28),

(A9a)

${\sum }_k^{\phantom{2}}=\mathbb{E}\left[u_k^{\phantom{†}}{u}_k^{†}|{\mathcal{X}}_{k-1}\right]-v_{k-1}^{\phantom{†}}{v}_{k-1}^{†}$

By (18b), the second term in (A9a) converges to zero almost surely, as $k\to \infty$. It also converges to zero in expectation, since $\nabla D\left(\theta \right)$ is uniformly bounded for $\theta$ in the compact set ${\mathsf{\Theta }}^{*}$. For the first term in (A9a), since the $x_k$ are i.i.d. with distribution $P_{{\theta }^{*}}$, it follows that

(A9b)

$\mathbb{E}\left[u_k^{\phantom{†}}{u}_k^{†}|{\mathcal{X}}_{k-1}\right]=E_{{\theta }^{*}}\left[u({\theta }_{k-1},x)u^{†}({\theta }_{k-1},x)\right]$

Since $u(\theta ,x)$ is a continuous vector field on $\mathsf{\Theta }$ for each $x\in X$, and ${\theta }_{k-1}$ converge to ${\theta }^{*}$ almost surely, it follows that $u({\theta }_{k-1},x)$ converge to $u({\theta }^{*},x)$ for each $x\in X$, almost surely. On the other hand, it follows from (u2) that the functions under the expectation $E_{{\theta }^{*}}$ in (A9b) have bounded second order moments, so they are uniformly integrable [20]. Therefore,

(A9c)

$lim{E}_{{\theta }^{*}}\left[u({\theta }_{k-1},x)u^{†}({\theta }_{k-1},x)\right]=E_{{\theta }^{*}}\left[u({\theta }^{*},x)u^{†}({\theta }^{*},x)\right]={\sum }^{*}$

almost surely, by the definition (4a) of ${\sum }^{*}$. It now follows from (A9a), (A9b), and (A9c) that the following limit holds

(A9d)

$lim{\sum }_k={\sum }^{*}\mathrm{almost}\mathrm{surely}$

To obtain (A8a) it is enough to note, as already stated in the verification of Condition (27b), that the ${\sum }_k$ are uniformly bounded in the Frobenius matrix norm. Thus, (A9d) implies (A8a), by the dominated convergence theorem. □

Appendix C. Background on the Information Metric

Let $D\left(\theta \right)$ be the Kullback–Leibler divergence (2a) or any other so-called $\alpha$-divergence [19]. Assume the Riemannian metric $⟨·,·⟩$ of $\mathsf{\Theta }$ coincides with the information metric of the model P. Then, for any local coordinates $({\tau }^{\alpha };\alpha =1,\dots ,d)$, with origin at ${\theta }^{*}$, the following relation holds, by definition of the information metric (see [19], Page 54),

(A10)

${\left.\frac{{\partial }^{2}D}{\phantom{(}\partial {\tau }^{\alpha }\partial {\tau }^{\beta }}\right|}_{{\tau }^{\alpha }=0}={⟨\frac{\partial }{\partial {\tau }^{\alpha }},\frac{\partial }{\partial {\tau }^{\beta }}⟩}_{{\theta }^{*}}$

where $\frac{\partial }{\phantom{(}\partial {\tau }^{\alpha }}$ denote the coordinate vector fields of the local coordinates ${\tau }^{\alpha }$. It is also possible to express (A10) in terms of the Riemannian distance $d(·,·)$, induced by the information metric $⟨·,·⟩$. Precisely,

(A11)

$D\left(\theta \right)=\frac{1}{2}{d}^2(\theta ,{\theta }^{*})+o\left(d^2(\theta ,{\theta }^{*})\right)$

This follows immediately from the second-order Taylor expansion of $D\left(\theta \right)$, since ${\theta }^{*}$ is a minimum of $D\left(\theta \right)$, by using (A10). Formula (A11) shows that the divergence $D\left(\theta \right)$ is equivalent to half the squared Riemannian distance $d^2(\theta ,{\theta }^{*})$, at $\theta ={\theta }^{*}$.

The scalar products appearing in (A10) form the components of the information matrix $I^{\tau }$ of the coordinates ${\tau }^{\alpha }$,

$I_{\alpha \beta }^{\tau }={\left.\frac{{\partial }^{2}D}{\phantom{(}\partial {\tau }^{\alpha }\partial {\tau }^{\beta }}\right|}_{{\tau }^{\alpha }=0}$

In any change of coordinates, these transform like the components of a $(0,2)$ covariant tensor [16]. That is, if $({\theta }^{\alpha };\alpha =1,\dots ,d)$ are any local coordinates defined at ${\theta }^{*}$,

$I_{\alpha \beta }^{\tau }={\left(\frac{\partial {\theta }^{\gamma }}{\partial {\tau }^{\alpha }}\right)}_{{\theta }^{*}}{I}_{\gamma \kappa }^{\theta }{\left(\frac{\partial {\theta }^{\kappa }}{\partial {\tau }^{\beta }}\right)}_{{\theta }^{*}}$

where the subscript ${\theta }^{*}$ indicates the derivative is evaluated at ${\theta }^{*}$, and where $I_{\gamma \kappa }^{\theta }$ are the components of the information matrix $I^{\theta }$ of the coordinates ${\theta }^{\alpha }$.

The recursive estimates ${\theta }_n$ are said to be asymptotically efficient, if they are asymptotically efficient in any local coordinates ${\tau }^{\alpha }$, with origin at ${\theta }^{*}$. That is, according to the classical definition of asymptotic efficiency [13,14], if the following weak limit of probability distributions is verified [20],

(A12)

$\mathcal{L}\left\{\left(n^{1/2}{\tau }_n^{\alpha }\right)\right\}⟹N_d\left(0,{\sum }^{\tau }\right){\sum }^{\tau }={\left(I^{\tau }\right)}^{-1}$

where $\mathcal{L}\{\dots \}$ denotes the probability distribution of the quantity in braces, ${\tau }_n^{\alpha }={\tau }^{\alpha }\left({\theta }_n\right)$ are the coordinates of the recursive estimates ${\theta }_n$, and $N_d\left(0,{\sum }^{\tau }\right)$ denotes a centred d-variate normal distribution with covariance matrix ${\sum }^{\tau }$.

It is important to note that asymptotic efficiency of the recursive estimates ${\theta }_n$ is an intrinsic geometric property, which does not depend on the particular choice of local coordinates ${\tau }^{\alpha }$, with origin at ${\theta }^{*}$. This can be seen from the transformation rule of the components of the information matrix, described above. In fact, since these transform like the components of a $(0,2)$ covariant tensor, the components of ${\sum }^{\tau }$ transform like those of a $(2,0)$ contravariant tensor, which is the correct transformation rule for the components of a covariance matrix.

Figures

Figure 1. Fast non-asymptotic rate of convergence

Figure 2. Asymptotic efficiency (optimal rate of convergence)