A Geodesic-Based Riemannian Gradient Approach to Averaging on the Lorentz Group

Abstract

In this paper, we propose an efficient algorithm to solve the averaging problem on the Lorentz group $O(n,k)$. Firstly, we introduce the geometric structures of $O(n,k)$ endowed with a Riemannian metric where geodesic could be written in closed form. Then, the algorithm is presented based on the Riemannian-steepest-descent approach. Finally, we compare the above algorithm with the Euclidean gradient algorithm and the extended Hamiltonian algorithm. Numerical experiments show that the geodesic-based Riemannian-steepest-descent algorithm performs the best in terms of the convergence rate.

1. Introduction

Averaging measured data on matrix manifolds is one of the most frequently arising problems in many fields. For example, one may need to determine an average eye of a set of optical systems [1]. In [2], the authors concern the statistics of covariance matrices for biomedical engineering applications. The current research mainly focuses on the compact Lie groups, such as the special orthogonal group $SO\left(n\right)$ [3] and the unitary group $U\left(n\right)$ [4], or other matrix manifolds, such as the special Euclidean group $SE(n,\mathbb{R})$ [5], Grassmann manifold $Gr(n,k)$ [6,7], the Stiefel manifold $St(n,k)$ [8,9] and the symmetry positive-definite matrix manifold $SPD\left(n\right)$ [10,11]. However, there is little research on averaging over the Lorentz group $O(n,k)$.

In [12], Clifford algebra, as a generalization of quaternions, is applied to approximate the average on the special Lorentz group $SO(1,2)$. In engineering, the Lorentz group plays a fundamental role, e.g., in the analysis of motion of charged particles in electromagnetism [13]. In physics, the Lorentz group forms a basis for the transformation of parabolic catadioptric image features [14]. According to [15], the Lorentz group $O(3,1)$ is equal to the linear transformations of the projective plane and used to characterize the catadioptric fundamental matrices that are applied to study the catadioptric cameras.

This paper aims at investigating the averaging on the Lorentz group with a left-invariant metric that plays the role of Riemannian metric so that the geodesic is given in closed form. The considered averaging optimization problem could be solved numerically via a geodesic-based Riemannian-steepest-descent algorithm (RSDA). Furthermore, the devised method RSDA is compared with the line-search algorithm, Euclidean gradient algorithm (EGA) and the second-order learning algorithm, extended Hamiltonian algorithm (EHA), proposed in [16,17].

The remainder of the paper is organized as follows. In Section 2, we summarize the geometric structures of the Lorentz group. The averaging optimization problem, the geodesic-based Riemannian-steepest-descent algorithm and the extended Hamiltonian algorithm on the Lorentz group are presented in Section 3. In Section 4, we show the numerical results of the RSDA and compare the convergence rate of the RSDA with the other two algorithms. Section 5 presents the conclusions of the paper.

2. Geometry of the Lorentz Group

In this section, we review the foundations of Lorentz group and describe a Riemannian metrization in which the geodesic could be written in closed forms.

Let $M(n,\mathbb{R})$ be the set of $n\times n$ real matrices and $GL(n,\mathbb{R})$ be its subset containing only nonsingular matrices. It is well-known that $GL(n,\mathbb{R})$ is a Lie group, i.e., a group which is also a differentiable manifold and for which the operations of group multiplication and inverse are smooth.

Definition 1.

The Lorentz group $O(n,k)$ as a Lie subgroup of $GL(n+k,\mathbb{R})$ is defined by

$$\begin{array}{c}\hfill O(n,k):=\{A\in GL(n+k,\mathbb{R})\mid A^{\top }{I}_{n,k}A=I_{n,k}\},I_{n,k}:=\left(\begin{array}{cc}{I}_n& {\mathbf{0}}_{n\times k}\\ {\mathbf{0}}_{k\times n}& -I_k\end{array}\right),\end{array}$$

where I is identity matrix, and $\mathbf{0}$ is $n\times k$ zero matrix.

For any $A\in O(n,k)$, the following identities hold:

$$det\left(A\right)=\pm 1,I_{n,k}^2=I_{n+k},I_{n,k}^{\top }=I_{n,k}^{-1}=I_{n,k},A^{\top }=I_{n,k}{A}^{-1}{I}_{n,k},A^{-\top }=I_{n,k}A{I}_{n,k}.$$

The tangent space $T_{A}O(n,k)$, the Lie algebra $\mathfrak{o}(n,k)$ and the normal space $N_{A}O(n,k)$ associated with the Lorentz group $O(n,k)$ can be characterized as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}\hfill T_{A}O(n,k)& =\{A{I}_{n,k}\mathsf{\Omega }\mid \mathsf{\Omega }\in {\mathbb{R}}^{(n+k)\times (n+k)},{\mathsf{\Omega }}^{\top }=-\mathsf{\Omega }\},\hfill \\ \hfill \mathfrak{o}(n,k)& =\{I_{n,k}\mathsf{\Omega }\mid \mathsf{\Omega }\in {\mathbb{R}}^{(n+k)\times (n+k)},{\mathsf{\Omega }}^{\top }=-\mathsf{\Omega }\},\hfill \\ \hfill N_{A}O(n,k)& =\{I_{n,k}AS\mid S\in {\mathbb{R}}^{(n+k)\times (n+k)},S^{\top }=S\}.\hfill \end{array}\end{array}\right.\end{array}$$

We employ the left-invariant metric as the Riemannian metric on $T_{A}O(n,k)$ at $A\in O(n,k)$

$${\langle X,Y\rangle }_A:={\langle A^{-1}X,A^{-1}Y\rangle }_{I_{n+k}}=tr({\left(A^{-1}X\right)}^{\top }\left(A^{-1}Y\right)),X,Y\in T_{A}O(n,k).$$

(1)

Let $\gamma :[0,1]\to O(n,k)$ be a smooth curve in $O(n,k)$, and its length associated with the Riemannian metric (1) is defined as

$$L\left(\gamma \right)={\int }_0^1\sqrt{{\langle \dot{\gamma }\left(t\right),\dot{\gamma }\left(t\right)\rangle }_{\gamma \left(t\right)}}dt.$$

(2)

The geodesic distance between two points $A_1$ and $A_2$ in $O(n,k)$ is the infimum of lengths of curves connecting them, namely,

$$d(A_1,A_2):=inf\left\{L\left(\gamma \right)\right|\gamma :[0,1]\to O(n,k),\gamma \left(0\right)=A_1,\gamma \left(1\right)=A_2\}.$$

Applying the variational method from Equation (2), we can find that the geodesic equation satisfies

$$\ddot{\gamma }+{\Gamma }_{\gamma }(\dot{\gamma },\dot{\gamma })=0,$$

where $\Gamma$ denotes the Christoffel symbol, the over-dot and the double over-dot denote first-order and second-order derivation with respect to parameter t, respectively. When the values $\gamma \left(0\right)=A\in O(n,k)$ and $\dot{\gamma }\left(0\right)=X\in T_{A}O(n,k)$ are specified, we shall denote the solution as ${\gamma }_{A,X}\left(t\right)$, which is associated with the Riemannian exponential map ${exp}_A:T_{A}O(n,k)\to O(n,k)$ defined as

$${exp}_A\left(X\right)={\gamma }_{A,X}\left(1\right).$$

Let $f:O(n,k)\to \mathbb{R}$ be a differentiable function and ${\nabla }_{A}f\in T_{A}O(n,k)$ denotes the Riemannian gradient of function f with respect to the metric (1), which is defined by the compatibility condition

$${\langle {\nabla }_{A}f,X\rangle }_A={\langle {\partial }_{A}f,X\rangle }^E,\forall X\in T_{A}O(n,k),$$

(3)

where ${\partial }_{A}f$ denotes the Euclidean gradient of function f and ${\langle ·,·\rangle }^E$ denotes the Euclidean metric.

The exponential of a matrix A in $M(n,\mathbb{R})$ is defined by

$$e^A=\sum _{k=0}^{\infty }\frac{1}{k!}{A}^k.$$

We remark that $e^{A+B}=e^A{e}^B$ when $AB=BA$ for A and B in $M(n,\mathbb{R})$.

For the Frobenius norm $\parallel ·\parallel$, when $\parallel A-I_n\parallel <1$, the logarithm of A is given by

$$logA=\sum _{k=1}^{\infty }{(-1)}^{k+1}\frac{{(A-I_n)}^k}{k}.$$

In general, $log\left(AB\right)\ne logA+logB$. We here recall the important fact that

$$log\left(B^{-1}AB\right)=B^{-1}log\left(A\right)B,B\in GL(n,\mathbb{R}).$$

Lemma 1.

Suppose that $f\left(A\right)$ is a scalar function of $n\times n$ matrix A in $M(n,\mathbb{R})$ [18]. If

$$df\left(A\right)=tr(WdA)$$

holds, then

$${\partial }_{A}f=W^{\top },$$

where $W^{\top }$ denotes the derivative of $f\left(A\right)$.

The structure of the Riemannian gradient associated with the Riemannian metric (1) is given by the following result.

Theorem 1.

The Riemannian gradient of a sufficiently regular function $f:O(n,k)\to \mathbb{R}$ associated with the Riemannian metric (1) satisfies

$${\nabla }_{A}f=\frac{1}{2}A{I}_{n,k}(I_{n,k}{A}^{\top }{\partial }_{A}f-{\partial }_A^{\top }fA{I}_{n,k}).$$

(4)

Proof of Theorem 1.

According to equality (3), the gradient ${\nabla }_{A}f$ is computed as

$${\langle {\partial }_{A}f-A^{-\top }{A}^{-1}{\nabla }_{A}f,X\rangle }^E=0,\forall X\in T_{A}O(n,k),$$

as X is arbitrary, the condition above implies that ${\partial }_{A}f-A^{-\top }{A}^{-1}{\nabla }_{A}f\in N_{A}O(n,k)$, hence ${\partial }_{A}f-A^{-\top }{A}^{-1}{\nabla }_{A}f=I_{n,k}AS$, with $S^{\top }=S$, and thus we have

$${\nabla }_{A}f=A{A}^{\top }{\partial }_{A}f-A{I}_{n,k}S.$$

Since ${\nabla }_{A}f\in T_{A}O(n,k)$, we have

$${\left({\nabla }_{A}f\right)}^{\top }{I}_{n,k}A+A^{\top }{I}_{n,k}{\nabla }_{A}f=0.$$

Inserting the expression ${\nabla }_{A}f$ into the equation above, we have

$$S=\frac{1}{2}(I_{n,k}{A}^{\top }{\partial }_{A}f+{\partial }_A^{\top }fA{I}_{n,k}).$$

Substituting S into the expression of ${\nabla }_{A}f$, we finish the proof. ☐

In the text, we would like to give credit to [19] for the explicit expression of geodesic over the general linear group endowed with the natural left-invariant Riemannian metric (1). Under the same Riemannian metric, it is clear that the geodesic on the Lorentz group can be written in closed form (cf. Theorem 2.14 in [19]).

Proposition 1.

The geodesic ${\gamma }_{A,X}:[0,1]\to O(n,k)$, with $A\in O(n,k)$, $X\in T_{A}O(n,k)$ corresponding to the Riemannian metric (1) has expression

$${\gamma }_{A,X}\left(t\right)=A{e}^{t{(A^{-1}X)}^{\top }}{e}^{t[A^{-1}X-{(A^{-1}X)}^{\top }]}.$$

(5)

Proof of Proposition 1.

In fact, by the variational formula $\delta {\int }_0^1{\langle \dot{\gamma },\dot{\gamma }\rangle }_{\gamma }dt=0$, the geodesic equation in Christoffel form reads

$$\ddot{\gamma }+\dot{\gamma }{\dot{\gamma }}^{\top }{\gamma }^{-\top }-\gamma {\dot{\gamma }}^{\top }{\gamma }^{-\top }{\gamma }^{-1}\dot{\gamma }-\dot{\gamma }{\gamma }^{-1}\dot{\gamma }=0.$$

By the definition that $H\left(t\right):={\gamma }^{-1}\left(t\right)\dot{\gamma }\left(t\right)$, we see that the geodesic satisfies the Euler–Lagrange equation

$$\dot{H}\left(t\right)=H^{\top }\left(t\right)H\left(t\right)-H\left(t\right)H^{\top }\left(t\right).$$

We can verify that the curve $\gamma \left(t\right)=A{e}^{t{(A^{-1}X)}^{\top }}{e}^{t[A^{-1}X-{(A^{-1}X)}^{\top }]}$ emanating from A with a velocity X is the solution of the Euler–Lagrange equation above. ☐

Remark 1.

From the proof of Proposition 1 and the fact $A^{-\top }=I_{n,k}A{I}_{n,k}$, $\forall X\in T_{A}O(n,k)$, we conclude that

$$\begin{array}{cc}\hfill {\Gamma }_A(X,X)& =X{X}^{\top }{A}^{-\top }-A{X}^{\top }{A}^{-\top }{A}^{-1}X-X{A}^{-1}X\hfill \\ & =X{X}^{\top }{I}_{n,k}A{I}_{n,k}-A{X}^{\top }{I}_{n,k}A{I}_{n,k}{A}^{-1}X-X{A}^{-1}X.\hfill \end{array}$$

3. Optimization on the Lorentz Group

This section aims at investigating the averaging optimization problem on the Lorentz group, and finding the average value out of a set of Lorentz matrices. Section 3.1 and Section 3.2 recall the geodesic-based Riemannian-steepest-descent algorithm and the extended Hamiltonian algorithm on the Lorentz group, respectively.

3.1. Riemannian-Steepest-Descent Algorithm on the Lorentz Group

According to the expression of geodesic (5) with $A=\gamma \left(0\right)$, $B=\gamma \left(1\right)\in O(n,k)$ and $X=\dot{\gamma }\left(0\right)\in T_{A}O(n,k)$, we have

$$A^{-1}B=e^{{(A^{-1}X)}^{\top }}{e}^{A^{-1}X-{(A^{-1}X)}^{\top }},$$

(6)

and, furthermore, we have

$$log(A^{-1}B)=log[e^{{(A^{-1}X)}^{\top }}{e}^{A^{-1}X-{(A^{-1}X)}^{\top }}].$$

(7)

Note that the geodesic distance $d(A,B)=\sqrt{{\langle X,X\rangle }_A}=\parallel A^{-1}X\parallel$. From Equation (6), it is difficult to give an explicit expression of X in A and B. In [20], the authors show that the geodesic distance connecting A and B could be obtained by solving the following problem

$$\underset{U\in \mathfrak{o}(n,k)}{min}{\parallel e^{U^{\top }}{e}^{U-U^{\top }}-A^{-1}B\parallel }^2.$$

(8)

Then, the geodesic distance between A and B, $d(A,B)=\parallel U^{*}\parallel$, where $U^{*}$ is a minimizer of Equation (8). In order to compute $U^{*}$, a descent procedure can be applied so we need to find the gradient first. The differential with respect to U of the function $J\left(U\right)=\parallel e^{U^{\top }}{e}^{U-U^{\top }}-A^{-1}{B\parallel }^2$ is given by

$$\begin{array}{cc}\hfill dJ\left(U\right)& =2d(tr({(e^{U^{\top }}{e}^{U-U^{\top }}-A^{-1}B)}^{\top }(e^{U^{\top }}{e}^{U-U^{\top }}-A^{-1}B)))\hfill \\ & =2tr({(e^{U^{\top }}{e}^{U-U^{\top }}-A^{-1}B)}^{\top }{e}^{U^{\top }}dU{e}^{U-U^{\top }})\hfill \\ & =2tr(e^{U-U^{\top }}{(e^{U^{\top }}{e}^{U-U^{\top }}-A^{-1}B)}^{\top }{e}^{U^{\top }}dU),\hfill \end{array}$$

and thus, according to Lemma 1, we can get the Euclidean gradient of $J\left(U\right)$

$${\partial }_{U}J=2e^U(e^{U^{\top }}{e}^{U-U^{\top }}-A^{-1}B)e^{U^{\top }-U}.$$

(9)

On the other hand, recall the Goldberg exponential formula in [21]

$$log\left(e^X{e}^Y\right)=X+Y+\sum _{n=2}^{\infty }\sum _{\left|\omega \right|=n}{g}_{\omega }\omega ,$$

(10)

where $\omega$ denotes a word in the letters X and Y, and the inner sum is over all words $\omega$ with length $\left|\omega \right|=n$. The symbol $g_{\omega }$ denotes Goldberg coefficient on the word $\omega$, a rational number. If the word $\omega$ begins with X, say

$$\omega ={\omega }_{X,Y}=X^{r_1}{Y}^{r_2}{X}^{r_3}\cdots {(X\bigvee Y)}^{r_m},$$

where $r_1,\cdots ,r_m$ are positive integers satisfying $r_1+\cdots +r_m=n$, and ${(X\bigvee Y)}^{r_m}$ denotes $X^{r_m}$ if m is odd and $Y^{r_m}$ if m is even. ${\omega }_{Y,X}$ is similarly defined.

The Goldberg coefficient $g_{\omega }$ is

$$g_{\omega }={\int }_0^1{t}^{m^{\prime }}{(t-1)}^{m^{\prime \prime }}{G}_{r_1}\left(t\right)\cdots G_{r_m}\left(t\right)dt,$$

in which $m^{\prime }=[m/2]$ and $m^{″}=[(m-1)/2]$, $\left[\right]$ denotes the greatest integer of the enclosed number, and the polynomials $G_r\left(t\right)$ are defined recursively by $G_1\left(t\right)=1$ and $G_r\left(t\right)=r^{-1}(d/dt)\left\{r(r-1)G_{r-1}\left(t\right)\right\}$ for $r=2,3,\dots$.

According to Theorem 1 in [22], we can obtain the following result.

Proposition 2.

If $\parallel U\parallel <\frac{1}{2}$, then

$$\parallel log\left(A^{-1}B\right)-U\parallel \le \frac{{4\parallel U\parallel }^2}{1-2\parallel U\parallel }.$$

(11)

Proof of Proposition 2.

According to the Goldberg exponential formula (10), Equation (7) can be recast as

$$\begin{array}{c}\hfill log\left(A^{-1}B\right)=log(e^{U^{\top }}{e}^{U-U^{\top }})=U+\sum _{n=2}^{\infty }\sum _{\left|\omega \right|=n}{g}_{\omega }\omega ,\end{array}$$

(12)

in which $\omega$ denotes a word in the letters $U^{\top }$ and $U-U^{\top }$.

By using properties of the Goldberg coefficient $g_{\omega }$ (see [23]), we get

$$|g_{\omega }|\le 2^{-(n-m)}(m{C}_{m-1}^{m^{\prime }}),$$

where $C_n^m=n!/m!(n-m)!$ is the usual binomial coefficient. Furthermore, the number of words ${\omega }_{U^{\top },U-U^{\top }}$ of length $\left|\omega \right|=n$ is just $C_{n-1}^{m-1}$. Thus, the sum of $|g_{\omega }|$ extended over all words $\omega$ of length n and involved m parts is

$$2\times C_{n-1}^{m-1}\times 2^{-(n-m)}{(m{C}_{m-1}^{m^{\prime }})}^{-1}.$$

Therefore, due to $\parallel U-U^{\top }\parallel \le 2\parallel U\parallel$, and for $n\ge 2$, $\parallel \omega \parallel \le 2^{n-1}{\parallel U\parallel }^n$, we have

$$\begin{array}{cc}\hfill \parallel log\left(A^{-1}B\right)-U\parallel & \le \sum _{n=2}^{\infty }\sum _{\left|\omega \right|=n}|g_{\omega }|\parallel \omega \parallel \hfill \\ & \le \sum _{n=2}^{\infty }\sum _{m=1}^{n}2\times C_{n-1}^{m-1}\times 2^{-(n-m)}{(m{C}_{m-1}^{m^{\prime }})}^{-1}\times 2^{n-1}{\parallel U\parallel }^n\hfill \\ & \le \sum _{n=2}^{\infty }\sum _{m=1}^n{2}^m{C}_{n-1}^{m-1}{\left(2^{m-1}\right)}^{-1}{\parallel U\parallel }^n\hfill \\ & =2\sum _{n=2}^{\infty }\sum _{m=1}^n{C}_{n-1}^{m-1}{\parallel U\parallel }^n\hfill \\ & =\sum _{n=2}^{\infty }{2}^n{\parallel U\parallel }^n\hfill \\ & =\frac{{4\parallel U\parallel }^2}{1-2\parallel U\parallel }.\hfill \end{array}$$

Here, in the last step, the equality holds if and only if $\parallel U\parallel <\frac{1}{2}$. ☐

Remark 2.

If the condition $\parallel U-U^{\top }\parallel \le ϵ\parallel U\parallel$, $0<ϵ\le 1$ is added, then, when $\parallel U\parallel <1$,

$$\parallel log\left(A^{-1}B\right)-U\parallel \le \frac{{2ϵ\parallel U\parallel }^2}{1-\parallel U\parallel }.$$

Remark 3.

Proposition 2 implies that we can replace the geodesic distance $d(A,B)=\parallel U\parallel$ by $\parallel log\left(A^{-1}B\right)\parallel$ where the latter one is much easier to handle.

In order to compare different Lorentz matrices, the following discrepancy $D:O(n,k)\times O(n,k)\to \mathbb{R}$ is defined

$$D(A,B):=\parallel log\left(A^{-1}B\right)\parallel =\sqrt{tr\left({log}^{\top }\left(A^{-1}B\right)log\left(A^{-1}B\right)\right)},A,B\in O(n,k).$$

(13)

For the Lorentz group, the criterion function $f:O(n,k)\to \mathbb{R}$ to measure a collection $\{B_1,B_2,\dots ,B_N\}$ of samples randomly generated is defined by

$$f\left(A\right):=\frac{1}{N}\sum _q{\parallel log\left(B_q^{-1}A\right)\parallel }^2,A,B_q\in O(n,k),q=1,2,\cdots ,N,$$

(14)

which appears as a specialization, to the Lorentz group, of the Karcher criterion [24].

Definition 2.

The average value of $\{B_1,B_2,\dots ,B_N\}$ is defined as

$$\overline{A}:=arg\underset{A\in O(n,k)}{min}f\left(A\right).$$

(15)

Based on the case in the Euclidean space of the line search method with the negative Euclidean gradient, in order to achieve the average value, a geodesic-based Riemannian-steepest-descent algorithm can be expressed as [25]

$$A_{p+1}={\gamma }_{A_p,-{\nabla }_{A_p}f}\left(t_p\right),$$

(16)

where $p\ge 0$ is an iteration step-counter, the initial guess $A_0\in O(n,k)$ is arbitrarily chosen, and $t_p$ denotes a (possibly variable) stepsize schedule.

Note that the Riemannian gradient of the criterion function is crucial to solving the optimization problem by the geodesic-based Riemannian-steepest-descent algorithm (16).

Theorem 2.

The Riemannian gradient of the criterion function (14) is

$$\begin{array}{c}\hfill {\nabla }_{A}f=\frac{2}{N}\sum _{q}Alog\left(B_q^{-1}A\right).\end{array}$$

(17)

Proof of Theorem 2.

The differential with respect to A of the criterion function (14) is

$$\begin{array}{cc}\hfill df\left(A\right)& =\frac{1}{N}\sum _{q}d(tr({log}^{\top }\left(B_q^{-1}A\right)log\left(B_q^{-1}A\right))\hfill \\ & =\frac{2}{N}\sum _{q}tr({log}^{\top }\left(B_q^{-1}A\right)dlog\left(B_q^{-1}A\right))\hfill \\ & =\frac{2}{N}\sum _{q}tr({log}^{\top }\left(B_q^{-1}A\right)A^{-1}dA)\hfill \\ & ={〈{\left(\frac{2}{N}\sum _q({log}^{\top }\left(B_q^{-1}A\right)A^{-1}\right)}^{\top },dA〉}^E,\hfill \end{array}$$

and, furthermore, according to Lemma 1, we get the Euclidean gradient

$${\partial }_{A}f=\frac{2}{N}\sum _q{A}^{-\top }log\left(B_q^{-1}A\right).$$

(18)

Substituting the Euclidean gradient (18) into the expression of the Riemannian gradient (4), we have

$$\begin{array}{cc}\hfill {\nabla }_{A}f& =\frac{1}{N}\sum _q(Alog\left(B_q^{-1}A\right)-A{I}_{n,k}{log}^{\top }\left(B_q^{-1}A\right)I_{n,k})\hfill \\ & =\frac{1}{N}\sum _q(Alog\left(B_q^{-1}A\right)-Alog(I_{n,k}{\left(B_q^{-1}A\right)}^{\top }{I}_{n,k}))\hfill \\ & =\frac{1}{N}\sum _q(Alog\left(B_q^{-1}A\right)-Alog{\left(B_q^{-1}A\right)}^{-1})\hfill \\ & =\frac{2}{N}\sum _{q}Alog\left(B_q^{-1}A\right).\hfill \end{array}$$

 ☐

Associated with the expressions of the geodesic (5) and the Riemannian gradient (17) of the criterion function, the geodesic-based Riemannian-steepest-descent algorithm (16) is recast as

$$A_{p+1}=A_p{e}^{t_p{U}_p^{\top }}{e}^{t_p(U_p-U_p^{\top })},$$

(19)

where $U_p:=-A_p^{-1}{\nabla }_{A_p}f=-\frac{2}{N}{\sum }_{q}log\left(B_q^{-1}{A}_p\right)$.

In order to compute a nearly-optimal stepsize, the function $f\left(A_{p+1}\right)$ considered as a function of the stepsize $t_p$ can be expanded in Taylor series at the point $t_p=0$ as

$$f\left(A_{p+1}\right)=\frac{1}{2}{f}_{2,p}{t}_p^2+f_{1,p}{t}_p+f\left(A_p\right)+o\left(t_p^2\right),$$

where the $f_{1,p}$ and $f_{2,p}$ as the coefficients of the Taylor expansion are computed by

$$f_{1,p}=\frac{df\left(A_{p+1}\right)}{d{t}_p}{|}_{t_p=0},f_{2,p}=\frac{d^{2}f\left(A_{p+1}\right)}{d{t}_p^2}{|}_{t_p=0}.$$

Furthermore, the nearly-optimal stepsize value that minimizes the $f\left(A_{p+1}\right)$ is ${\widehat{t}}_p:=-\frac{f_{1,p}}{f_{2,p}},$ provided that $f_{1,p}\le 0$, $f_{2,p}>0$, hence ${\widehat{t}}_p\ge 0$.

Thus, to obtain the $f_{1,p}$ and $f_{2,p}$ , the function $f\left(A_{p+1}\right)$ is recast as

$$f\left(A_{p+1}\right)=\frac{1}{N}\sum _{q}tr({log}^{\top }(B_q^{-1}{A}_p{e}^{t_p{U}_p^{\top }}{e}^{t_p(U_p-U_p^{\top })})log(B_q^{-1}{A}_p{e}^{t_p{U}_p^{\top }}{e}^{t_p(U_p-U_p^{\top })})).$$

Calculation leads to the first-order derivative of the function $f\left(A_{p+1}\right)$

$$\begin{array}{c}\hfill \frac{df\left(A_{p+1}\right)}{d{t}_p}=\frac{2}{N}\sum _{q}tr({log}^{\top }(B_q^{-1}{A}_p{e}^{t_p{U}_p^{\top }}{e}^{t_p(U_p-U_p^{\top })})e^{t_p(U_p^{\top }-U_p)}{U}_p{e}^{t_p(U_p-U_p^{\top })}),\end{array}$$

and setting $t_p=0$ yields the $f_{1,p}$

$$f_{1,p}=\frac{2}{N}\sum _{q}tr\left({log}^{\top }\left(B_q^{-1}{A}_p\right)U_p\right)=-tr\left(U_p^{\top }{U}_p\right)\le 0.$$

Under the first-order derivative of the function $f\left(A_{p+1}\right)$ above, direct calculations lead to the second-order derivative of the function $f\left(A_{p+1}\right)$

$$\begin{array}{cc}\hfill \frac{d^{2}f\left(A_{p+1}\right)}{d{t}_p^2}& =\frac{2}{N}\sum _{q}tr(e^{t_p(U_p^{\top }-U_p)}{U}_p^{\top }{U}_p{e}^{t_p(U_p-U_p^{\top })}+{log}^{\top }(B_q^{-1}{A}_p{e}^{t_p{U}_p^{\top }}{e}^{t_p(U_p-U_p^{\top })})\hfill \\ & e^{t_p(U_p^{\top }-U_p)}(U_p^{\top }{U}_p-U_p{U}_p^{\top })e^{t_p(U_p-U_p^{\top })}),\hfill \end{array}$$

and, setting $t_p=0,$ we get

$$\begin{array}{cc}\hfill f_{2,p}& =\frac{2}{N}\sum _{q}tr\left(U_p^{\top }{U}_p+{log}^{\top }\left(B_q^{-1}{A}_p\right)(U_p^{\top }{U}_p-U_p{U}_p^{\top })\right)\hfill \\ & =2tr\left(U_p^{\top }{U}_p\right)-tr\left(U_p^{\top }(U_p^{\top }{U}_p-U_p{U}_p^{\top })\right)\hfill \\ & =2tr\left(U_p^{\top }{U}_p\right)\ge 0.\hfill \end{array}$$

Hence, the nearly-optimal stepsize of the geodesic-based Riemannian-steepest-descent algorithm (19) is

$${\widehat{t}}_p=-\frac{f_{1,p}}{f_{2,p}}=\frac{1}{2}.$$

(20)

Therefore, the iteration algorithm (19) with the stepsize (20) may be recast explicitly as

$$A_{p+1}=A_p{e}^{\frac{1}{N}{\displaystyle \sum _q}{log}^{\top }\left(A_p^{-1}{B}_q\right)}{e}^{\frac{1}{N}{\displaystyle \sum _q}(log\left(A_p^{-1}{B}_q\right)-{log}^{\top }\left(A_p^{-1}{B}_q\right))}.$$

(21)

3.2. Extended Hamiltonian Algorithm on the Lorentz Group

References [16,17] introduced a general theory of extended Hamiltonian (second order) learning on Riemannian manifold, especially, as an instance of learning by constrained criterion function optimization on the matrix manifolds. For the Lorentz group, the extended Hamiltonian algorithm can be expressed by

$$\left\{\begin{array}{c}\dot{A}=X,\hfill \\ \dot{X}=-{\Gamma }_A(X,X)-{\nabla }_{A}f-\mu X,\hfill \end{array}\right.$$

(22)

where $A\in O(n,k)$, $X\in T_{A}O(n,k)$ denotes the instantaneous learning velocity, $f:O(n,k)\to \mathbb{R}$ denotes a criterion function, and the constant $\mu >0$ denotes a viscosity parameter.

Inserting the the expressions of the Christoffel matrix function ${\Gamma }_A$ in Remark 1 and the Riemannian gradient ${\nabla }_{A}f$ (4) into the general extended Hamiltonian system (22), calculations lead to the following expressions:

$$\left\{\begin{array}{c}\dot{A}=AH,\hfill \\ \dot{H}=H^{\top }H-H{H}^{\top }+\frac{1}{2}{I}_{n,k}({\partial }_A^{\top }fA{I}_{n,k}-I_{n,k}{A}^{\top }{\partial }_{A}f)-\mu H.\hfill \end{array}\right.$$

Hence, the second equation may be integrated numerically by a Euler stepping method, while the first one may be integrated via the geodesic. Namely, system (22) may be implemented by

$$\left\{\begin{array}{c}{A}_{p+1}=A_p{e}^{\eta H_p^{\top }}{e}^{\eta H_p-\eta H_p^{\top }},\hfill \\ H_{p+1}=\eta (H_p^{\top }{H}_p-H_p{H}_p^{\top })+(1-\eta \mu )H_p+\frac{1}{2}\eta I_{n,k}({\partial }_{A_p}^{\top }f{A}_p{I}_{n,k}-I_{n,k}{A}_p^{\top }{\partial }_{A_p}f),\hfill \end{array}\right.$$

(23)

where $\eta >0$ denotes the learning rate. The effectiveness of the algorithm is ensured if and only if

$$\sqrt{2{\lambda }_{max}}<\mu <1/\eta ,$$

where ${\lambda }_{max}$ denotes the largest eigenvalue of the Hessian matrix of the criterion function $f\left(A\right)$ (see [17]).

4. Numerical Experiments

In the present section, we give two numerical experiments on the Lorentz group $O(3,1)$ to demonstrate the effectiveness and performance of the proposed algorithms. The Lorentz group as a homogeneous space has four connected components. Now, we study the optimization on a connected component $S{O}_0(3,1)$ containing the identical matrix $I_4$.

In order to emphasize the behavior of the optimization methods RSDA, the EHA and the EGA, in the following experiments, $\eta =0.5$ and $\mu =1.5$ are used.

4.1. Numerical Experiments on Averaging Two Lorentz Matrices

In Section 3.1, it is noticed that the geodesic distance $d(A,B)$ and the discrepancy $D(A,B)$ between two points A and B are different. Now, given $B_1$, $B_2\in O(3,1)$, numerical experiment results are shown to compare $d(B_1,B_2)$ and $D(B_1,B_2)$. The following two points are sought:

$$B_1=\left(\begin{array}{cccc}cosh\left(0.5\right)& 0& 0& sinh\left(0.5\right)\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ sinh\left(0.5\right)& 0& 0& cosh\left(0.5\right)\end{array}\right),B_2=\left(\begin{array}{cccc}cosh\left(0.1\right)& 0& 0& sinh\left(0.1\right)\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ sinh\left(0.1\right)& 0& 0& cosh\left(0.1\right)\end{array}\right).$$

It can be verified that $B_i^{\top }{I}_{3,1}{B}_i=I_{3,1}$, $B_i\in O(3,1)$, $i=1,2$. The optimization problem to compute the geodesic distance $d(B_1,B_2)$ can be solved by means of a descent procedure along the negative Euclidean gradient (9) with a line search strategy. In Figure 1, the results show that the Frobenius norm of U is decreasing steadily when the criterion function $J\left(U\right)$ tends to be constant 0 during iteration. By the numerical experiments, we can get

$$U^{*}=\left(\begin{array}{c}\begin{array}{cccc}0.0000& 0.0000& 0.0000& -0.4000\\ 0.0000& -0.0000& -0.0000& 0.0000\\ 0.0000& -0.0000& -0.0000& 0.0000\\ -0.4000& 0.0000& 0.0000& -0.0000\end{array}\end{array}\right).$$

It is interesting that the geodesic distance between $B_1$ and $B_2$ is $d(B_1,B_2)=\parallel U^{*}\parallel \approx 0.5657$, and the discrepancy between $B_1$ and $B_2$ is $D(B_1,B_2)=\parallel log\left(B_1^{-1}{B}_2\right)\parallel \approx 0.5657$.

In the following experiment, we consider the average value of two points $B_1$ and $B_2$. Figure 2 displays the results of running the RSDA, the EHA, and the EGA. Note that the iteration sequence $\left\{A_p\right\}$ is expected to converge to a matrix that locates amid the samples $\left\{B_q\right\}$, and the condition $\parallel B_q^{-1}{A}_p-I_4\parallel <1$ holds. The initial point is set as $A_0=B_1$, which satisfies $\parallel B_2^{-1}{A}_0-I_4\parallel \approx 0.5921<1$. More initial points can be chosen in a neighbor of the matrix $B_1$ by the following rule (25). It can be found that the averaging algorithms converge steadily and the RSDA has the fastest convergence rate among three algorithms and needs one iteration to obtain the average value as follows:

$$\overline{A}=\left(\begin{array}{c}\begin{array}{cccc}1.0453& 0.0000& 0.0000& 0.3045\\ 0.0000& 1.0000& 0.0000& 0.0000\\ 0.0000& 0.0000& 1.0000& 0.0000\\ 0.3045& 0.0000& 0.0000& 1.0453\end{array}\end{array}\right).$$

However, the EHA and the EGA need 13 and 18 iterations to realize the same accuracy, respectively. It is interesting to compare the discrepancy $D(\overline{A},B_1)\approx 0.2828$ with the discrepancy $D(\overline{A},B_2)\approx 0.2828$, which confirms that the computed average value is truly a midpoint (i.e., it is located at approximately equal discrepancy from the two matrices).

4.2. Numerical Experiments on Averaging Several Lorentz Matrices

The numerical experiments rely on the availability of a way to generate pseudo-random samples on the Lorentz group. Given a point $B\in O(n,k)$ which is referred to as “center of mass” or simply center of the random distribution, it is possible to generate a random sample $B_q\in O(n,k)$ in a neighbor of a matrix B by the rule [26,27]

$$B_q=B{e}^{B^{-1}X},$$

(24)

where the center matrix B can be generated by a curve departing from the identical matrix, $B=e^V$ with V generated randomly matrix in $\mathfrak{o}(n,k)$.

According to the structure of the tangent space $T_{B}O(n,k)$, expression (24) can be written

$$B_q=B{e}^{I_{n,k}(C-C^{\top })},$$

(25)

where $C\in {\mathbb{R}}^{(n+k)\times (n+k)}$ is any unstructured random matrix.

Figure 3 displays the results on optimization over $O(3,1)$ with N = 50 Lorentz matrices. The initial point $A_0$ is randomly generated,

$$A_0=\left(\begin{array}{c}\begin{array}{cccc}0.8412& 0.2117& 0.5487& -0.2351\\ -0.0045& 1.0783& -0.2160& -0.4577\\ -0.5410& 0.2998& 0.8456& -0.3126\\ -0.0203& -0.5454& -0.2508& 1.1665\end{array}\end{array}\right).$$

Figure 3a shows the values of the criterion function (14) through the RSDA, the EHA and the EGA. Numerical experiments show that the average value of 50 Lorentz matrices is

$$\overline{A}=\left(\begin{array}{c}\begin{array}{cccc}0.9127& -0.1134& 0.3933& 0.0235\\ 0.0317& 0.9853& 0.2035& -0.1154\\ -0.4088& -0.1675& 0.8980& -0.0407\\ 0.0342& -0.1087& -0.0504& 1.0077\end{array}\end{array}\right),$$

which can be obtained after seven iterations using the RSDA, while the EHA and the EGA need 15 iterations. Figure 3b gives the norm $\parallel A^{-1}{\nabla }_{A}f\parallel$ of the Riemannian gradient (17) during iteration by the RSDA, from which we can obtain that the norm of the Riemannian gradient tends to be constant 0 after iteration 6. Figure 3c,d give the values of the discrepancy $D(A,B_q)$ before iteration and after iteration, respectively. In particular, note that all the discrepancy from the samples $B_q$ and the matrix A are decreased substantially. The results show that the convergence rate of the RSDA is still the fastest among three algorithms.

5. Conclusions

In this paper, we consider the averaging optimization problem for the Lorentz group $O(n,k)$, namely, to measure the average of a set of Lorentz matrices. In order to tackle the related optimization problem, a geodesic-based Riemannian-steepest-descent algorithm is presented on the Lorentz group endowed with a left-invariant metric (Riemannian metric). The devised averaging algorithm is compared with the extended Hamiltonian algorithm and the Euclidean gradient algorithm. The results of numerical experiments over the Lorentz group $O(3,1)$ show that the convergence rate of the geodesic-based Riemannian-steepest-descent algorithm is the best one among these three algorithms.