Application of gradient descent algorithms based on geodesic distances

TL;DR

This paper applies Riemannian and natural gradient algorithms to optimize control and compute the Karcher mean on the manifold of positive definite Hermitian matrices using geodesic distances, demonstrating their effectiveness through simulations.

Contribution

It introduces novel gradient algorithms based on geodesic distances for control and averaging problems on positive definite matrix manifolds.

Findings

Riemannian gradient algorithm effectively controls matrix systems.

Natural gradient algorithm efficiently computes the Karcher mean.

Simulations confirm the algorithms' computational advantages.

Abstract

In this paper, the Riemannian gradient algorithm and the natural gradient algorithm are applied to solve descent direction problems on the manifold of positive definite Hermitian matrices, where the geodesic distance is considered as the cost function. The first proposed problem is control for positive definite Hermitian matrix systems whose outputs only depend on their inputs. The geodesic distance is adopted as the difference of the output matrix and the target matrix. The controller to adjust the input is obtained such that the output matrix is as close as possible to the target matrix. We show the trajectory of the control input on the manifold using the Riemannian gradient algorithm. The second application is to compute the Karcher mean of a finite set of given Toeplitz positive definite Hermitian matrices, which is defined as the minimizer of the sum of geodesic distances. To obtain more efficient iterative algorithm compared with traditional ones, a natural gradient algorithm is proposed to compute the Karcher mean. Illustrative simulations are provided to show the computational behavior of the proposed algorithms.

Equations

$$⟨A,B⟩=\mathrm{tr}\left(A^{H}B\right),$$

$$\Vert A\Vert =⟨A,A{⟩}^{\frac{1}{2}}.$$

$$g_A(X,Y):=⟨X,Y{⟩}_A=⟨A^{-1}X,A^{-1}Y{⟩}_I=\mathrm{tr}\left(A^{-1}X{A}^{-1}Y\right),$$

$$\ell (\gamma (t)):={\int }_0^1\sqrt{⟨\dot{\gamma }(t),\dot{\gamma }(t){⟩}_{\gamma (t)}}dt={\int }_0^1\sqrt{\mathrm{tr}{\left({\gamma }^{-1}(t)\dot{\gamma }(t)\right)}^2}dt.$$

$$d(A,B):=\inf \left\{\ell (\gamma )\mid \gamma :[0,1]\to P(n)\text{ with }\gamma (0)=A,\gamma (1)=B\right\}.$$

$$\gamma (t)=A^{\frac{1}{2}}{\left(A^{-\frac{1}{2}}B{A}^{-\frac{1}{2}}\right)}^t{A}^{\frac{1}{2}}\in P(n)$$

$$A\circ B=A^{\frac{1}{2}}{\left(A^{-\frac{1}{2}}B{A}^{-\frac{1}{2}}\right)}^{\frac{1}{2}}{A}^{\frac{1}{2}}$$

$$d(A,B)=\Vert \ln (A^{-\frac{1}{2}}B{A}^{-\frac{1}{2}}){\Vert }_F={\left(\sum _{i=1}^n\ln ({\lambda }_i{)}^2\right)}^{\frac{1}{2}},$$

$$u_{\ast }=\arg \underset{u}{\min }J(u),$$

$$J(u)=d^2(A(u),B).$$

$$\mathcal{M}=\left\{A(u)\mid u=(u^1,u^2,\dots .,u^m)\in \Theta \subset {\mathbb{R}}^m\right\}$$

$$A(u_{k+1})=A^{\frac{1}{2}}(u_k)\exp \left\{-{\eta }_k\ln \left(A^{-\frac{1}{2}}(u_k)B{A}^{-\frac{1}{2}}(u_k)\right)\right\}{A}^{\frac{1}{2}}(u_k)$$

$$\left({\nabla }_{A}J\right)(u_k)=A^{\frac{1}{2}}(u_k)\ln \left(A^{-\frac{1}{2}}(u_k)B{A}^{-\frac{1}{2}}(u_k)\right)A^{\frac{1}{2}}(u_k).$$

$${\exp }_A\{X\}=A^{\frac{1}{2}}\exp \left\{A^{-\frac{1}{2}}X{A}^{-\frac{1}{2}}\right\}{A}^{\frac{1}{2}},$$

$$A(u_{k+1})$$

$$=A^{\frac{1}{2}}(u_k)\exp \left\{-{\eta }_k{A}^{-\frac{1}{2}}(u_k)\left({\nabla }_{A}J\right)(u_k)A^{-\frac{1}{2}}(u_k)\right\}{A}^{\frac{1}{2}}(u_k)$$

$$=A^{\frac{1}{2}}(u_k)\exp \left\{-{\eta }_k\ln \left(A^{-\frac{1}{2}}(u_k)B{A}^{-\frac{1}{2}}(u_k)\right)\right\}{A}^{\frac{1}{2}}(u_k).$$

$${\theta }_{k+1}={\theta }_k-{\eta }_k{G}^{-1}\nabla L({\theta }_k),$$

$$\nabla L(\theta )=\left(\frac{\partial }{\partial {\theta }^1}L(\theta ),\frac{\partial }{\partial {\theta }^2}L(\theta ),\dots ,\frac{\partial }{\partial {\theta }^m}L(\theta )\right).$$

$$\frac{d}{dt}\mathrm{tr}\left(X^T(t)X(t)\right)=2\mathrm{tr}\left(X^T(t)\frac{d}{dt}X(t)\right),$$

$$\frac{d}{dt}\mathrm{tr}\left(\ln X(t)\right)=\mathrm{tr}\left(X^{-1}(t)\frac{d}{dt}X(t)\right),$$

$$\frac{d}{dt}\mathrm{tr}\left(AX(t)B\right)=\mathrm{tr}\left(A\frac{d}{dt}X(t)B\right).$$

$$u_{k+1}=u_k-{\eta }_k{G}^{-1}\nabla J(u_k),$$

$$\frac{\partial }{\partial u_k^i}J(u_k)=2\mathrm{tr}\left(B^{-\frac{1}{2}}\ln \left(B^{-\frac{1}{2}}A(u_k)B^{-\frac{1}{2}}\right)B^{\frac{1}{2}}{A}^{-1}(u_k)\frac{\partial }{\partial u_k^i}A(u_k)\right),i=1,2,\dots ,m.$$

$$u_{k+1}=u_k-{\eta }_k{G}^{-1}\nabla J(u_k),$$

$$\frac{\partial }{\partial u_k^i}J(u_k)$$

$$=2\mathrm{tr}\left(B^{-\frac{1}{2}}\ln \left(B^{-\frac{1}{2}}A(u_k)B^{-\frac{1}{2}}\right)B^{\frac{1}{2}}{A}^{-1}(u_k)\frac{\partial }{\partial u_k^i}A(u_k)\right),i=1,2,\dots ,m.$$

$$A(u)=\left(\begin{array}{cc}{u}^1& i{u}^3\\ -i{u}^3& u^2\end{array}\right),u^1>0,u^1{u}^2-(u^3{)}^2>0.$$

$$\mathcal{M}_{1}=\left\{A(y;u)=\begin{pmatrix}u^{1}&iu^{3}\\ -iu^{3}&u^{2}\end{pmatrix}~{}\Big{|}~{}u^{1}>0,u^{1}u^{2}-(u^{3})^{2}>0\right\}.$$

$$B=\left(\begin{array}{cc}55& 2i\\ -2i& 45\end{array}\right),$$

$$A(u)=\left(\begin{array}{cc}{u}^1& 0\\ 0& u^2\end{array}\right),u^1,u^2>0,$$

$$\mathcal{M}_{2}=\left\{A(u)=\begin{pmatrix}u^{1}&0\\ 0&u^{2}\end{pmatrix}~{}\Big{|}~{}u^{1},u^{2}>0\right\}.$$

Full text

Application of gradient descent algorithms based on geodesic distances

Xiaomin Duan

*School of Science, Dalian Jiaotong University, Dalian 116028, P.R. China

E-mail: dxmhope@djtu.edu.cn*

Huafei Sun

*School of Mathematics and Statistics, Beijing Institute of Technology, Beijing 100081, P.R. China

E-mail: huafeisun@bit.edu.cn*

Linyu Peng

*Waseda Institute for Advanced Study, Waseda University, Tokyo 169-8050, Japan

E-mail: l.peng@aoni.waseda.jp* Corresponding author

Abstract. In this paper, the Riemannian gradient algorithm and the natural gradient algorithm are applied to solve descent direction problems on the manifold of positive definite Hermitian matrices, where the geodesic distance is considered as the cost function. The first proposed problem is control for positive definite Hermitian matrix systems whose outputs only depend on their inputs. The geodesic distance is adopted as the difference of the output matrix and the target matrix. The controller to adjust the input is obtained such that the output matrix is as close as possible to the target matrix. We show the trajectory of the control input on the manifold using the Riemannian gradient algorithm. The second application is to compute the Karcher mean of a finite set of given Toeplitz positive definite Hermitian matrices, which is defined as the minimizer of the sum of geodesic distances. To obtain more efficient iterative algorithm compared with traditional ones, a natural gradient algorithm is proposed to compute the Karcher mean. Illustrative simulations are provided to show the computational behavior of the proposed algorithms.

Keywords: Riemannian gradient algorithm, natural gradient algorithm, system control, Karcher mean, Toeplitz positive definite Hermitian matrix

MSC: 26E60, 53B20, 93A30, 22E60

1 Introduction

Gradient adaptation is commonly applied to minimize a cost function by adjusting the parameters. Although it is often easy to implement, convergence speed of the gradient adaptation can be slow when the slope of the cost function varies widely for a small change of the parameters. To overcome the weakness of slow convergence, Amari et al. ([2, 3]) proposed the natural gradient algorithm which defines the steepest descent direction in Riemannian spaces based on the Riemannian structure of the parameter spaces. Amari also proved that the natural gradient is asymptotically Fisher-efficient for the maximum likelihood estimation, implying that it has almost the same performance as the optimal batch estimation of the parameters. The natural gradient algorithm has been widely applied into, for instance neural network, optimal control, offering a new way to solve such problems more effectively, cf. [17, 23, 25, 26].

Although the natural gradient algorithm defines the steepest descent direction, iteration trajectory of the parameters is not necessary the shortest, not to mention the difficulty to computer inverse of the metric. These problems are solved by using the Riemannian gradient algorithm in particular for matrix manifolds, separately introduced by Barbaresco [5] and Lenglet et al. [16] with wide applications, e.g. [8, 20]. It is realized that the iterative path of each parameter is along its geodesic, though the descent speed of the algorithm is not the fastest in some cases and the scope of the application is sometimes limited.

In this paper, the set of $n\times n$ positive definite Hermitian matrices is defined as a manifold $P(n)$, whose geodesic connecting two matrices was studied in Moakher [19]. Noting that the geodesic distance represents the infimum about length functions of the curves connecting two matrices, we apply both the Riemannian gradient algorithm and the natural gradient algorithm to solve the descent direction problems taking the geodesic distance as a cost function. The first problem is control of positive definite Hermitian matrix systems on manifold $P(n)$ using different gradient algorithms. Supposing the output is only determined by the control input, we take the geodesic distance as the measure of the output matrix and the target matrix. Controller to adjust the control input is shown, such that the output matrix is as close as possible to the target matrix. Trajectory of the control inout is also obtained. Second, both gradient algorithms are used to computer the Karcher mean of a finite set of given Toeplitz positive definite Hermitian matrices when sum of geodesic distances between any two matrices is viewed as the cost function. The examples show that convergence rate of the natural gradient algorithm is faster than that of the Riemannian gradient algorithm.

2 Riemannian metric and geodesics on manifold $P(n)$

Let $M(n,\mathbb{C})$ be the set of $n\times n$ complex matrices and $GL(n,\mathbb{C})$ be its subset containing only non-singular matrices. It is well known that $GL(n,\mathbb{C})$ is a Lie group, roughly speaking a group on which a differentiable manifold can also be defined. Its Lie algebra is denoted by $gl(n,\mathbb{C})$. In $M(n,\mathbb{C})$, one has the Euclidean inner product, known as the Frobenius inner product defined by

[TABLE]

where $\operatorname{tr}$ stands for the trace and the superscript $A^{\operatorname{H}}$ denotes the conjugate transport of matrix $A$. The associated norm is defined as

[TABLE]

With the above defined inner product, $M(n,\mathbb{C})$ is flat.

It is well known that the set $P(n)$ of all $n\times n$ positive definite Hermitian matrices is an $n^{2}$-dimensional manifold. Let us denote the space of all $n\times n$ Hermitian matrices by $H(n)$. The exponential map from $H(n)$ to $P(n)$ is one-to-one and onto. As $P(n)$ is an open subset of $H(n)$, for each $A\in P(n)$ we identify the set $T_{A}P(n)$ of tangent vectors to $P(n)$ at $A$. Moreover, the Riemannian metric on $P(n)$ is given by

[TABLE]

where $I$ denotes the identity element of $P(n)$ and $X,Y\in T_{A}P(n)$. The positive definiteness of this metric is a consequence of the positive definiteness of the Frobenius inner product.

Let $[0,1]$ be a closed interval in $\mathbb{R}$, and $\gamma:[0,1]\rightarrow P(n)$ be a sufficiently smooth curve on manifold $P(n)$. The length of $\gamma(t)$ is

[TABLE]

The geodesic distance between two matrices $A$ and $B$ on manifold $P(n)$ is the minimal length of curves connecting them:

[TABLE]

It transpires that length-minimizing smooth curves are geodesics, thus the infimum of (2.5) is achieved by geodesic curves. The Hopf–Rinow theorem [13] implies that $P(n)$ is geodesically complete. This means that the interval $[0,1]$ can be extended to $(-\infty,+\infty)$ and hence, for any given pair $A,B$, we can find a geodesic curve ${\gamma}(t)$ such that ${\gamma}(0)=A$ and ${\gamma}(1)=B$, namely by taking the initial velocity as $\dot{\gamma}(0)=A^{\frac{1}{2}}\ln\left(A^{-\frac{1}{2}}BA^{-\frac{1}{2}}\right)A^{\frac{1}{2}}$. Note that the length $\ell(\gamma(t))$ is invariant under congruent transformation $\gamma(t)\mapsto C\gamma(t)C^{\operatorname{H}}$, for $\forall C\in GL(n,\mathbb{C})$. As $\frac{\operatorname{d}}{\operatorname{d}\!t}\gamma^{-1}(t)=-\gamma^{-1}(t)\dot{\gamma}(t)\gamma^{-1}(t)$, one can readily see that this length is also invariant under inversion.

Let the geodesic curve $\gamma(t)$ be

[TABLE]

with $\gamma(0)=A,\gamma(1)=B$ and $\gamma^{\prime}(0)=\ln\left(A^{-\frac{1}{2}}BA^{-\frac{1}{2}}\right)\in H(n)$. Then the midpoint of $A$ and $B$, denoted as $A\circ B$, is given by

[TABLE]

and the geodesic distance $d(A,B)$ can be computed explicitly by

[TABLE]

where $\lambda_{i}$ are eigenvalues of $A^{-\frac{1}{2}}BA^{-\frac{1}{2}}$. Since $\lambda_{i}$ are also eigenvalues of $A^{-1}B$, one can compute the distance $d(A,B)$, in practice, without invoking the matrix square root $A^{-\frac{1}{2}}$.

3 Control for positive definite Hermitian matrix systems

One of the purposes in control theory is to design the control input so that the output approximates the target. Many algorithms for specific approximation problems have been proposed ([9, 15]). Among them, Zhang et al. [25] proposed a steepest descent algorithm based on the natural gradient to design the controller of an open-loop stochastic distribution control system of multi-input and single output with a stochastic noise. In biomedicine field, Hughes et al. [12] developed a control law that can anticipate meals given a probabilistic description of the patient’s eating behavior in the form of a random meal profile.

In this section, the geodesic distance on manifold $P(n)$ is taken as the cost function to solve control problems of positive definite Hermitian matrix systems. We assume that the output matrix in $P(n)$ is only determined by the control input through a system, allowing us to define the output as a matrix $A(u)\in P(n)$ as a function of the input $u=(u^{1},u^{2},\ldots,u^{m})$. Our purpose is to design the control input $u$, such that $A(u)$ is as close as possible to another target positive definite Hermitian matrix $B$ (see Fig. 3.1).

The key points for designing an algorithm are then as followings:

1. To define a distance function to measure the difference between the system output and the target.

2. To computer the trajectory of input $u$ so that the output system approximates the target, as close as possible.

In order to make the matrix $A(u)$ to be as close as possible to the given target matrix $B$, we use the geodesic distance (2.8) to measure the difference between the matrices $A(u)$ and $B$. Then we are going to design a controller and obtain the $u_{*}$ such that

[TABLE]

where the cost function $J(u)$ is defined by

[TABLE]

Let the system be well defined such that

[TABLE]

is a submanifold of manifold $P(n)$ where the control input $u$ plays the role of local coordinates.

In the following, the Riemannian gradient descent algorithm and the natural gradient descent algorithm will be used to solve this control problem, respectively. Moreover, we will analyze the suitability of two algorithms in the illustrative examples. In fact, when the target matrix $B$ lies on $\mathcal{M}$, both gradient descent algorithms are applicable. One difference is that the Riemannian gradient descent algorithm realizes the optimisation of the input trajectory while the natural gradient algorithm converges faster than the former. When the target matrix $B$ does not lie on $\mathcal{M}$, however, only the natural gradient algorithm is applicable.

3.1 The Riemannian gradient descent algorithm

Now we consider how to solve the control problem proposed above using the Riemannian gradient descent algorithm, in the case that the target matrix $B$ belongs to the output submanifold $\mathcal{M}$.

Since both the output $A(u)$ and the target matrix $B$ lie on submanifold $\mathcal{M}$, we make use of the geodesic equation to derive the trajectory and the negative gradient of the cost function $J(u)$ about $A(u)$ as the direction to give the iterative formula.

Theorem 3.1**.**

For the control input $u=(u^{1},u^{2},\ldots,u^{m})$ of a given positive definite Hermitian matrix system, the iterative formula is given by

[TABLE]

where $\eta_{k}$ is the learning rate at time $k$.

Proof.

If the gradient of $J(u)$ about $A(u)$ is denoted by $\left(\nabla_{A}J\right)(u)$, then (cf. [5])

[TABLE]

Recall that the Riemannian exponential map $\exp_{A}$ on manifold $P(n)$ is defined by

[TABLE]

where $X$ is in the tangent space $T_{A}P(n)$. Then we obtain the iterative formula as

[TABLE]

This finishes the proof. ∎

Now, we give the Riemannian gradient descent algorithm for the proposed control problem for positive definite Hermitian matrix systems.

Algorithm 3.1**.**

*For the control input $u=(u^{1},u^{2},\ldots,u^{m})$ on a given positive definite Hermitian matrix system, the iteration algorithm is

1. Set $u_{0}=(u_{0}^{1},u_{0}^{2},\ldots,u_{0}^{m})$ as an initial input. Choose a fixed learning rate $\eta$ for simplicity and a desired tolerance $\varepsilon>0$.

2. At time $k$, calculate $A(u_{k})$ using (3.4) and $d({A(u_{k})},B)$.

3. If $d({A(u_{k})},B)<\varepsilon$ then stop. Otherwise, move to step $4$.

4. Increase $k$ by one and go back to step $2$.*

Remark 3.1**.**

The initial output matrix $A(u_{0})$ will converge to the final output matrix $A(u_{*})$ along the geodesic connecting them, hence this algorithm realizes the trajectory optimisation of input $u$.

Remark 3.2**.**

If the target matrix $B$ is not on submanifold $\mathcal{M}$, it is impossible to find a geodesic on submanifold $\mathcal{M}$ such that it connects the output $A(u)$ and the target $B$. Thus, at this time, the Riemannian gradient algorithm needs be improved.

3.2 Natural gradient descent algorithm

The ordinary gradient, commonly used in learning methods on Euclidean spaces, does not give the steepest direction of a cost function on manifold, but the natural gradient does. Next, we will first introduce an important lemma about the natural gradient and then propose the natural gradient algorithm for the control problem.

Let $L(\theta)$ be a function defined in a Riemannian manifold parametrized by $\theta\in\mathbb{R}^{m}$.

Lemma 3.2** ([2]).**

The natural gradient algorithm on a Riemannian manifold is given by

[TABLE]

where $G^{-1}=(g^{ij})$ is the inverse of the Riemannian metric $G=(g_{ij})$, $L(\theta)$ is the cost function and

[TABLE]

In this subsection, we will give the natural gradient descent algorithm for the considered system from the viewpoint of information geometry. This algorithm can be applied no matter whether the target matrix $B$ is on the output submanifold $\mathcal{M}$. The following lemma is useful for computing the gradient of cost function.

Lemma 3.3** ([24]).**

Let $X(t)$ be a function-valued matrix of the real variable $t$ and let $A,B$ be constant matrices. We assume that, for all $t$ in its domain, $X(t)$ is an invertible matrix which does not have eigenvalues on the closed negative real line. Then

[TABLE]

[TABLE]

[TABLE]

Let $u=(u^{1},u^{2},\ldots,u^{m})$ be a parameter space on which a cost function $J(u)$ is defined, we get the following theorem.

Theorem 3.4**.**

The iterative process on manifold $P(n)$ is given by

[TABLE]

where the component of gradient ${\nabla}J(u_{t})$ satisfies

[TABLE]

Proof.

According to Lemma 3.2, we can get the iterative process as

[TABLE]

where the Fisher metric matrix $G$ is obtained by (2.3). Let $X(u_{k})=\ln\left(A^{-\frac{1}{2}}(u_{k})BA^{-\frac{1}{2}}(u_{k})\right)$. It is easy to show that $X(t)$ is symmetric. Using Lemma 3.3, we have the fact that

[TABLE]

This completes the proof of Theorem 3.4. ∎

From the above discussion, we formulate the natural gradient algorithm for the considered system as follows:

Algorithm 3.2**.**

*For the control of the input $u=(u^{1},u^{2},\ldots,u^{m})$ on the considered Hermitian positive definite matrix system, we have the steps that

1. Set $u_{0}=(u_{0}^{1},u_{0}^{2},\ldots,u_{0}^{m})$ as an initial input. Choose a fixed learning rate $\eta$ and a desired tolerance $\varepsilon>0$.

2. At time $k$, calculate $u_{k}$ using (3.13) and ${\nabla}J(u_{k})$.

3. If $\|{\nabla}J(u_{k})\|_{F}<\varepsilon$, stop. Otherwise, move to step $4$.

4. Increase $k$ by one and go back to step $2$.*

3.3 Simulations

From the following examples, we will show the efficiency of the two proposed algorithms, where the tolerance is $\varepsilon=10^{-15}$. In the first example, the target matrix $B$ lies in the submanifold $\mathcal{M}$, while it does not in the second example.

Example 3.1. We assume the target matrix $B$ is a point of the output submanifold, so both algorithms can be used. We choose a 3-dimensional input $u=(u^{1},u^{2},u^{3})$ and define the matrix system as

[TABLE]

The output submanifold is then

[TABLE]

We take $u_{0}^{1}=1,u_{0}^{2}=2,u_{0}^{3}=1$ as the initial input $u_{0}$ and give the target matrix $B$ by

[TABLE]

so the coordinates of the target point $B$ is $(55,45,2)$. Using Algorithm 3.1 and Algorithm 3.2, the trajectories of the input $u$ from the initial state $A(u_{0})$ to the target matrix $B$ are obtained efficiently (see Fig. 3.3).

It is easy to see that the trajectory of the input $u$ given by the Riemannian gradient algorithm is along the geodesic connecting the initial value $u_{0}$ and the target $B$ so that the path is the shortest one. In addition, although the trajectory of the input $u$ given by the natural gradient algorithm is not optimal, the convergence is faster than the Riemannian gradient algorithm (see Fig. 3.3).

Example 3.2. Now, we consider when the target matrix $B$ does not belong to the output submanifold. In this case, only the natural gradient Algorithm 3.2 is applicable to solve the control problem. Setting the input $u$ to be a 2-dimensional vector $(u^{1},u^{2})$ and the output matrix to be that

[TABLE]

so output submanifold is

[TABLE]

Let us take $u_{0}^{1}=1,u_{0}^{2}=4$ as the coordinates of the initial state $A(u_{0})$ and give the target matrix $B$ by

[TABLE]

Then, using Algorithm 3.2 to simulate the control process, we obtain the trajectory of the input $u$ from the initial state $A(u_{0})$ to the approximate matrix $A(u_{*})$ of the target matrix $B$ efficiently. The coordinate of $A(u_{*})$ is $u_{*}=(44.721,35.777)$ which can be taken as the geodesic projection of the target $B$ onto submanifold $\mathcal{M}_{2}$ (see Fig. 3.5). Furthermore, when we set the learning rate $\eta=0.1,0.2,0.5$ respectively, the efficiency and the convergence of Algorithm 3.2 are shown by Fig. 3.5.

4 Karcher mean of Toeplitz positive definite Hermitian matrices

Mean of matrices plays an important role in many fields, such as numerical analysis, probability and statistics, engineering, biological and social sciences (cf. [1, 6, 18, 22]). Many algorithms have been developed to computer such means, see e.g. ([7, 11, 21]). For electromagnetic or acoustic sensors, and more especially for radar, lidar or echography, it is necessary to consider the spatial complex data for the array processing or the time complex data for the Doppler processing: $Z_{n}=[z_{1},z_{2},\ldots,z_{n}]^{\operatorname{T}}$(cf. [5]). The covariance matrices of these complex data $R_{n}=E[Z_{n}Z_{n}^{\operatorname{H}}]$ are the Toeplitz positive definite Hermitian matrices

[TABLE]

where $r_{k}=E[z_{n}\overline{z}_{n-k}]$ and $Z^{\operatorname{H}}R_{n}Z>0,$ for $\forall Z\in\mathbb{C}^{n}$.

In this section, we consider the complex circular multivariate Gaussian distribution of zero mean with probability density function

[TABLE]

Let us denote the set of all $n\times n$ Toeplitz positive definite Hermitian matrices by $Sym(n,\mathbb{C})$ which is obviously an $(n^{2}-n+1)$-dimensional submanifold of manifold $P(n)$. In the Riemannian sense, the mean $\overline{R}$ of $N$ given positive definite Hermitian matrices $R^{1},R^{2},\ldots,R^{N}$ is defined as [10]

[TABLE]

which is called the Karcher mean. For $N$ distributions $p(\cdot|R^{i}),i=1,2,\ldots,N$, we denote their covariance matrices by $R^{k}$ and define the cost function by

[TABLE]

Note that the local curvature of complex circular multivariate Gaussian distribution of zero mean is a non-positive constant, so the Karcher mean is unique [14].

In [4], it was shown that the Jacobi field for the Karcher mean vanishes. In order to estimate the Doppler ambiance, Barbaresco [5] computed the Jacobi field and proposed the Riemannian gradient algorithm to compute the Karcher Mean by

[TABLE]

with $\eta$ the learning rate.

In the following, we will propose the natural gradient algorithm to compute the Karcher mean of $N$ Toeplitz positive definite Hermitian matrices $R^{i}(i=1,2,\ldots,N)$ followed with simulations.

4.1 Natural gradient descent algorithm

Let $\theta=(\theta^{1},\theta^{2},\ldots,\theta^{m})$ be a parameter space on which a function $L(\theta)$ is defined. Analogously to the proof of Theorem 3.4, we have the following theorem:

Theorem 4.1**.**

The iterative process on manifold $Sym(n,\mathbb{C})$ is given by

[TABLE]

where

[TABLE]

and the component of gradient ${\nabla}L(\theta)$ satisfies

[TABLE]

Now we are ready to formulate the natural gradient algorithm.

Algorithm 4.1**.**

*The natural gradient algorithm to computer the Karcher mean of $N$ matrices of the manifold $Sym(n,\mathbb{C})$ is as follows:

1. Take the arithmetic mean $\frac{1}{N}\sum_{i=1}^{N}R^{i}$ as the initial point $\theta_{0}$. Choose a learning rate $\eta$ and a desired tolerance $\varepsilon>0$.

2. At time $k$, calculate $\theta_{k}$ using (4.6) and ${\nabla}L(\theta_{k})$.

3. If $\|{\nabla}L(\theta_{k})\|<\varepsilon$, stop. Otherwise, move to step $4$.

4. Increase $k$ by one and go back to step $2$.*

4.2 Simulations

In this subsection, we use the two gradient descent algorithms mentioned above to compute Karcher mean of $N$ given Toeplitz positive definite Hermitian matrices, where the tolerance is again $\varepsilon=10^{-15}$. From the following examples, it is shown that the natural gradient algorithm is more efficient than the algorithm (4.5).

For simplicity, we choose a 2-dimensional spatial complex data for the array processing or the time complex data for the Doppler processing: $Z_{2}=[z_{1},z_{2}]^{\operatorname{T}}$. Then, the covariance matrices of these complex data can be written as

[TABLE]

Example 4.1. We computer the Karcher mean of $R^{1},R^{2}$ on manifold $Sym(2,\mathbb{C})$ as a first example, where

[TABLE]

In fact, it is easy to know their Karcher mean $R^{1}\circ R^{2}$ corresponds to (2.7). To show the algorithm’s efficiency, the natural gradient algorithm and the algorithm (4.5) are used respectively. Choose the arithmetic mean $\frac{1}{2}\sum_{i=1}^{2}R^{i}$ as the initial point $\theta_{0}$. On manifold $Sym(2,\mathbb{C})$, the adjustment of the coordinate vector $\theta$ is given by (4.5) and (4.6). From (2.6), the geodesic between $R^{1}$ and $R^{2}$ is computed (see Fig. 4.1). We denote the Karcher mean by a red pentacle. It is shown that both algorithms converge to the red pentacle. It is also shown that the natural gradient algorithm (4.6) converges faster than the algorithm (4.5) (see Fig. 4.2).

Example 4.2. Here, we consider three $2\times 2$ Toeplitz positive definite Hermitian matrices

[TABLE]

Using (2.6), we can get the geodesics between each two of the three points $R^{1},R^{2}$ and $R^{3}$ on $Sym(2,\mathbb{C})$, which form a geodesic triangle (see Fig. 4.4). The midpoint of each geodesic is obtained using (2.7). Thus, each median connects a vertex with the midpoint of its opposing side. In Euclidean spaces, these centerlines always meet in a single point which is the Karcher mean. However, in curved spaces, as Fig. 4.4 shows, this is no longer true. In fact, when we change the look angle to Fig. 4.4, it is shown that three midlines are very likely non-coplanar (see Fig. 4.4).

Now we compute the Karcher mean $p(\cdot|\overline{R})$ of 3 distributions $p(\cdot|R^{i}),i=1,2,3$ using these two algorithms. Similarly to Example 4.1, first we still take the arithmetic mean $\frac{1}{3}\sum_{i=1}^{3}R^{i}$ as the initial point $\theta_{0}$ . On manifold $Sym(2,\mathbb{C})$, the adjustment of the coordinate vector $\theta$ is given by (4.5) and (4.6). Finally, both algorithms converge to the Karcher mean (see Fig. 4.5)

[TABLE]

Again the natural gradient algorithm is faster than the algorithm (4.5) (see Fig. 4.6).

5 Conclusions

In this paper, we applied the Riemannian gradient algorithm and the natural gradient algorithm to the control of positive Hermitian matrix systems as well as the computation of Karcher mean of Toeplitz positive definite Hermitian matrices. Their behaviors were also compared.

For the control system, when the target matrix belongs to the output submanifold, the Riemannian gradient descent algorithm and the natural gradient algorithm are both applicable. It was shown that the former realizes the optimisation of input trajectory while the latter has preferable convergence. We also proposed the natural gradient algorithm to compute the Karcher mean of matrices in the submanifold $Sym(n,\mathbb{C})$, taking the sum of geodesic distances as the cost function. The simulations showed that convergence rate of the natural gradient algorithm is faster than the Riemannian gradient algorithm, which has been widely used to solve such optimisation problems during the last decade.

Acknowledgements

X.D. is supported by the National Natural Science Foundation of China (No. 61401058) and the Natural Science Foundation of Liaoning Province (No. 20180550112). H.S. is partially supported by the National Natural Science Foundation of China (No. 61179031, No. 10932002). L.P. is supported by JSPS Grant-in-Aid for Scientific Research (No. 16KT0024), the MEXT “Top Global University Project”, Waseda University Grant for Special Research Projects (No. 2019C-179, No. 2019E-036) and Waseda University Grant Program for Promotion of International Joint Research.