Salem Said
email: salem.said@u-bordeaux.fr

Nicolas Le Bihan
email: nicolas.le-bihan@gipsa-lab.grenoble-inp.fr

Jonathan Manton
email: jmanton@unimelb.edu.au

Hidden Markov chains and fields with observations in Riemannian manifolds

Keywords: Riemannian manifold, hidden Markov model, Markov field, EM algorithm

Hidden Markov chain, or Markov field, models, with observations in a Euclidean space, play a major role across signal and image processing. The present work provides a statistical framework which can be used to extend these models, along with related, popular algorithms (such as the Baum-Welch algorithm), to the case where the observations lie in a Riemannian manifold. It is motivated by the potential use of hidden Markov chains and fields, with observations in Riemannian manifolds, as models for complex signals and images.

Hidden Markov chains and fields with observations in Riemannian manifolds

Salem Said * Nicolas Le Bihan * * Jonathan H. Manton * * *

INTRODUCTION

The present work is concerned with hidden Markov chain and Markov field models, with observations in Riemannian manifolds.

It introduces a general formulation of hidden Markov chain models, whose observations lie in a Riemannian manifold, and derives an expectation-maximisation algorithm, to estimate the parameters of these models -Sections 2 and 3, respectively.

It also describes a general hidden Markov field model, with observations in a Riemannian manifold, and discusses briefly the estimation of this model, using the expectationmaximisation approach -Section 5.

Hidden Markov chain and field models have been highly influential in signal and image analysis [START_REF] Elliot | Hidden Markov Models : Estimation and Control[END_REF] [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. For example, the Baum-Welch algorithm, for estimating hidden Markov chain models, is a staple of applications such as speech recognition and protein sequencing [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF] [START_REF] Durbin | Biological sequence analysis[END_REF]. Also, hidden Markov fields are central in several approaches to image restoration and segmentation [START_REF]Markov Random Field Modeling in Image Analysis[END_REF].

While quite extensive, existing literature is almost entirely focused on hidden Markov chain and field models with observations in a Euclidean space. The present work aims to provide a statistical framework, which may be used to extend current algorithms, dealing with observations in a Euclidean space, to observations in more general Riemannian manifolds.

For example, the expectation-maximisation algorithm, derived in Section 3, extends the popular Baum-Welch algorithm to hidden Markov chains with observations in a Riemannian manifold.

The motivation for considering hidden Markov chains and fields with observations in Riemannian manifolds comes from the important role which data in Riemannian manifolds can play in the study of complex signals and images (or even three-dimensional fields, such as turbulent flows [START_REF] Zare | Colour of Turbulence[END_REF]).

To understand this, consider the case of a non-stationary, multivariate signal, observed over a long time interval. This is typical of a phased-array radar, or brain-computer interface signal [START_REF] Jeuris | The Kähler mean of block Toeplitz matrices with Toeplitz structured blocks[END_REF] [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF].

While this signal cannot be entirely described by a single tractable model, it is still possible to "zoom in" on short time windows, which admit of a complete, tractable description, very often in terms of a covariance matrix (for array radar signals, this has an additional block-Toeplitz structure [START_REF] Jeuris | The Kähler mean of block Toeplitz matrices with Toeplitz structured blocks[END_REF]). Now, the original complex signal appears as a time series of covariance matrices, or other descriptors, one for each short time window, and it seems possible to extract useful information, by studying the statistical distribution of these descriptors.

This idea can also be applied to a complex image, or three-dimensional field, by zooming in on small connected regions [START_REF] Tuzel | Pedestrian detection via classification on Riemannian manifolds[END_REF] [START_REF] Zare | Colour of Turbulence[END_REF]. Recently, it has generated significant interest in the problem of estimating an unknown distribution, underlying a population of descriptors (e.g., covariance matrices, principal subspaces, histogrammes) [START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF][11] [START_REF] Banerjee | Clustering on the unit hypersphere using von Mises-Fisher distributions[END_REF].

Mathematically, this problem was framed as the problem of estimating an unknown probability distribution on a Riemannian manifold. Indeed, in several applications, this Riemannian approach has become standard, due to its superior performance (e.g., brain-computer interface [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF]).

Existing work on this problem has always started from the assumption that descriptors, obtained from different time windows of a signal (or regions of an image), are sampled independently, all from the same underlying distribution. This assumption obscures the dynamics of the signal (or the spatial dependence structure within the image), potentially leading to the loss of crucial information.

The statistical framework introduced in the present work captures this information, in the form of a hidden Markov dependence structure, which controls the distribution of the descriptors.

HIDDEN MARKOV CHAINS

Hidden Markov chains are popular models for many kinds of signals (neuro-biological, speech, etc. [1][2]). Here, these models are generalised to manifold-valued signals.

In a hidden Markov chain model, one is interested in understanding some hidden, time-varying, finitely-valued, process, say (q t ; t = 1, 2, . . .) which takes values in a finite set S. When q t = α for some α ∈ S, one says that the process is in state α at time t.

Moreover, it is assumed that (q t ) is a time-stationary Markov chain, so there exists a so-called transition matrix, (P αβ ; α, β ∈ S), which gives the conditional probabilities

P(q t+1 = β|q t = α) = P αβ
(1) In particular, if π t (α) = P(q t = α) is the distribution of q t , then [START_REF] Norris | Markov Chains[END_REF] π

t+1 (β) = α∈S π t (α)P αβ (2)
describes the transition from time t to time t + 1.

The states (q t ) are hidden, so they can only be observed through random outputs (y t ) which have their values in a Riemannian manifold M .

Moreover, these outputs are generated independently from each other, and each y t depends only on

q t [1][2].
It is assumed M is a homogeneous Riemannian manifold, and y t is distributed according to a location-scale model, p(y

t |q t = α) = f (y t |ȳ α , σ α ) (3)
where p(y t |q t = α) denotes the conditional density, with respect to the Riemannian volume measure of M .

The location parameters ȳα ∈ M , and scale parameters σ α ≥ 0 are unknown, but the function f (y t |ȳ α , σ α ) is known, and of the form [START_REF] Said | Warped Riemannian metrics for location-scale models[END_REF] f (y|ȳ, σ) = exp [ η(σ)D(y , ȳ) -ψ(η(σ))] (4) where D : M × M → R, and where η(σ) is the so-called natural parameter, while ψ(η) is a strictly convex function.

Examples of location-scale models of the form (4) include the von Mises-Fisher model, and the Riemannian Gaussian model, which are briefly recalled, in Remark 1, at the end of the present section.

From the definition of π t and from (3), it follows that the probability density of y t is a time-varying mixture density

p(y t ) = α∈S π t (α) f (y t |ȳ α , σ α ) (5)
Having access only to the observations (y t ; t = 1, 2, . . .), one hopes to recover as much information about the Markov chain (q t ) as possible : its transition matrix (P αβ ) and the parameters (ȳ α , σ α ) which govern its outputs y t .

The problem of estimating (P αβ ) and (ȳ α , σ α ) is addressed in Section 3. A further, equally interesting problem, will not be considered, for lack of space. This is the problem of estimating the invariant distribution of the chain (q t ) [START_REF] Norris | Markov Chains[END_REF] : if the transition matrix (P αβ ) is irreducible and aperiodic, then π t (α) → π(α), the invariant distribution, as t → ∞. This means that, in time, the observation y t will tend to sample from a mixture distribution of the same form ( 5), but with the invariant weights π(α) instead of π t (α).

Remark 1 : location-scale models of the form (4) were introduced in [START_REF] Said | Warped Riemannian metrics for location-scale models[END_REF]. For these models, M is any homogeneous Riemannian manifold [START_REF] Besse | Einstein manifolds[END_REF] : a Riemannian manifold with a group of isometries G which acts transitively.

Isometry means each g ∈ G is a mapping g : M → M which preserves the Riemannian metric of M . Moreover, transitive action means that for each y, z ∈ M there exists g ∈ G such that g • y = z (here, g • y = g(y)).

The function D : M ×M → R, appearing in ( 4) is required to satisfy the group-invariance property

D(y , ȳ) = D(g • y , g • ȳ) (6) for any g ∈ G and y, ȳ ∈ M .
This property guarantees that f (y|ȳ, σ) is indeed a true probability density. In particular, the function ψ(η) is strictly convex, because it is the cumulant generating function of the statistic ∆ = D(y , ȳ),

E (ȳ,σ) exp(η ∆) = e ψ(η)
where the expectation is with respect to the density f (y|ȳ, σ) -this identity holds for any real η as long as the expectation is finite.

Examples of location-scale models of the form (4) include the von Mises-Fisher model, and the Riemannian Gaussian model.

von Mises-Fisher model [START_REF] Mardia | Directional statistics[END_REF] : for this model, Here, G acts on M by g • y = gyg † where † denotes transpose. Moreover, this action preserves the Riemannian distance,

M = S d-1 and G = O(d) (7)
d 2 (y, z) = tr log(y -1 z) 2
where tr denotes the trace, and log the symmetric matrix logarithm. The function

D : M × M → R is then given by D(y , ȳ) = d 2 (y, ȳ) (11
) and an expression of ψ(η), in terms of a multivariate integral, is given in [START_REF] Said | Riemannian Gaussian distributions in the space of symmetric positive definite matrices[END_REF].

The model is obtained by setting η(σ) = -1 2σ 2 in (4). This gives the Riemannian Gaussian density

f (y|ȳ, σ) = Z -1 (σ) exp - d 2 (y, ȳ) 2σ 2 (12)
where Z(σ) = e ψ(η(σ)) .

THE EM ALGORITHM

Here, an EM (expectation-maximisation) algorithm will be introduced, which addresses the estimation problem defined in Section 2.

To see why the expectation-maximisation approach is a natural choice, assume it were possible to access the socalled complete data x t = (y t , q t ) at times t = 1, . . . , T . The resulting log-likelihood function would be,

(x|θ) = log π 1 (q 1 )
T -1

t=1
Pq t q t+1 + log

T t=1 f (y t |ȳq t , σq t ) ( 13
)
where θ = (P αβ , ȳα , σ α ). However, this could also be written1 ,

(x|θ) = α,β∈S N αβ (q) log P αβ + α∈S T t=1 δα(q t ) log f (y t |ȳα , σα) (14)
where N αβ (q) is the number of times that (q t , q t+1 ) = (α, β), and where δ α (q t ) = 1 if q t = α and = 0 otherwise.

Estimating θ from the complete data, by maximising the log-likelihood function [START_REF] Said | Warped Riemannian metrics for location-scale models[END_REF], is relatively straightforward, as explained in Remark 3 below.

Unfortunately, since the Markov chain (q t ) is hidden, one has access only to the observations (y t ; t = 1, . . . , T ).

The log-likelihood function for these observations is much more complicated than ( 14),

(y|θ) = log   q 1 ∈S . . . q T ∈S π 1 (q 1 ) T -1 t=1 Pq t q t+1 T t=1 f (y t |ȳq t , σq t )   (15)
Instead of directly maximising [START_REF] Besse | Einstein manifolds[END_REF], the expectationmaximisation approach reduces the problem of estimating θ from the observations y t to a sequence of steps which only involve maximising a function of the form [START_REF] Said | Warped Riemannian metrics for location-scale models[END_REF].

Specifically, the k-th iteration of the EM algorithm computes estimates

θ k = ( P k αβ , ȳk α , σ k α )
which are guaranteed to increase the log-likelihood [START_REF] Besse | Einstein manifolds[END_REF], in the sense that (y|θ k ) ≥ (y|θ k-1 ). Each iteration is made up of two steps [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]

[2] • E step : compute the function Q(θ|θ k-1 ) = E θ k-1 [ (x|θ)|y] ( 16
) where E θ k-1 means the expectation is computed under the current estimate θ k-1 . This is called the expectation step.

• M step : maximise Q(θ|θ k-1 ) to obtain the new estimate

θ k = argmax θ Q(θ|θ k-1 ) (17)
In principle, the algorithm cycles through these two steps, until the estimates θ k converge. In practice, this means the algorithm is designed to stop when the difference between θ k and θ k-1 falls short of some pre-assigned accuracy.

There is no theoretical guarantee that the θ k will converge to the global maximum of the log-likelihood [START_REF] Besse | Einstein manifolds[END_REF]. In fact, convergence typically takes place only to a local maximum, depending on initialisation.

However, the EM algorithm remains quite advantageous, from a computational point of view, since it boils down to relatively simple operations, with little computational complexity, in comparison with alternative methods. Indeed, replacing ( 14) into ( 16), yields Q(θ|θ k-1 ) = Q(θ),

Q(θ) = α,β∈S ν αβ (y) log P αβ + α∈S T t=1 ω t α (y) log f (y t |ȳα , σα) (18)
with ν αβ (y) and ω t α (y) computed from the observations y t ,

ν αβ (y) = T -1 t=1 P θ k-1 (q t = α , q t+1 = β|y) (19a) ω t α (y) = P θ k-1 (q t = α|y) (19b
) where P θ k-1 means the probability is computed under the current estimate θ k-1 .

Thus, the E step ( 16) amounts to implementing formulae (19a) and (19b). This is efficiently realised using Levinson's forward-Backward algorithm, recalled in Remark 2.

On the other hand, the M step [START_REF] Said | Riemannian Gaussian distributions in the space of symmetric positive definite matrices[END_REF] amounts to maximising the function [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], which is of the form ( 14), and can be maximised as explained in Remark 3.

Remark 2 : in order to implement formulae [START_REF] Devijver | Baum's forward-backward algorithm revisited[END_REF], it is helpful to use the forward variables Φ t (α) and backward variables B t (α), which satisfy

ν αβ (y) = T -1 t=1 Φ t (α) (P αβ f (y t+1 |ȳ β , σ β )) B t+1 (β) (20a) ω t α (y) = Φ t (α) B t (α) (20b)
and can be computed using Levinson's forward-backward algorithm, in terms of the forward recursion

Φ t+1 (β) = Φ -1 t α∈S Φ t (α)P αβ f (y t+1 |ȳ β , σ β ) (21a)
and of the backward recursion,

B t (α) = Φ -1 t+1 β∈S P αβ B t+1 (β) f (y t+1 |ȳ β , σ β ) (21b)
where Φ t is a normalising factor,

Φ t = β∈S α∈S Φ t (α)P αβ f (y t+1 |ȳ β , σ β ) (21c)
and where (21a) is taken with the initial condition

Φ 1 (α) = π 1 (α)f (y 1 |ȳ α , σ α ) α∈S π 1 (α)f (y 1 |ȳ α , σ α ) (21d)
and (21b) is taken with the terminal condition B T (β) = 1.

The recursions (21a) and (21b) are efficient alternatives to Baum's forward-backward algorithm, introduced in [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF].

They have the same computational complexity (roughly, the order of T |S| 2 where |S| is the number of elements in S). On the other hand, they do not suffer from the same numerical instability, which makes Baum's algorithm rather difficult to use [START_REF] Devijver | Baum's forward-backward algorithm revisited[END_REF].

The proof of (20) and of ( 21) is discussed in Appendix A.

Remark 3 : consider the task of maximising a function of the form [START_REF] Said | Warped Riemannian metrics for location-scale models[END_REF],

Q(θ) = α,β∈S ν αβ log P αβ + α∈S T t=1 ω t α log f (y t |ȳα , σα) (22)
with ν αβ , ω t α > 0 and f (y|ȳ, σ) given by ( 4), where the natural parameter η(σ) is negative -if η(σ) is positive, as in the von Mises-Fisher model, it is possible to absorb a minus sign into the function D.

Then, any global maximiser of [START_REF] Grimmett | Probability on graphs : stochastic processes on graphs and lattices[END_REF], say θ = ( Pαβ , ŷα , σα ), must satisfy Pαβ = Pαβ (ν) ; ŷα = ŷα (ω) ; η(σ α ) = ηα (ω)

where the following definitions are used

Pαβ (ν) = ν αβ β∈S ν αβ (23a) ŷα (ω) = argmin y∈M T t=1 ω t α D(y t , y) (23b) ηα (ω) = (ψ ) -1 T t=1 ω t α D(y t , ŷα ) T t=1 ω t α (23c)
Here, ψ is the derivative of ψ and (ψ ) -1 its reciprocal function. In particular, (ψ ) -1 is well-defined since ψ is strictly convex.

The proof of Formulae ( 23) is discussed in Appendix B.

For the von Mises-Fisher model, the solution of the minimisation problem (23b) is given by [START_REF] Mardia | Directional statistics[END_REF] ŷα (ω) = U For the Riemannian Gaussian model, (23b) reads, after replacing [START_REF] Chevallier | Kernel density estimation on the Siegel space with application to Radar processing[END_REF],

ŷα (ω) = argmin y∈M T t=1 ω t α d 2 (y t , y) (24b)
This means that ŷα (ω) is the weighted Riemannian centre of mass of the observations y t , and can be computed using any of the standard methods for computing Riemannian centres of mass [START_REF] Said | Riemannian Gaussian distributions in the space of symmetric positive definite matrices[END_REF][10].

Formula (23c) is easily evaluated. The function ψ (η) being known, ηα (ω) is the unique solution of the equation

ψ (η α (ω)) = T t=1 ω t α D(y t , ŷα ) T t=1 ω t α (25)
Which can be found, to any accuracy, in exponential time, using the Newton-Raphson method, since the function ψ(η) is strictly convex.

It is also possible to find ηα (ω) directly, by performing a search in a table which contains the values of the function ψ (η). If the table is sufficiently large, this approach is successful within constant time, but has limited accuracy.

COMPUTER EXPERIMENT

Here, a numerical experiment is considered, illustrating the hidden Markov chain model of Section 2, and testing the validity of the EM algorithm of Section 3.

In this experiment, the hidden Markov chain (q t ) takes values in a three-element set S = {1, 2, 3}, and has initial distribution π 0 (α) = δ 1 (α), where δ 1 (α) = 1 if α = 1 and = 0 otherwise, and transition matrix (P αβ ),

(P αβ ) = 0.4 0.3 0.3 0.2 0.6 0.2 0.1 0.1 0.8 (26)
The outputs (y t ) were generated from a Riemannian Gaussian model in the Poincaré disk. In other words, each (y t ) has its values in M = {z ∈ C : |z| < 1}, and the conditional density (3) takes on the form [START_REF] Banerjee | Clustering on the unit hypersphere using von Mises-Fisher distributions[END_REF], where [START_REF] Vinberg | Geometry II : Spaces of constant curvature[END_REF] d(y, z)

= acosh 1 + 2|y-z| 2 (1-|y| 2 )(1-|z| 2 ) (27a) Z(σ) = (2π) 3/2 σ e σ 2 2 erf σ √ 2 (27b)
erf being the error function [START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF].

The Poincaré disk is considered here, rather than a space P d of positive-definite matrices, as in Section 2, since it is much easier to visualise (indeed, it is just the inside of a unit disk in the plane). On the other hand, the Poincaré disk, with Riemannian distance (27a) is isometric to the space of 2 × 2 symmetric positive-definite matrices, with unit determinant [START_REF] Vinberg | Geometry II : Spaces of constant curvature[END_REF].

Some background on the Riemannian geometry of the Poincaré disk is recalled in Remark 4.

The following sets of location parameters ȳα and scale parameters σ α were used ȳ1 = 0 ; ȳ2 = 0.82 i + 0.29 ; ȳ3 = 0.82 i -0.29 (28a)

σ 1 = 0.1 ; σ 2 = 0.4 ; σ 3 = 0.4 (28b) where i = -1.
With the hidden Markov model defined in this way, the observations (y t ; t = 1, . . . , T ) where T = 10000 appeared as a "face" pattern, within the Poincaré disk, shown in Figure 1.

If the scale parameters σ α are all multiplied by 5, the face becomes more blurred, since each conditional density (3) has its dispersion multiplied by 5. This can be seen in Figure 2.

In both figures 1 and 2, the "eyes" (centred at ȳ2 and ȳ3 ) appear smaller than the "mouth", although are they have scale parameters σ 2 and σ 3 four times larger (see (28b) !).

This apparent contradiction is due to the fact that the Riemannian metric of the Poincaré disk is stretched near the boundary of the disk -see Formula (29), Remark 4. Thus, while they appear smaller to us, objects near the boundary are, in fact, bigger when measured by this Riemannian metric.

Consider now the application of the EM algorithm of Section 3, to the estimation of the transition matrix (26) and location and scale parameters (28). Remark 4 : the Poincaré disk is the set of complex numbers z with |z| < 1, equipped with a conformal Riemannian metric, which turns it into a space of constant negative curvature [START_REF] Vinberg | Geometry II : Spaces of constant curvature[END_REF]. This Riemannian metric measures the length of a plane vector w (identified with a complex number), attached at the point z, using the Riemannian norm

w 2 z = 4 κ |w| 2 (1 -|z| 2 ) 2 (29)
where κ > 0 is a scaling factor, and |w| 2 = ww * is the squared modulus, or squared Euclidean norm (the star denotes the conjugate).

It is called a conformal metric, since the Riemannian norm w z is proportional to the Euclidean norm |w|. The sectional curvature associated to this metric is constant, and equal to -κ.

In (29), the conformal factor (1-|z| 2 ) -2 goes to infinity as |z| approaches 1. Accordingly, the Riemannian norm w z of a vector w, with fixed Euclidean norm |w|, becomes arbitrarily large, when z approaches the boundary of the disc.

In other words, the Riemannian metric of the Poincaré disk is stretched near the boundary.

This algorithm was run N mc = 20 times, each time for N em = 1000 iterations. An individual run took about one hour, on a standard laptop.

All of these runs were carried out with the same initial guess θ 1 , but the observations (y t ; t = 1, . . . , T ) were generated anew for each run, leading to a different (in fact, value of the final estimate

θ Nem = ( P Nem αβ , ȳNem α , σ Nem α )
Then, the mean value of ( P Nem αβ ) was found to be Here, E mc is the "Monte Carlo expectation", equal to the empirical average, taken over N mc runs.

The dispersion of ( P Nem αβ ), about its mean value, was measured by the maximum empirical variance of the entries P Nem αβ . This was found to be max Building on this remark, it is possible to infer that the discrepancy, between these mean values and the true values of the transition probabilities P αβ , and location and scale parameters ȳα and σ α , is due to the insufficient number of observations, T = 10000, rather than to a slow convergence of the EM algorithm. This is supported by the fact that each run of the EM algorithm was observed to converge to the final estimate θ Nem well before N em = 1000 iterations, producing nearly constant estimates θ k ≈ θ Nem not later than k = 800.

When the number of observations was increased to T = 50000, an individual run, of N em = 1000 iterations of the algorithm, took much longer, over five hours.

For this reason, only one run of the algorithm was carried out. This was found to produce final estimates significantly closer to the true values of the transition matrix and location and scale parameters. However, for lack of time, a more detailed statistical analysis was not possible, in this case.

HIDDEN MARKOV FIELDS

Similar to hidden Markov chain models, hidden Markov field models involve a hidden stochastic process2 (q z ), which takes its values in a finite set S, observed through random outputs (y z ), which belong to a homogeneous Riemannian manifold M . However, in the case of a hidden Markov field, the index z is not a time variable (running through t = 1, 2, . . .), but a vertex of a finite unidrected graph Z. Then, the stochastic process (q z ; z ∈ Z) is called a Markov field if the following two properties are satisfied [START_REF] Grimmett | Probability on graphs : stochastic processes on graphs and lattices[END_REF] P(q z = α|q w ; w = z) = P(q z = α|q w ; w ∼ z)

where w ∼ z means w and z are adjacent, and

P (∩ z∈Z {q z = ς(z)}) > 0 (30b)
for any function ς : Z → S. Such a function is called a configuration of the field.

Property (30a) says that the state q z at a vertex z depends on states at the other vertices q w only through those w which are adjacent to z.

On the other hand, Property (30b) is required for the Hammersley-Clifford theorem [START_REF] Grimmett | Probability on graphs : stochastic processes on graphs and lattices[END_REF]. This theorem implies that the probability distribution of a Markov field is a Gibbs distribution.

For simplicity, it is here assumed that the graph Z is a square grid, where each vertex is adjacent to its immediate neighbors (this is a simple model of an image). In this case, the field (q z ) is said to have a Gibbs distribution, if the probability mass function of the states q z takes on the form p(q z ; z ∈ Z) ∝ z∈Z exp V (q z ) + 1 2 w∼z J(q z , q w ) (31)

Here, ∝ indicates a missing normalising factor W (V,J) .

Computing W (V,J) requires summing over all possible configurations of the field q -often, an impracticable task.

The observed outputs (y z ) depend on the hidden states (q z ) exactly as in Section 2. That is,

p(y z |q z = α) = f (y z |ȳ α , σ α ) (32)
where the function f (y z |ȳ α , σ α ) is of the form (4). Then, having access only to these outputs (y z ), one hopes to recover the dependence structure of the Hidden field, given by the functions V : S → R and J : S × S → R, as well as the parameters (ȳ α , σ α ) which govern the (y z ).

This problem generalises the hidden Markov chain model problem of Section 2, which was addressed in Section 3.

It is a problem of parameter estimation, for the unknown parameter θ = (V α , J αβ , ȳα , σ α ), where V α = V (α) and J αβ = J(α, β) for α, β ∈ S.

In the special case where the outputs y z belong to a Euclidean space M = R d , several EM algorithms were introduced in the image analysis literature, in order to address this problem [23][24].

The extension of these algorithms, to observed outputs y z in any homogeneous Riemannian manifold M , is here discussed briefly.

Each iteration of the EM algorithm includes an E step and an M step. Here, these take on the following form, where the notation is the same as in ( 16) and ( 17) of Section 3unfortunately, due to lack of space, the following discussion does not include any proofs.

• E step : compute the function Moreover, the maximum in (36a) is unique, since Ψ(V, J) is a strictly convex function of V = (V α ) and J = (J αβ ), (when these are considered as V ∈ R |S| and J ∈ R |S|×|S| ).

Q(θ|θ k-1 ) = Q(ȳ α , σ α ) + Q(V, J) (33
Remark 5 : in its form specified by (34) and (36), the EM algorithm is not practically applicable. Indeed, both the normalising constant W (V, J), and the conditional probabilities and expectations (35), are impossible to evaluate, outside of trivial situations.

To circumvent this (fundamental) difficulty, various approximate methods, for evaluating W (V, J) and (35), have been proposed, such as the ones based on mean-field approximations [START_REF] Celeux | EM procedures using mean-field-like approximations for Markov-modelbased image segmentation[END_REF] [START_REF] Celeux | EM-based image segmentation using Potts models with external field[END_REF].

It will be the aim of the authors' future work to adapt and apply these methods, to the general case where the observed outputs y z belong to a homogeneous Riemannian manifold.

  where S d-1 ⊂ R d is the (d -1)-dimensional unit sphere, and O(d) is the group of orthogonal transformation of R d . The function D : M × M → R is given by D(y , ȳ) = y, ȳ R d (8) where •, • R d is the Euclidean scalar product in R d , and ψ(η) has the following expression (where ν = d/2) e ψ(η) = (2π) ν η 1-ν I ν-1 (η) (9) where I ν-1 is the modified Bessel function of order ν -1. Finally, the von Mises-Fisher model is obtained by setting η(σ) = σ in (4) -then, η is called the concentration parameter. Riemannian Gaussian model [17][10] : for this model M = P d and G = GL(d) (10) where P d is the space of d × d symmetric positive-definite matrices, and GL(d) the group of invertible d×d matrices.

  with U : R d -{0} → R d the mapping U (y) = y/ y R d where • R d denotes the Euclidean norm.

Fig. 1 .Fig. 2 .

 12 Fig. 1. 10 4 observations ; location-scale parameters (28)

1 ) 1 ) = 0. 5 E mc (ȳ Nem 2 ) 2 ) = 0. 3 E

 115223 = -0.02 i + 0.41 ; E mc (σ Nem = 0.83 i + 0.30 ; E mc (σ Nem mc (ȳ Nem 3 ) = 0.83 i -0.29 ; E mc (σ Nem3 ) = 0.3

1 where

 1 V mc is the "Monte Carlo variance", equal to the empirical variance, taken over N mc runs.These numerical results show that the value of the final estimate θ Nem is rather stable, and is not affected very strongly by the randomness of the observations. Indeed, the dispersions of P Nem αβ or of ȳNem α and σ Nem α appear small in comparison to their mean values.

  ) in terms of the following formulaeQ(ȳ α , σ α ) = α∈S z∈Z ω z α (y) log f (y z |ȳ α , σ α ) (34a) Q(V, J) = -log W (V, J) + α∈S ω α (y)V α + 1 2 α,β∈S ν αβ (y)J αβ (34b)which involve conditional probabilities and expectations,ω z α (y) = P θ k-1 (q z = α|y) ; ω α (y) = z∈Z ω z α (y) (35a) ν αβ (y) = w∼z P θ k-1 (q z = α , q w = β|y)(35b)• M step : maximise Q(θ|θ k-1 ) to obtain new estimates, θ k = ( Vα , Ĵαβ , ŷα , σα ). These are given by,( V , Ĵ) = argmax (V,J) ω, V + ν, J -Ψ(V, J) (36a)and, in the notation of Remark 3, ŷα = argmin ŷ∈M z∈Z ω z α (y) D(y z , ŷ) (36b)η(σ α ) = (ψ ) -1 z∈Z ω z α (y) D(y z , ŷα ) z∈Z ω z α (y)(36c)Here, Ψ(V, J) = log W (V, J), and ω, V = α∈S ω α (y)V α ; ν, J = 1 2α,β∈S ν αβ (y)J αβ

the term depending on π1 is discarded, since it is statistically nonsignificant when T is large.

a stochastic process is any family of jointly defined random variables[START_REF] Kallenberg | Foundations of modern probability[END_REF].

Appendix A. FORWARD-BACKWARD VARIABLES

Define the forward variables Φ t (α) and backward variables B t (α) by [START_REF] Devijver | Baum's forward-backward algorithm revisited[END_REF] -for convenience, write P instead of P θ k-1 Φ t (α) = P (q t = α|y t , . . . , y 1 ) (A.1) B t (α) = P (y T . . . , y t+1 |q t = β) p(y T . . . , y t+1 |y t , . . . , y 1 ) (A.2)

To see that these satisfy [START_REF] Vinberg | Geometry II : Spaces of constant curvature[END_REF], the key is to make use of the fact that x t = (y t , q t ) is a Markov chain, with transition probabilities [START_REF] Elliot | Hidden Markov Models : Estimation and Control[END_REF][2] p(y t+1 , q t+1 = β|y t , q t = α)

For example, (20a) can be obtained by using Bayes rule to express each term in the sum (19a)

where the numerator is the joint probability density of y T , . . . , y 1 on the event (q t , q t+1 ) = (α, β) -this is a density with respect to the Riemannian volume measure of the product manifold

This joint density can be written as a product, p(y T , . . . , y t+2 | y t+1 , q t+1 = β , q t = α , y t , . . . , y 1 ) × p(y t+1 , q t+1 = β | q t = α , y t , . . . , y 1 ) × P(q t = α| y t , . . . , y 1 ) × p(y t , . . . , y 1 ) by repeatedly applying (again!) Bayes rule. However, by the Markov property of x t = (y t , q t ), this product becomes p(y T , . . . , y t+2 | y t+1 , q t+1 = β) × p(y t+1 , q t+1 = β | q t = α) × P(q t = α| y t , . . . , y 1 ) × p(y t , . . . , y 1 ) but, using (A.1), (A.2) and (A.3), this is the same as

Therefore, replacing back into (A.4), it follows each term in the sum (19a) is given by

Thus, summing over t = 1 , . . . , T -1 yields (20a).

An analogous, and simpler, reasoning can be used to obtain (20b).

In addition to [START_REF] Vinberg | Geometry II : Spaces of constant curvature[END_REF], it should be proved that the forward and backward variables are given by [START_REF] Kallenberg | Foundations of modern probability[END_REF].

To obtain (21a), note first that the initial condition (21d) follows directly from (A.1). In addition, also from (A.1) Φ t+1 (β) ∝ P (q t+1 = β|y t+1 , . . . , y 1 )

∝ P (y t+1 , q t+1 = β|y t , . . . , y 1 )

where ∝ indicates a missing normalising factor. Thus, using the Markov property of x t = (y t , q t ), it is seen that, up to normalisation, Φ t+1 (β) is equal to α∈S p(y t+1 ; q t+1 = β|q t = α) × P(q t = α|y t , . . . , y 1 )

This immediately gives (21a), after noting that the missing normalising factor can be found by summing over β ∈ S, as in (21c).

The proof of (21b) is similar to that of (21a), and is here omitted to avoid repetition.

Appendix B. PROOF OF FORMULAE (23)

Recall the function Q(θ) of ( 22),

Proof of (23a) : only the first term in the expression of Q(θ) depends on (P αβ ). Therefore, this term can be maximised separately.

Denote this first term by Q(P ). Then, maximising Q(P ) is equivalent to maximising

where Pαβ (ν) is given by (23a).

Let Pαβ = Pαβ (ν). For each α ∈ S, let Pα = ( Pαβ ; β ∈ S) and P α = (P αβ ; β ∈ S). Both Pα and P α are probability distributions on S. If KL( Pα P α ) is the Kullback-Leibler divergence between Pα and P α , then [START_REF] Cover | Elements of Information Theory[END_REF] Q

where H( Pα ) = β∈S Pαβ log( Pαβ ) is the entropy of Pα (which does not depend on (P αβ )).

Therefore, Q(θ) reaches its maximum when each each KL( Pα P α ) reaches its minimum. But this happens when KL( Pα P α ) = 0 at P αβ = Pαβ , as required.

Proof of (23b) and (23c) : denote the second term in the expression of Q(θ) by Q(ȳ, σ) and consider the task of maximising Q(ȳ, σ).

This can be carried out by maximising each one of the functions

with respect to ȳα and σ α . From (4),

By first maximising over ȳα and recalling that η(σ α ) is negative, it is seen than any maximiser ŷα must be equal to ŷα (ω) given by (23b).

Then, it remains to maximise (over σ α ) the function

This is a strictly concave function of η(σ α ), since ψ(η) is a strictly convex function. Thus, it is possible to maximise by differentiating with respect to η(σ α ) and setting the derivative to zero. This yields the equation which is the same as equation ( 25), whose unique solution is η(σ α ) = ηα (ω) given by (23c).


Jeanne Lefevre
email: jeanne.lefevre@gipsa-lab.fr

Jonathan H Manton
email: jmanton@unimelb.edu.au

Nicolas Le Bihan
email: nicolas.le-bihan@gipsa-lab.fr

An optimisation geometry framework for the Rayleigh quotient

Keywords: global optimisation, optimisation geometry, homotopy method, Rayleigh quotient

This paper briefly introduces optimisation geometry, a method based on family of functions that proposes to solve complex optimisation problems with continuation methods and pre-computed points. As an illustration, the problem of tracking eigenvectors is presented, based on the Rayleigh quotient, and conditions under which the proposed approach is operational are detailed.

I. INTRODUCTION

The central idea of optimisation geometry is to consider optimisation problems through families of functions rather than individual cost functions [START_REF] Manton | Optimisation Geometry[END_REF]. In many applications such as inference or tracking, optimisation problems depend on an observable (some data) which may vary while the core structure of the problem remains unchanged. Classical optimisation would treat every occurrence of a new set of data as a once-off optimisation problem. We argue that the common structure of a problem of optimisation involving different observables should be used to reduce the complexity of the overall problem: for observations θ 1 , θ 2 , . . . , θ n compute local minima x * 1 , x * 2 , . . . , x * n . For instance, a classical inference problem is to fit a given distribution to some data, denoted hereafter θ 0 . This usually involves maximizing the likelihood L(x; θ 0 ) where x is a vector parametrizing the distribution. Suppose the problem is to find the best parameters for θ 0 and θ 1 two independent set of data. We assume independence between the observation as we do not treat the overlap with tracking methods, where it is assumed that θ 0 and θ 1 are close. In the previous question, we are left with two distinct optimisation problems: maximizing L(x; θ 0 ) and L(x; θ 1 ) with respect to x. In optimisation geometry, we propose to see the problem as depending on both x and θ, so that a solution found for θ 0 can be used to determine a solution for θ 1 . Suppose for instance that the structure of L is such that L(x; θ 1 ) = L(h(x); θ 0 ) where h is an application between vectors of parameters. Then it is obvious that solving the problem in θ 1 can be achieved by finding the antecedent of a maxima x * 0 of L(x; θ 0 ) by the function h. The goal is to exploit the properties of the function h, even though it has no explicit algebraic formulation. This approach allows to relocate the complexity of the problem from the graph of a function in x to the graph of a function in θ. We propose to show in this paper how optimisation geometry method is compatible with an homotopy method for the problem of estimating and tracking eigenvectors thanks to the Rayleigh quotient [START_REF] Chu | A note on the homotopy method for linear algebraic eigenvalue problems[END_REF]. The proposed approach is particularly well suited for functions defined over manifolds and could be generalized to a wide variety of optimisation problems in signal processing, i.e. source separation [START_REF] Bouchard | Riemannian optimization and approximate joint diagonalization for blind source separation[END_REF] or array processing [START_REF] Choi | Adaptive antenna array for directionof-arrival estimation utilizing the conjugate gradient method[END_REF] for example.

II. OPTIMISATION GEOMETRY

This section introduces some general results on optimisation geometry applied to a generic function f , before some illustrations depicting how this can lead to using the same numerical methods as these used by the homotopy community [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF]. Let f be a function parametrized by two quantities: x and θ. We refer to x as the parameter variable, while θ is the observation or data. The set of functions {f θ = f (•; θ)|θ ∈ Θ} where Θ is any set of possible observations, defines a family of functions of x parametrized by θ.

A. Optimisation as a tracking problem

The idea of optimisation geometry is to use a critical point x * 0 of the function f θ0 to reach a critical point x * 1 of f θ1 . A naive approach would be to use x * 0 as the initialization for an iterative method on f θ1 . This idea is generally supported by perturbation theory results, such as in [START_REF] Blum | Complexity and real computation[END_REF], which shows that under regularity assumptions on f , if θ 0 and θ 1 are in a small enough neighborhood of each other, then x * 0 and x * 1 are also within a small radius of each other. Note however that in the proposed framework, no assumptions are made on the distance between θ 0 and θ 1 . They can be very far away from each other, and so can the critical points of f θ0 and f θ1 .

Our main assumption will be the local existence of a function h so that f (x, θ 1 ) = f (h(x), θ 0 ). Then, even if h is not explicitly available, it is possible to iteratively follow its graph, initialized at a critical point (x * 0 , θ 0 ) of f θ0 up to (x * 1 , θ 1 ) a critical point of f θ1 . The next assumptions are used to prove the local existence of h and the tractability of a tracking problem on the graph of h. Let f : X × Θ → R be the cost function defined earlier. We suppose that X is a closed compact manifold and that Θ is a connected Riemannian manifold. Let θ 0 , θ 1 be two sets of data in Θ, then there is a smooth path θ: [0, 1] → Θ that satisfies θ(0) = θ 0 and θ(1) = θ 1 because Θ is smooth and connected [START_REF] Lee | Smooth manifolds[END_REF]. 955 978-9-0827-9705-3

EUSIPCO 2020

Furthermore, it is always possible to choose θ such that the derivative along the curve never vanishes. The geodesic between θ 0 and θ 1 for instance satisfies this condition. The following definition brings up the condition on a subset of a family of functions {f θ } for h to be nice. They are basically asked to be the "nicest" possible functions on a manifold, i.e. Morse functions.

Definition II.1. Let θ be a smooth path in Θ between θ 0 and θ 1 . Let

f |θ : X × θ([0, 1]) → R x , θ → f (x, θ) = f θ (x)
be the restriction of f to the image of θ. f |θ is said to be fibre wise-Morse if it is a smooth function in (x, θ) and for every θ ∈ θ([0, 1]) , f θ has only non-degenerate critical points. A critical point of f θ is called a fibre-wise critical point of f |θ .

We borrowed the term "fibre" to the fibre-bundles terminology because it is the natural extension for our framework. Similarly, we call a fibre-wise critical point of the function f a point (x * , θ) such that x * is a critical point of f θ . With the notion of fibre-wise Morse function and fibre-wise critical points in hand, we can now state the main result of this paper.

Theorem II.1. If θ is a smooth diffeomorphic path with nonvanishing derivative and f |θ is fibre-wise Morse, then the set of fibre-wise critical points Ñ of f is a smooth topologically closed submanifold of dimension one parametrizable by t.

Proof. Let f • θ: X × [0, 1] → R x , t → f θ(t) (x)
The derivative of f •θ with respect to x is:

∇ x f •θ(x, t) = ∇ f θ(t) (x).
The derivative of a function h on a manifold is the application from M to the tangent bundle T M that sends x on the vector ∇ h(x) in the tangent space at T x M [START_REF] Lee | Manifolds and differential geometry[END_REF]. Note that fibre-wise critical points denoted (θ * (t), x * ) are the points satisfying

∇ f θ * (t) (x * ) = 0 x *
where 0 x * is the 0 of the tangent plane T x * X. The set of fibre-wise critical points is then exactly the pre-image of the zero section 0

T X = {0 x ∈ T X|x ∈ X} by the map ∇ x f • θ.
According to the pre-image theorem for transverse maps [START_REF] Guillemin | Differential topology[END_REF], the set

(∇ x f •θ) -1 (0 T X ) is a closed submanifold of X×[0, 1] with no boundaries. Note that ∇ x f • θ is transverse to 0 T X because f θ is fibre-wise Morse. Furthermore, given (dx, dt) ∈ T x X × [0, 1], the kernel of the derivative D [∇ x f • θ] (x, t).(dx, dt) = Hess x f • θ(x, t)dx + ∂d ∂dt ∇ x f • θ(x, t).dt
is made of tangent vectors of the shape (0, dt), which shows that Ñ is one-dimensional and parametrizable by t.

A direct consequence of theorem II.1 is that the shape of Ñ is constrained both by f θ being fibre-wise Morse and by X t 0 1

Fig. 1. The domain X is a circle (represented on the left) and the t axis is at the bottom. Every section of the cylinder represents the space X. Every blue point a line above an fixed value of t is a potential critical point of the function f |θ(t) . In this example from left to right: components of Ñ merge, have a turning point with respect to t, cross each other,have a parametrization whose derivative goes to zero and do not admit a parametrization with respect to t (the blue circle). All this situation are prevented if f θ is fibre-wise Morse.

X t 0 1

• x 0 3 * • x 1 3 * • x 0 2 * • x 1 2 * • x 0 1 * • x 1 1 * γ 1 γ 2 γ 3 Fig. 2.
A possible shape for the set of the fibre-wise critical points Ñ of a fibre-wise Morse function f θ : the domain space X is a circle, and each t on the horizontal axis can be sent to a unique θ. Ñ is represented by the red lines, it has three connected components in this example. It can be seen that the critical points of f θ 0 and

f θ 1 are related in X × [0, 1] through Ñ
the topology of X. Figure 1 shows an example of pathological situations that can not happen under our assumptions, while figure 2 shows a possible realization for the set of fibre-wise critical points of a fibre-wise Morse function. Because Ñ is a manifold, it is made of finitely many disjoint components of dimension 1. Because it is one-dimensional, we call a component in Ñ a path. Two paths can not merge or cross. As Ñ is closed and without boundaries, a path can not stop, and because it is parametrizable by t, a tangent to a path can never be orthogonal to the axis [0, 1]. In particular, the paths can have no turning points with regard to t.

B. Homotopy methods

Homotopy methods have been known since Poincaré but their interest really came into view in the eighties, with the introduction of the fifth generation of computer whose most important feature was to allow parallel computing [START_REF] Wilson | Mathematical conversations: selections from The mathematical intelligencer[END_REF], [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF]. Indeed, homotopy methods first consisted in finding the zeros of polynomials by following curves linking the zeros of a trivial polynomial to the zeros of a more complex polynomial. By definition, the method can only be well-behaved if all curves are distinct, such that different starting points would lead to different endpoints. It is naturally easy to parallelize. On the opposite, with iterative methods such as the Newton descent it is not possible to predict which zero of the polynomial will be reached. As a consequence, in order to find all the roots of a polynomial, it is necessary to run sequentially the Newton descent, excluding at each step the roots that were already found. It should be noted that the numerical methods adapted to tracking the zeros of polynomials through a homotopy method are not restricted to polynomials and can instead be used as long as the path following problem has some nice properties such as:

1) The path γ can be continuously parametrized by a real variable t ∈ [0, 1].

2) The derivative of a parametrization γ(t) of the path never vanishes.

3) The path γ coincides with the set of attractor of an iterative method in a neighborhood containing γ. We showed in the previous paragraph that the set of fibrewise critical points Ñ of a fibre-wise Morse function f |θ satisfies all three conditions. It is made of several disconnected components, that do not merge or collapse into each other. Hence, one can find separate open neighborhoods around each one of them. Furthermore, the set of fibre-wise critical points are the zeros for the gradient ∇ x f θ , it is therefore an attractor for any iterative root-finding method. Hence, point 3) is satisfied. Furthermore, if θ is parametrized by t ∈ [0, 1] then so is Ñ . Smoothness comes from the manifold structure of Ñ , and non-vanishing derivatives of the parametrization is, as seen previously, a consequence of the regular value theorem. This shows how optimisation geometry naturally leads to the same continuation methods as the ones used in homotopy frameworks.

III. OPTIMISATION GEOMETRY APPLIED TO THE RAYLEIGH QUOTIENT

Eigenpairs of matrices are commonly approximated through iterative methods such as power and inverse iteration, QR algorithms or others [START_REF] Demmel | Applied numerical linear algebra[END_REF]. For a general matrix with no structure, these methods provide a linear rate of convergence. Algorithms for structured matrices however can reach cubic convergence [START_REF] Demmel | Applied numerical linear algebra[END_REF]. For Hermitian matrices, an efficient method is the Rayleigh quotient iteration, which possesses cubic convergence rate to one of the eigenpairs [START_REF] Helmke | Optimization and dynamical systems[END_REF]. This section recalls the definition and interest of the Rayleigh quotient and shows how optimisation applied to it leads to manifolds as the natural domain space for computations.

A. Three definitions of the Rayleigh quotient

Let A be a symmetric real matrix with eigenvalues (λ 1 , . . . , λ n ). The Rayleigh quotient with respect to A, denoted R(•; A), is the smooth function:

R: R n \{0} × M n (R) → R x ; A → x T Ax x T x
The following theorem is a special case of the Courant-Fischer minimax theorem characterizing the eigenvalues of a real symmetric n × n matrix [START_REF] Golub | Matrix computations[END_REF].

Theorem III.1. The critical points (resp. critical values) for R(•; A) are the eigenvectors (resp. eigenvalues) of A.

Note that the Rayleigh quotient is scale-invariant, that is:

∀λ ∈ R, R(x; A) = R(λx; A)
In practice, this property raises a problem for numerical methods. For instance, the Newton iteration at x k on a scaleinvariant function yields the Newton iterate N (x k ) = 2x k [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. As a consequence, the Newton method applied to the Rayleigh quotient does not converge to the eigenvectors of the matrix A, unless its initial point is actually an eigenvector. Because of its homogeneity, R is also ill-suited to an optimisation geometry framework. Let A a symmetric matrix be the data, and x the vector-parameter on which we optimize. It is obvious that for any A, R(•; A) is not a Morse function because if x * is a critical point, then the whole line defined by λx * , λ ∈ R is also made of critical points. Hence critical points of the Rayleigh quotient for any matrix A are not isolated, they are necessarily degenerate. In a sense, it is the excess of symmetry in R that generates a problem. To reduce the redundant information, one can change the domain space of the Rayleigh quotient. There are two ways to do so that are briefly introduced hereafter.

B. Rayleigh quotient on the sphere

Rather than defining R over R, where an infinity of points are sent to the same number, one can define R on the sphere S n-1 . By abuse of notation, we consider:

R s : S n-1 ×S n → R x, A → x T Ax.
The product x T Ax is not defined if x is a point on the sphere, but this expression should be understood as: φ(x) T Aφ(x) where φ is an embedding from the sphere to the set of vector with norm 1 in R n . The eigenvectors of A can be found by φ(x * ) where x * is a critical point of R. Similarly, we define -x the antipodal element to x on the sphere. Then, if φ(x * ) is an eigenvector of A, x * and -x * are critical points of R s . Furthermore, R s and R coincide on points where they are both defined: for all x in the image of the embedding φ, R s (φ -1 (x)) = R(x).

C. Rayleigh quotient on the projective space

Similarly, because R sends every line in R n to the same number, it is tempting to redifine the domain of R by quotienting R n by the relation x ∼ λx, ∀λ ∈ R. This gives the real projective space RP n-1 known as the set of directions in R n . It is a compact smooth manifold which, contrarily to the sphere, is not embedded in R n . Imagine building the projective space by gluing every point on a sphere with its antipodal. This operation in R n leads to a self-intersection, showing it can not be a faithful representation (an embedding) of RP n-1 . We define:

R p : RP n-1 × S n → R x , A → x T Ax
where x denotes an equivalence class in R n : x = {λx ; λ ∈ R}. Again, note that R p and R evaluated on the same element (modulo the equivalence relation) coincide. Furthermore, if x * and -x * are both critical points of R s A , x * is the only critical point of R p A . These restrictions of the domain of R now enable to give the conditions under which R s A or R p A are Morse functions.

Theorem III.2. Let A be a symmetric matrix in S n (R). The functions R s A and R p A are Morse if and only if all the eigenvalues of A are distinct.

Proof. For η a tangent vector to the sphere at x, the secondderivative of R s

A is [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]:

D 2 (R s A )(x).η x = 2(Id -xx T )(Aη x -η x x T Ax) if v is an eigenvector of A with eigenvalue λ: D 2 (R s A )(v).η v = 2(Id -vv T )(Aη v -λ v η v )
the last expression cancels if and only if η v is an eigenvector of A for the eigenvalue λ v . Because all eigenvalues are distinct, it cancels only for vectors proportional to v. As the tangent plane at v on S n-1 contains only vectors orthogonal to v in R n , the second derivative D 2 (R s A )(v) is strictly injective. This is true for any eigenvector of A, hence the critical points of R s A are not singular. Similarly, for η x a tangent vector of P R n-1 at x, we have

D 2 R p A (x).η x = 2P h x (Aη x -λ i η x ) where P h x = I -1 x 2 xx T is a non-vanishing matrix. Hence, D 2 R p A (x).η x = 0 ⇐⇒ η x = x which is not possible, as the tangent space T x RP n-1 = {y ∈ R n , x T y = 0} does not contain x.
The previous theorem shows that the Rayleigh quotient fits into our optimisation geometry framework as long as there exists a path Γ(t) between A and B where all the matrices Γ(t) have distinct eigenvalues. In that case only, the function R s,p Γ is fibre-wise Morse. We discuss in next section about the possibility to find such a path.

D. A path in the set of symmetric matrices

As mentioned earlier, homotopy methods are essentially path-following method, they rely on the definition of a homotopy map [START_REF] Allgower | Numerical continuation methods: an introduction[END_REF] parametrised by a variable t such that H(0, •) is a function whose zeros are known and H(1, •) is the function whose zeros are seeked. For instance, a classical map that is built to solve the problem of finding the eigenpairs of a matrix A is [START_REF] Li | Homotopy algorithm for symmetric eigenvalue problems[END_REF], [START_REF] Chu | A note on the homotopy method for linear algebraic eigenvalue problems[END_REF]:

H(t, x, λ) = ((1 -t)D + tA) x -λx, x T x -1 2
Observe that H (t, (x, λ)) = 0 ⇐⇒ (λ, x) is an eigenpair of (1 -t)D + tA with x = 1. In particular H (1, (x, λ)) = Ax, x T x-1 2 so that the set H (1, (•, •)) -1 ({0}) is exactly made of the eigenpairs of A for eigenvectors of norm 1.

Similarly, for t = 0, H (0, (x, λ)) = Dx, x T x-1 2 so that the antecedent of 0 are the eigenpairs of the matrix D. The eigenpairs of D are the starting points of the continuation method used to solve this homotopy problem, they should therefore be known in advance. This can be achieved through precomputations or more simply by choosing D diagonal which is the approach adopted in [START_REF] Li | Homotopy algorithm for symmetric eigenvalue problems[END_REF], [START_REF] Chu | A note on the homotopy method for linear algebraic eigenvalue problems[END_REF]. Note that the second term of the function is to control the norm of x, so that if H is defined on the sphere or the projective space, it's not needed anymore. It was already noticed in [START_REF] Li | Homotopy algorithm for symmetric eigenvalue problems[END_REF] that the homotopy map is well suited to continuation method if and only if all the matrices in the set {(1 -t)A + tD | t ∈ [0, 1]} have distinct eigenvalues. We find the same conditions as with the optimisation geometry framework, which is logical as the numerical method for homotopy method and optimisation geometry are interchangeable. The following result is shown in [START_REF] Chu | A note on the homotopy method for linear algebraic eigenvalue problems[END_REF]:

Proposition III.3. Let E be the set of n-tuples in C n such that D = diag(d 1 , •, d n ) satisfies: (1 -t)A + tD have distinct eigenvalues for all t ∈ [0, 1]}. Then, E is a dense open subset of C n of full Lebesgue measure in C n .
From a practical point of view, the previous theorem says that the homotopy method can be initialized with a random diagonal matrix in diag(C n ) and it gives a probability one that the homotopy map from D to A built from this matrix is well-behaved. By well-behaved, we mean that the continuation method used to numerically solve the problem is proved to converge in a finite number of steps. Interestingly, looking into the proof of this theorem shows that it is crucial to choose randomly D in C n . The set E of real n-tuples satisfying the condition of the proposition III.3 has Lebesgue measure 0 in C n . In optimisation geometry, we wish to use the precomputed points from previous iteration of an optimisation algorithm, meaning that the starting points are not picked at random. But as we saw, if the starting point is real, it should be chosen carefully to maintain the eigenvalue separation throughout the path. Given the poor reliability of the convex combination in the real case, we propose to examine a path following method based on the Riemannian geodesic between two matrices. Let A and B be two invertible positive symmetric matrices. Then, the riemannian geodesic between A and B is parametrised by: Γ

(t) = A 1 2 (A -1 2 BA -1 2 ) t A 1 2
. Note that positiveness was not a criteria until there, but because the set of invertible symmetric matrices is not connected it is convenient to choose A and B positive from here. Figure 3 compares in the 2-dimensional case A, B ∈ S + 2 (R) the Riemmannian geodesic and the convex combination. Figure 4 displays the resulting paths that would be followed by a continuation method in the cylinder S 1 × [0, 1]. It is visible on the figures that the choice of the path in S + 2 (R) greatly influences the path in the set of eigenvectors. This has its importance, as convergence of continuation methods are more or less fast given elements such as the variation of the derivatives on this path, and the absolute value of the derivative. In this precise example, a high derivative results in a curve with almost vertical sections. Hence our prediction is that for the two matrices chosen in figure 3, the straight line gives better numerical results than the Riemannian geodesic, as we can see the variations on the grey curves are smoother than the variations on the red curve. This result however can not be generalized to every possible couple of matrices.

IV. CONCLUSION

In this paper we briefly presented the idea of optimisation geometry, and showed how it compares with the work of the homotopy community on the special example on the Rayleigh quotient. We highlighted prioritary targets for future work, such as determining properties on the separate eigenvalue condition along a Riemannian geodesic in the set of symmetric positive matrices. Future work will consist in bringing more insight and theoretical results on the relation between the choice of the paths in a particular instance of optimisation geometry, and the convergence properties of optimisation methods. Rayleigh quotient is a simple example of how knowledge on the structure of the function, and in particular the data space, is necessary to adapt the path strategy. If probabilistic results are available, deterministic approaches are still an open problem, for the Rayleigh quotient as for many other functions.

Fig. 3 .

 3 Fig.3. Three dimensional illustration of paths in the set of positive semidefinite matrices: paths go from A (red dot) at t = 0 to B (blue dot) at t = 1. Black line is a straight line between the two while red line is the Riemmannian geodesic. The wireframe in the background is the limit set of positive semi-definite matrices.

Fig. 4 .

 4 Fig. 4. Paths in the set of one-dimensional eigenvectors: black and red lines are the critical point along respectively the straight line and the Riemmannian geodesic between A and B of the Rayleigh quotient. Parametrization of the path is shown both as a color gradient and on the x-axis: left-wite for t = 0, right-darker for t = 1.