# Mutual space-frequency distribution of Gaussian signal

###### Abstract

Mutual space-frequency distribution is proposed and it is shown that Wigner
and Weyl distribution functions are only particular cases of these
distribution. Mutual distribution for Gaussian signal is analytically
obtained. The simple connection between Wigner and Weyl distributions is
established. It is shown that Wigner distribution forms as the rotational
displacement of Weyl distribution on informational diagram of conjugate
coordinates $(x;p)$ on an angle proportional to the mutual parameter $t$. The
results of direct calculations of mutual distribution for Gaussian signal in
the mutual domain are presented.

Yu. M. Kozlovskii

##  1 Introduction

The main aim of space-frequency analysis is elaboration of distributions in
order to get the information about signal simultaneously in coordinate and
frequency domains.As a rule the Fourier transform is used to receive the
signal frequency spectrum. This well-known transformation is good tool for the
analysis of signal intensity distribution in the frequency plane. Such
analysis foresees the calculation of Fourier-spectrum for constant
coordinates. Practically, we have to deal with a certain momentary value of
coordinates for which the signal and their Fourier-image are simultaneously
determinated. With the coordinate variation the appropriate conversion of
Fourier-spectrum also takes place and at once the problem of signal analysis
occurs; the latter contains frequency components that are variable in
accordance with coordinate. In such a case it is important to know the value
of coordinate at which the corresponding transformation of frequency spectrum
takes place. In order to investigate the variation of signal spectrum with the
variation of its coordinate as far back as in 60-80-ies of the previous
century was propound a new approach. It unites the information about
coordinate and frequency constituent of the signal in so-called space-
frequency representations. In such representations under consideration is a
certain mutual function of coordinate and frequency. The idea of construction
of mutual representations originates in works of E.Wigner (1932) [1], D.Gabor
(1946) [2] and J.Weyl (1932) [3]. Before the 80-ies of the previous century
tens of space-frequency representations of such case were taken under
consideration [7, 8, 10]. However, the Wigner and Weyl distributions that are
the most used for the present day remained the prerogative of quantum
mechanics and they have not have precisely expressed use. Only in 1980
T.Klasen and B. Mecklenbrauker worked out the theory of application of Wigner
distribution for space-frequency analysis of signals. Its main results were
published in the series of works under the title ”Wigner distribution - the
instrument for space-frequency analysis of signals” [4, 5, 6]. The successful
use of the Wigner distribution in the theory of signals was stipulated by its
”good” mathematical characteristics, especially by its representative
characteristics that is basic in the restoration of signal intensity
distribution.

Within the framework of the given investigation among the variety of space-
frequency representations we single out two basic distributions by Wigner and
Weyl that are widely used by the theory of signals in solving inverse physical
problems [8, 11]. The investigation of the signal characteristics takes place
on the basis of comparison with its displaced analogues. The shift within a
time results in subtraction of a specific value from the signal argument

$\rightarrow x_{\tau}(t)=x(t+\tau).$ (1)

The suitable displacement in accordance with frequency results in displacement
of argument of the Fourier-spectrum signal, what equals to multiplication by
phase multiplier in coordinate plane.

${\bf X_{\omega}}\rightarrow X_{\omega}(t)=x(t)e^{i\omega t}.$ (2)

Similar correlations are well-known from classical analysis [11]. Within the
limits of space-frequency analysis we are interested in the signals displaced
simultaneously with time and frequency, namely

$x_{-\frac{\tau}{2},-\frac{\omega}{2}}=x\left(t-\frac{\tau}{2}\right)e^{-i\omega
t/2},$ (3)
$x_{\frac{\tau}{2},\frac{\omega}{2}}=x\left(t+\frac{\tau}{2}\right)e^{i\omega
t/2}.$ (4)

We can easily calculate the value of displacement between the signals

$d({\bf x_{-\frac{\tau}{2},-\frac{\omega}{2}}},{\bf
x_{\frac{\tau}{2},\frac{\omega}{2}}})^{2}=2||x||^{2}-2\Re\left\\{{\cal
A}_{xx}(\omega,\tau)\right\\}.$ (5)

${\cal A}_{xx}(\omega,\tau)$ in this correlation plays a part of the distance
and is called a time-frequency autocorrelation function or the ambiguity
function

${\cal A}_{xx}(\omega,\tau)=\int\limits
x^{*}\left(t-\frac{\tau}{2}\right)x\left(t+\frac{\tau}{2}\right)\exp{(-i\omega
t)}dt.$ (6)

In accordance to the Parceval theory we can rewrite this correlation by
Fourier-images of displaced signals

${\cal A}_{xx}(\omega,\tau)=\frac{1}{2\pi}\int\limits
X\left(\nu-\frac{\omega}{2}\right)X^{*}\left(\nu+\frac{\omega}{2}\right)\exp{(i\omega\tau)}d\omega.$
(7)

The function in such a form was for the first time set by J. Weyl [3] and for
the present time is known in the theory of signals under the name of Weyl
distribution or the ambiguity function. Having realized the direct and inverse
Fourier-transform we get the value

${\cal
W}_{xx}(t,\nu)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}{\cal
A}_{xx}(\omega,\tau)e^{i\omega t}e^{-i\nu\tau}d\omega d\tau.$ (8)

that in explicite form has the following notation

${\cal W}_{xx}(t,\nu)=\int\limits
x^{*}\left(t-\frac{\tau}{2}\right)x\left(t+\frac{\tau}{2}\right)\exp{(-i\nu\tau)}d\tau,$
(9)

or by the Fourier-image of the function ${x}(t)$

${\cal W}_{xx}(t,\nu)=\frac{1}{2\pi}\int\limits
X\left(\nu-\frac{\omega}{2}\right)X^{*}\left(\nu+\frac{\omega}{2}\right)\exp{(i\omega
t)}d\omega.$ (10)

The function was firstly introduced by E.Wigner and is named after him - the
Wigner function of distribution or just Wigner distribution [1].

Joint space-frequency represenations are widely used not only in the theory of
signals. They have a number of practical use in different fields of physics,
geology, seismology, etc. Within the limits of the given research we are
interested in the use of such distributions in the region of representations
treatment and recognition of images. This field of the physics imposes a set
of demands that the generalized distribution of signals have to meet. For the
efficient use in the theory of representatives of space-frequency distribution
they have to be characterized by representative property; for the limiting
values of the variable $t$ they have to develop into known distributions; to
have high distributive capacity in the field of Wigner distribution as well as
in Weyl distribution; to take positive values.

Representations that the most precisely meet the demands stated in the theory
of image are basic distributions by Wigner and Weyl. These distributions also
have their own peculiarities. Unfortunately, up to the present there exists no
simple deduction about the expediency of use of this or that distribution.
There exist a number of signals for which Wigner and Weyl distributions
proceed into the negative region. In such cases these distributions are
interpreted as quasiprobable. In spite of the external resemblance of the
properties of Wigner and Weil distributions they have the peculiarity of
principle in the mechanism of renewal of the entrance signal according to the
known distribution. Wigner formalism allows renewing the signal according to
the so-called marginal distribution

$|x(t)|^{2}=\int\limits_{-\infty}^{\infty}{\cal W}_{xx}(t,\nu)d\nu.$ (11)

In the frame of Weyl formalism signal reconstruction takes place using signal
restoration scheme

$|x(t)|^{2}=\int\limits_{-\infty}^{\infty}{\cal A}_{xx}(0,\omega)e^{i\omega
t}d\omega.$ (12)

Sometimes working with the scheme of renewal according to the Weyl
distribution is to the great extent more easily (due to the integrating of its
crosscut) than in the case of marginal distribution. For the present Wigner
distribution is more commonly applied as it uses marginal distributions that
can be measured by experiment. Though, both approaches have the right to
existence.

Within the framework of the given experiment we investigate uninterrupted
transition between these distributions by means of introduction of a
generalized common function of time and frequency that depends from a variable
$t$.

For present day various types of space-frequency distributions are
successfully used for the analysis of nonstationary signals [7, 8]. Many of
such distributions are characterized by advantages as well as by disadvantages
in use in various fields of physics. Wigner and Weyl distributions are widely
used in space-frequency analysis and in particular in optical information
processing systems. Well known is a fact, that the very distributions posess
such charateristics, that are successfully used for description of many
optical systems. The spectrum of appliance of these distributions is extremely
wide. They are used, in particular, in the theory of optical lens system,
theory of communication, hydrolocation and other fields [9, 10, 11]. The
researches of last year proved the efficiency of use of space-frequency
distributions in biology and medicine; especially Wigner distribution was
successfully used for renewal of volumetric structure of objects within the
framework of optical tomography [12, 13, 14]. One of the promising
investigation directions within the space-frequency processing of signals is
studying the properties of novel space-frequency representations of the
distributions, with the aim of their further applications in different areas
of physics and medicine. Unfortunately, it often happens that some space-
frequency distributions do not meet demands raised by one or another specific
application. In this relation, many of the existing distributions need
generalization or improvements when applied to a given problem. During the
second half of the past century and the beginning of this one, a cleartendency
has been observed towards generalization of different space-frequency
distributions. The first attempt of such a generalization has been due to
L.Cohen [15] as long ago as in 1966. The author has introduced a number of
quasiprobable distributions that provide proper quantum mechanical marginal
distributions. Within the limits of this research the Wigner distribution was
examined as a separate case. The next step has been done by N. De Brujin in
1973 [16]. His work has been devoted to elaboration of theory of generalized
functions, with application concerned with Wigner and Weil distributions.
Summarizing the results of numerous investigations L.Cohen [7, 10] has
suggested to has suggested a generalized distribution involving a certain
kernel

$\displaystyle C(x,\omega,\Phi)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}f\left(y+\frac{x_{0}}{2}\right)f^{*}\left(y-\frac{x_{0}}{2}\right)$
(13) $\displaystyle\times$ $\displaystyle\Phi(x,\omega)e^{-i(\omega
x_{0}-\omega_{0}x+\omega_{0}y)}dydx_{0}d\omega_{0}.$

Depending on the form of the kernel $\Phi(x,\omega)$ , this distribution
degenerates into one of the known distributions (Wigner, Weyl, Woodward,
Kirkwood, Page, Mark, etc.).

In the general case expression (13) describes the class of space-frequency
distribution later named Cohen’s class. Members of this class are known
distributions as well as a set of still unknown distributions, that also
satisfy all necessary requirements of existing distributions. The theory of
generalization of space-frequency distributions has been developed by also
famous specialist in the theory of signals by A. Mertins. In his monograph
”Signal Analysis”[11] he has singled out thistopic into a separate section
”General space-frequency distributions”. The author has stated that the Wigner
distribution serves as an excellent tool for space-frequency analysis as long
as a linear dependence is kept between the instantaneous coordinates
andfrequencies. Otherwise, a need in generalizing appears, whose general
principles are described in detail in the mentioned work.

Among numerous recent studies related to generalizing space-frequency
distributions, we should mention only the most typical ones. The PhD Thesis by
L. Durak ”Novel time-frequency analysis technique for deterministic signals”
[17] is one of such studies, where a close attention has been paid to
generalizing distributions and introducing their additional parameters.

Different types of generalizations of space-frequency distributions have been
thoroughly considered in the book by B. Boashash [8]. Among a number of
studies included in it, we would like to emphasize the works by R. Baraniuk
(p. 123), X.Xia (p. 223) and A. Papandreou-Suppappola (p. 643). Within the
mentioned collection, the work by G. Matz by F. Hlawatsch (p. 400) is of
particular interest in relation to the problem of distributions
generalization. It considers methodology for constructing generalized
distributions on the basis of both the Wigner distribution and the ambiguity
function (i.e. the Weyl distribution).

In the present work we try to use interlinks between the Wigner and Weyl
distributions with the aim of joining them into a single, more general
distribution. Up to date, it has been revealed that the two distributions are
related by a double Fourier transform. The results obtained by us allow
tracing transformation of one of thedistributions into the other, while
changing the distribution parameter $t$. This generalized distribution
generates a whole set of new distributions formed in the process of switching
between the basic distributions. The latter fact may be important from the
viewpoint of possible practical applications. For the present day a choice
between the Weyl and Wigner distributions remains ambiguous. Each of them has
its own scheme for reconstruction of signal intensity distribution. The scheme
adopted for the Wigner distribution includes calculating the marginal
distributions [7]. The Weyl distribution provides much simpler reconstruction
scheme,owing to simpler mathematical transformations [18, 22]. Traditionally,
the Wigner distribution has been used in a large majority of studies performed
within the field. Introduction of the mutual distribution would mean a
possibility for calculating ’mixed’states and determining the necessary
contributions of each of the limiting distributions. As stated above, there
appears a possibility for generalization of distributions concerning various
applied problems. However, only S. Chountasis has suggested the approach [19]
that enables transitions between the Wigner and Weyl distributions. Such
distributions play an important role in the analysis of phase space and,
moreover, can be immediately applied in the Wigner tomography [12, 14]..In
1999 S. Chountasis and co-authors [19, 21] have developed a general
distribution based on the Wigner formalism, which involves an additional
parameter $\theta$. This study has been performed in frame of quantum-
mechanical formalism. It allows passing the Wigner and Weyl distributions into
each other by means of changing the common parameter.

The problem of calculation of a classical analogue of this generalized
distribution remains urgent. It may be constructed based on the results [19]
or using the formalism of Weyl distribution, as has been done by the present
author when studying the properties of fractional Fourier transform [23].
Similarly to the works [18, 22], the author has employed peculiarities of
reconstruction of signal intensity based upon the Weyl distribution.
Meanwhile, it is just this reconstruction scheme is realized experimentally in
the real optical schemes [23].

##  2 Mutual distribution: basic relations

###  2.1 Theoretical statements

In the present work we propose to use a type of generalized distribution with
paramether $t$ based upon the Weyl distribution function. The use of Weyl
distribution has a peculiarity of principle in comparison with the function of
Wigner distribution, that consists of the possibility to renew the intensity
of distribution. The latter is registered experimentally at the output of the
optical system. Common distribution of two signals $f_{1}(x)$ and $f_{2}(x)$
may be written as follows [25]

$\displaystyle{\cal K}_{f_{1}f_{2}^{*}}^{(t)}(x;p)$ $\displaystyle=$
$\displaystyle\frac{C_{t}}{1+t}\int\int
dx_{0}d\omega_{0}\exp\left\\{i\left[x_{0}p-\omega_{0}x\right]\right\\}$ (14)
$\displaystyle\times$
$\displaystyle\exp\left\\{-i\frac{(x-x_{0})^{2}+(p-\omega_{0})^{2}}{tan(\theta/\
2)}\right\\}$ $\displaystyle\times$ $\displaystyle\int\limits
f_{1}\left(z+\frac{x_{0}}{2}\right)f_{2}^{*}\left(z-\frac{x_{0}}{2}\right)\exp{(-i\omega_{0}z)}dz.$

Constant $C_{t}$ and variable of the generalized distribution $t$ are
determined by the expressions

$C_{t}=\frac{2}{\pi}\frac{1}{1-\exp{i\theta}},\qquad t=\frac{\theta}{\pi}.$
(15)

Distribution (14) is called the mutual space-frequency distribution, or
concisely the mutual distribution. Limiting cases of ditribution (14) is Weyl
distribution (ambiguity function) (6) with the value of the parameter $t=0$
and Wigner distribution (9) with the value of parameter $t=1$. Accordingly,
two known distributions (6) and (9) have plenty of alternative distributions
and to each of them corresponds a specific value of parameter $t$.

The expression of the mutual distribution (14) may be also represented in the
simplified form using Weyl distribution (6)

$\displaystyle{\cal K}_{f_{1}f_{2}^{*}}^{(t)}(x;p)$ $\displaystyle=$
$\displaystyle\frac{C_{t}}{1+t}\int\int dx_{0}d\omega_{0}{\cal
A}_{f_{1}f_{2}^{*}}(x_{0};\omega_{0})$ (16) $\displaystyle\times$
$\displaystyle\exp\left\\{i\left[x_{0}p-\omega_{0}x\right]\right\\}\exp\left\\{-i\frac{(x-x_{0})^{2}+(p-\omega_{0})^{2}}{tan(\theta/\
2)}\right\\}.$

Performing the converted transformation we can render the Weyl distribution by
means of the above introduced function of the mutual distribution

$\displaystyle{\cal
A}_{f_{1}f_{2}^{*}}(x_{0}^{{}^{\prime}};\omega_{0}^{{}^{\prime}})$
$\displaystyle=$ $\displaystyle\frac{1+t}{C_{t}}\int\int dxdp{\cal
K}_{f_{1}f_{2}^{*}}^{(t)}(x;p)$ (17) $\displaystyle\times$
$\displaystyle\exp\left\\{-i\left[x_{0}^{{}^{\prime}}p-\omega_{0}^{{}^{\prime}}x\right]\right\\}\exp\left\\{i\frac{(x-x_{0}^{{}^{\prime}})^{2}+(p-\omega_{0}^{{}^{\prime}})^{2}}{tan(\theta/\
2)}\right\\}.$

Formula (17) constitutes the inverse connection between the simple ${\cal
A}_{f_{1}f_{2}^{*}}(x_{0};\omega_{0})$ and generalized ${\cal
K}_{f_{1}f_{2}^{*}}^{(t)}(x;p)$ Weyl distributions. This makes the possibility
of restoring the distribution of signal intensity according to the mutual
distribution what has not been established when using the generalized Wigner
distribution [19].

###  2.2 Representation in the terms of Wigner distribution

In order to make a comparison of the results with the existing analogues it is
indispensable to have possibility to precisely calculate the limiting cases of
the mutual distribution (17). Calculation of the limiting case $t=1$ by means
of the formula (16) may be conducted precisely, and in the case $t=0$ such
transition is not a trivia matter. To make the calculations simplier we
introduce the representation of the mutual space-frequency distribution by
means of Wigner distribution. We make use of the known identity

$\displaystyle{\cal
A}_{f_{1}f_{2}^{*}}(x_{0};\omega_{0})=\frac{1}{2\pi}\int\int d\xi d\eta{\cal
W}_{f_{1}f_{2}^{*}}(\eta;\xi)\exp(-i\omega_{0}\eta)\exp(i\xi x_{0}),$ (18)

that connects Weyl and Wigner distributions.

Placing (18) into the expression (16) and making a set of conversions we
receive a formula describing the mutual space-frequency distribution in the
terms of Wigner distribution

$\displaystyle{\cal{K}}_{f_{1}f_{2}^{*}}^{(t)}(x;p)$ $\displaystyle=$
$\displaystyle\widetilde{C}_{t}\int\int dx_{0}d\omega_{0}{\cal
W}_{f_{1}f_{2}^{*}}(x_{0};\omega_{0})\exp\left\\{-i\left[+x_{0}p+\omega_{0}x\right]\right\\}$
(19) $\displaystyle\times$
$\displaystyle\exp\left\\{i\frac{1}{4}tan\frac{\theta}{2}\left[(x-x_{0})^{2}+(p-\omega_{0})^{2}\right]\right\\},$

where the constants $\widetilde{C}_{t}$ are determined by means of the
correlation

$\widetilde{C}_{t}=\frac{-i}{\pi}\frac{1}{1-\exp{i\theta}}\tan\frac{\theta}{2}\frac{1}{1+t}.$
(20)

As it can be easily seen in the case of representation of the mutual
distribution by means of the Wigner distribution peculiarities in the point
$t=0$ and around it dissappear, however, there appear peculiarities around the
point $t=1$. Thereby the pair of representations: (16) and (19) complement one
another and completely describe mutual space-frequency distribution in the
region $t=[0,1]$.

###  2.3 Limiting cases

The aim of this work is to determine the mechanism of re-distribution between
Wigner and Weyl distributions and to investigate the peculiarities of mutual
distributions describing the region of values $t=[0,1]$. Consequently the
investigation of limiting cases (16) and (19) of mutual distribution is of
peculiar importance. Let us study the limiting cases.

_Case $t=1$._

To describe this case we make use of coordinate representation of the mutual
distribution (16). Placing $t=1$ (or $\theta=\pi$) under (16) we arrive to the
following result

$\displaystyle{\cal K}_{f_{1}f_{2}^{*}}^{t=1}(x;p)=\frac{1}{2\pi}\int\int
dx_{0}d\omega_{0}{\cal
A}_{f_{1}f_{2}^{*}}(x_{0};\omega_{0})\exp(ix_{0}p-i\omega_{0}x).$ (21)

Accordingly to (8) we have

$\displaystyle{\cal K}_{f_{1}f_{2}^{*}}^{t=1}(x;p)={\cal
W}_{f_{1}f_{2}^{*}}(x;p).$ (22)

The limiting case $t=1$ of mutual space-frequency distribution corresponds to
the function of Wigner distribution.

_Case $t=0$._

To describe this case we make use of coordinate representation of the mutual
distribution (19). Placing here $t=0$ (or $\theta=0$) we arrive to the
following result

$\displaystyle{\cal K}_{f_{1}f_{2}^{*}}^{t=0}(x;p)=\frac{1}{2\pi}\int\int
dx_{0}d\omega_{0}{\cal
W}_{f_{1}f_{2}^{*}}(x_{0};\omega_{0})\exp(-i\omega_{0}x+ix_{0}p).$ (23)

Accordingly to the formula connecting Wigner and Weyl distributions we obtain

$\displaystyle{\cal K}_{f_{1}f_{2}^{*}}^{t=0}(x;p)={\cal
A}_{f_{1}f_{2}^{*}}(x;p).$ (24)

Expression (24) is identical with the function of Weyl distribution. Hereby,
introduced by us distribution (14) in the values of the limiting cases of the
parameter $t$ describes known distributions (6) and (9). This distribution has
two equivalent representations: by means of the function of Wigner
distribution (19) or function of Weyl distribution (16). Having at least one
from the basic functions of distribution we can obtain the image of mutual
distribution. The object of the further investigation is to study the
peculiarities of the set of intermediate distributions when $(0<t<1)$. In
order to display the peculiarities of the mutual distribution and to visually
demonstrate the results we illustrate it on the example of Gaussian signal.
This is one of the simplest types of signals that allows us to make
calculations in the explicit form. We shall notice that similar calculations
may be also made with other signals (in peculiar with orthogonal impulse,
etc.). Choice of the Gaussian signal is related to explicit form of the common
distribution. Herewith we can trace the mechanism of transition of Weyl
distribution into Wigner distribution and vice versa.

##  3 Mutual distribution of Gaussian signal

###  3.1 Basic relations

In the work [24] have been investigated the main peculiarities of generalized
coordinate-frequency distribution on the example of Gaussian signal. Obtained
results were based on the frequent calculations of the proper distributions
that allow to properly evaluate processes of distribution. In the present work
we conduct an analytical calculation of the mutual space-frequency
distribution on the example of Gaussian signal that is represented in the
following way

$g(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp{\left(-\frac{x^{2}}{2\sigma^{2}}\right)}.$
(25)

Fourier image of these function is

$\hat{{\cal
F}}[g(x)]=\int\limits_{-\infty}^{\infty}g(x)e^{ixp}dx=\exp{\left(-\frac{p^{2}\sigma^{2}}{2}\right)}.$
(26)

Such selection form of Gaussian functions are stipulated by the fact that its
integrating according to all values $x$ results in value

$\int\limits_{-\infty}^{\infty}g(x)dx=1.$ (27)

As investigations of the mutual space-frequency distribution provides the
study of re-distribution between Wigner and Weyl distributions it is
reasonable to present the explicit form of these distributions for the case of
Gaussian signal [24].

_Weyl distribution for Gaussian signal_

${\cal
A}_{ff^{*}}(x_{0};\omega_{0})=\frac{1}{2\sqrt{\pi}\sigma}\exp{\left(-\frac{x_{0}^{2}}{4\sigma^{2}}-\frac{\omega_{0}^{2}\sigma^{2}}{4}\right)}.$
(28)

_Wigner distribution for Gaussian signal_

${\cal
W}_{ff^{*}}(x;\omega)=\frac{1}{\sqrt{\pi}\sigma}\exp{\left(-\frac{x^{2}}{\sigma^{2}}-\omega^{2}\sigma^{2}\right)}.$
(29)

It is known [7], that one of the basic characteristics of the signal is its
representativity that is the possibility to restoring the signal according to
its distribution. Schemes of renewing for Weyl and Wigner distributions are
well known. For the case of Gaussian signal they have the following form

_restoration scheme for Gaussian signal by Weyl distribution_

$|g(x)|^{2}=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}{\cal
A}_{ff^{*}}(0;\omega_{0})e^{i\omega_{0}x}d\omega_{0}=\frac{1}{2\pi\sigma^{2}}\exp{\left(-\frac{x^{2}}{\sigma^{2}}\right)}.$
(30)

_restoration scheme for Gaussian signal by Wigner distribution_

$|g(x)|^{2}=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}{\cal
W}_{ff^{*}}(x;\omega)d\omega=\frac{1}{2\pi\sigma^{2}}\exp{\left(-\frac{x^{2}}{\sigma^{2}}\right)}.$
(31)

Independently from the choice of distribution form Gaussian signal is
precisely renewed using both schemes of renewing.

We find explicit analytical form of the mutual distribution for the case of
Gaussian signal. For this reason we will place into the formula of mutual
distribution (16) the Weyl distribution of Gaussian signal (16). We obtain
expression

$\displaystyle{\cal
K}_{ff^{*}}^{t}(x;p)=\frac{C_{t}}{1+t}\frac{1}{2\sigma\sqrt{\pi}}\exp\left\\{-\frac{i}{T}(x^{2}+p^{2})\right\\}I_{1}(x,p)I_{2}(x,p),$
(32)

where

$\displaystyle I_{1}(x,p)=\int
dx_{0}\exp\left\\{ix_{0}p+\frac{i}{T}2xx_{0}-\frac{i}{T}x_{0}^{2}-\frac{x_{0}^{2}}{4\sigma^{2}}\right\\},$
(33) $\displaystyle I_{2}(x,p)=\int
d\omega_{0}\exp\left\\{-i\omega_{0}x+\frac{i}{T}2p\omega_{0}-\frac{i}{T}\omega_{0}^{2}-\frac{\omega_{0}^{2}\sigma^{2}}{4}\right\\}.$
(34)

Expressions (33) and (34) may be written in the following form

$\displaystyle
I_{1}(x,p)=\frac{\sqrt{\pi}}{\sqrt{a_{1}}}\exp\left\\{-\frac{1}{a_{1}}\left(\frac{x}{T}+\frac{p}{2}\right)^{2}\right\\},$
(35) $\displaystyle
I_{2}(x,p)=\frac{\sqrt{\pi}}{\sqrt{a_{2}}}\exp\left\\{-\frac{1}{a_{2}}\left(\frac{p}{T}-\frac{x}{2}\right)^{2}\right\\},$
(36)

where the following symboles are introduced

$\displaystyle a_{1}=\frac{T+i4\sigma^{2}}{4T\sigma^{2}},\qquad
a_{2}=\frac{T\sigma^{2}+4i}{4T},\qquad T=\tan\left(\frac{\theta}{2}\right).$
(37)

The mutual distribution of Gaussian signal (32) has the form

$\displaystyle{\cal K}_{ff^{*}}^{t}(x;p)$ $\displaystyle=$ $\displaystyle
C_{t}^{\sigma}\exp\left\\{-\frac{i}{T}\left(x^{2}+p^{2}\right)\right\\}$ (38)
$\displaystyle\times$
$\displaystyle\exp\left\\{-\frac{1}{a_{1}}\left(\frac{x}{T}+\frac{p}{2}\right)^{2}\right\\}\exp\left\\{-\frac{1}{a_{2}}\left(\frac{p}{T}-\frac{x}{2}\right)^{2}\right\\}.$

where constant

$C_{t}^{\sigma}=\frac{C_{t}}{1+t}\frac{\sqrt{\pi}}{2\sigma}\frac{1}{\sqrt{a_{1}a_{2}}}$
(39)

depends from dispersion of Gaussain distribution of signal $\sigma$ and from
the values $a_{1}$, $a_{2}$. From the values $a_{1}^{-1}$ and $a_{2}^{-1}$
depend also expressions in the index of the expression exponential curve (41).
We depict them in the form

$a_{1}^{-1}=a_{11}+ia_{12}:\qquad
a_{11}=\frac{4T^{2}\sigma^{2}}{16\sigma^{4}+T^{2}},\qquad
a_{12}=-\frac{16T\sigma^{4}}{16\sigma^{4}+T^{2}},$ (40)
$a_{2}^{-1}=a_{22}+ia_{21}:\qquad
a_{22}=\frac{4T^{2}\sigma^{2}}{16+T^{2}\sigma^{4}},\qquad
a_{21}=-\frac{16T}{16+T^{2}\sigma^{4}}.$ (41)

Placing (40) and (41) under (38) we find explicit analytical form of the
mutual distribution of the Gaussian signal (14)

$\displaystyle{\cal K}_{ff^{*}}^{t}(x;p)$ $\displaystyle=$
$\displaystyle\frac{C_{t}}{1+t}\frac{\sqrt{(}\pi)}{2\sigma}\frac{1}{\sqrt{a_{1}a_{2}}}\exp\left\\{-\frac{i}{T}\left(x^{2}+p^{2}\right)\right\\}$
(42) $\displaystyle\times$
$\displaystyle\exp\left\\{-a_{11}\left(\frac{x}{T}+\frac{p}{2}\right)^{2}-ia_{12}\left(\frac{x}{T}+\frac{p}{2}\right)^{2}\right\\}$
$\displaystyle\times$
$\displaystyle\exp\left\\{-a_{22}\left(\frac{p}{T}-\frac{x}{2}\right)^{2}-ia_{21}\left(\frac{p}{T}-\frac{x}{2}\right)^{2}\right\\}.$

We may check that expression (42) in limiting cases $t=0$ and $t=1$ changes
into Weyl and Wigner distributions respectively. In order to ascertain the
circumstance we investigate the conduct of the expression (42) in the limit
$t\rightarrow 1$. Parameter $t$ is determined according to (15) by value
$\theta$. We introduce small value $\alpha$ ($\alpha\ll 1$) and depict
$\theta$ in the form

$\theta=\pi-\alpha.$ (43)

Region of small values $\alpha$ corresponds to quasi-Wigner region of mutual
function of distribution. Then the value $C_{t}$ that is a part of
$C_{t}^{\sigma}$ in case of small values $\alpha$ can be depicted as follows

$C_{t}\approx\frac{1}{\pi}\left(1+i\frac{\alpha}{2}\right).$ (44)

Asymptotics of the values $(a_{1}a_{2})^{-1/2}$ may be easily find when use
its module $r_{12}$ and argument $\varphi_{12}=\varphi_{1}+\varphi_{2}$

$(a_{1}a_{2})^{-1/2}=r_{12}e^{i(\varphi_{1}+\varphi_{2})/2},$ (45)

where

$\displaystyle
r_{12}=\pm\frac{4T\sigma}{(16+T^{2}\sigma^{4})^{1/4}(16\sigma^{4}+T^{2})^{1/4}},$
$\displaystyle\varphi_{1}=\arctan\left(-\frac{4\sigma^{2}}{T}\right),\quad\varphi_{2}=\arctan\left(-\frac{4}{T\sigma^{2}}\right).$
(46)

In the case of small values $\alpha$ ($T\gg 1$) we have such approximate form
for the value $r_{12}$

$r_{12}=4\left(1-\frac{4}{T^{2}}\left(\sigma^{4}+\sigma^{-4}\right)\right),$
(47)

and for value $\varphi_{1}$ and $\varphi_{2}$ we find

$\varphi_{12}=\varphi_{1}+\varphi_{2}=-\frac{4}{T}\left(\sigma^{2}+\sigma^{-2}\right).$
(48)

Thereby, the constant $C_{t}^{\sigma}$ from (42) in quasi-Wigner region
$\alpha\ll 1$ has the form

$C_{t}^{\sigma(1)}=\frac{1}{\sigma\sqrt{\pi}}\left(1+\frac{i}{2}\alpha\right)\left(1-\frac{4}{T^{2}}\left(\sigma^{4}+\sigma^{-4}\right)\right)e^{-\frac{2i}{T}\left(\sigma^{2}+\sigma^{-2}\right)}.$
(49)

Index 1 in the value $C_{t}^{\sigma}$ denotes a condition $\alpha\ll 1$.
Naturally, that in the $\alpha\rightarrow 0$ we has the value

$C_{t}^{\sigma}=\frac{1}{\sigma\sqrt{\pi}},$

that precisely corresponds to the amplitude of the expression for Wigner
distribution (29).

Let us observe the coefficient $C_{t}^{\sigma}$ from (42) in the case of small
values $\theta$ $(t\rightarrow 0)$. We shall call the region of parameter $t$
values quasi-Weyl region as far as in the value $t=0$ we aquire Weyl
distribution. Similarly as in the case $t=1$, we find for small values the
following expressions

$C_{t}=\frac{i}{\pi T}(1-iT),\qquad
r_{12}=T,\qquad\varphi_{1}=\varphi_{2}=-\pi/2,\qquad
C_{t}^{\sigma(2)}=\frac{1}{2\sqrt{\pi}\sigma}.$ (50)

Such value of constant $C_{t}^{\sigma(2)}$ exactly corresponds to the
amplitude of the value from Weyl distribution of Gaussian signal (28). Thus,
the mutual distribution of Gaussian signal in limiting cases according to
amplitude precisely coincides with the known distributions of Wigner and Weyl.

Calculation of real and imaginary part of constant $C_{t}^{\sigma}$ depicted
on the Fig.1. for $C_{t}^{\sigma}$ provides an exact correspondence according
to amplitude with Wigner and Weyl distributions. In the common region appears
imaginary part $C_{t}^{\sigma}$ that is inherent only in mutual distribution.
In the limiting cases $t=0(\theta=0^{0})$ and $t=1(\theta=180^{0})$ imaginary
part dissappears what corresponds to the cases of basic distributions. Worth
mentioning is also peculiar conduct of the imaginary part of constant
$C_{t}^{\sigma}$ having maximum in the quasi-Wigner region.

Let us proceed the investigation of limiting cases of expessions placed in the
index of the exponent on a curve formula (42). As was shown above the
amplitude of mutual distribution in limiting cases coincides with the
amplitude of known distributions of Weyl and Wigner. Providing that the form
of these dirtibutions will also coincide functionally (as functions $x$ and
$p$), mutual distribution may be considered as generalization of well-known
distributions of Wigner and Weyl. Having made a number of mathematical
transformations the expression (42) aquires the following form

$\displaystyle{\cal K}_{rr^{*}}^{(t)}(x;p)$ $\displaystyle=$ $\displaystyle
C_{t}^{\sigma}\exp\left\\{-gx^{2}-fp^{2}-dxp\right\\},$ (51)

where coefficients $g$, $p$ and $d$ are certain complex values

$g=g_{0}+ig_{1},\quad f=f_{0}+if_{1},\quad d=d_{0}+id_{1}$ (52)

moreover

$\displaystyle g_{0}=\frac{\sigma^{2}}{\mu}(64+20T^{2}\sigma^{4}+T^{4}),\qquad
g_{1}=\frac{T}{\mu}(16[1-4\sigma^{4}]+T^{2}[\sigma^{4}-4]);$ $\displaystyle
f_{0}=\frac{\sigma^{2}}{\mu}(64\sigma^{4}+20T^{2}+T^{4}\sigma^{4}),\qquad
f_{1}=\frac{T\sigma^{4}}{\mu}(16[\sigma^{4}-4]+T[1-4\sigma^{4}]);$
$\displaystyle d_{0}=\frac{4T\sigma^{2}}{\mu}(1-\sigma^{4})(16-T^{2}),\qquad
d_{1}=16T^{2}\frac{1-\sigma^{8}}{\mu}.$ (53)

where the value $\mu$ has the form

$\displaystyle\mu=16^{2}\sigma^{4}+16T^{2}(1+\sigma^{8})+\sigma^{4}T^{4}$ (54)

Constant $C_{t}^{\sigma}$ from (39) may be represented as

$\displaystyle
C_{t}^{\sigma}=\frac{i}{\pi}\frac{1+e^{-i\theta}}{\sin\theta}\frac{1}{1+t}\frac{\sqrt{\pi}}{2\sigma}r_{12}e^{i(\varphi_{1}+\varphi_{2})/2},$
(55)

where

$\displaystyle r_{12}=\frac{4\sigma
T}{(16+T^{2}\sigma^{4})^{1/4}(16\sigma^{4}+T^{2})^{1/4}},$
$\displaystyle\varphi_{1}=\arctan\left(-\frac{4\sigma^{2}}{T}\right),\qquad\varphi_{2}=\arctan\left(-\frac{4}{T\sigma^{2}}\right)$
(56)

Representation (51) is an explicit form of the mutual distribution of Gaussian
signal. Let us consider asymptotic values of $g$, $f$ and $d$ in the limiting
cases $t=0$ and $t=1$.

In the case $t=0$ ($T\rightarrow 0$) from (53), (54) we find

$g_{0}=(4\sigma^{2})^{-1},\quad f_{0}=\sigma^{2}/4.$ (57)

The rest of coefficients turns into zero

$g_{1}=f_{1}=d_{0}=d_{1}=0.$ (58)

In the limiting case $t=1$ ($T\rightarrow\infty$) we receive

$g_{0}=\sigma^{-2},\quad f_{0}=\sigma^{2}.$ (59)

Other coefficients dissapear as in (58). We shall notice that imaginary part
of cross-expressions in (51) real part dissapears in $t\rightarrow 0$ as well
as in $t\rightarrow 1$.

What concerns the values $g$ and $f$, their imaginary parts also tend to zero
in the case $t=0$ and $t=1$. Thus, the very quadratical members are
responsible for forming of the distribution in limiting cases as they are
included into known Wigner and Weyl distributions, and cross-representations
that are inherent only in mutual distribution dissapear.

##  4 Signal restoration scheme in mutual distribution domain

As it is known, restoration of the signal according to Weyl distribution takes
place correspondingly to (30), and according to Wigner distribution -
accordning to (31). Above mentioned schemes of signal restoration may be
united into one formula using the mutual space-frequency distribution
suggested above. We shall introduce the value

$f_{\theta}(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\rho
w_{\theta}\left(x\sin(\theta/2)p\right)e^{ipx\cos(\theta/2)}.$ (60)

It can be easily seen, that when $\theta=0$ we arrive at Weyl renewal scheme,
and when $\theta=\pi$ \- we aquire Wigner renewal signal. Taking into
consideration (51), for the function $w_{\theta}$ we have

$\displaystyle
w_{\theta}\left(x\sin(\theta/2),p\right)=C_{t}^{\sigma}\exp\left\\{-g\sin^{2}(\theta/2)x^{2}-fp^{2}-d\sin(\theta/2)xp\right\\}.$
(61)

In the result of integration of (60) we aquire

$\displaystyle f_{\theta}(x)$ $\displaystyle=$
$\displaystyle\frac{C_{t}^{\sigma}}{2\pi}\frac{\sqrt{\pi}}{(f_{0}+if_{1})^{1/2}}e^{-x^{2}\sin\theta/2(g_{0}+ig_{1})}$
(62) $\displaystyle\times$
$\displaystyle\exp\left\\{\frac{x^{2}}{4}\frac{1}{(f_{0}+if_{1})}\left(\sin\frac{\theta}{2}(d_{0}+id_{1})-i\cos\theta/2\right)^{2}\right\\}.$

It can be easily seen that in limiting cases function $f_{\theta}(x)$
transforms into expressions (30) and (31). Taking into account expressions for
real and imaginary part of coefficients $g$, $f$ and $d$ we arrive at

$f_{\theta}(x)=\frac{1}{2\sqrt{\pi
r_{f}}}e^{-\frac{i}{2}\varphi_{f}}e^{-Gx^{2}},$ (63)

where

$r_{f}=\frac{\sigma^{2}}{\mu}\left[\left\\{64\sigma^{4}+20T^{2}+\sigma^{4}T^{4}\right\\}+T^{2}\sigma^{4}\left\\{16(\sigma^{4}-4)+T^{2}(1-4\sigma^{4})\right\\}^{2}\right]^{1/2},$
$\varphi_{f}=\arctan\left(\sigma^{2}T\frac{16(\sigma^{4}-4)+T^{2}(1-4\sigma^{4})}{64\sigma^{4}+20T^{2}+T^{4}\sigma^{4}}\right),$
$G=r_{g}\sin\frac{\theta}{2}e^{i\varphi_{g}}-\frac{1}{4}\frac{1}{r_{f}}\left\\{r_{\alpha}^{2}\sin^{2}\frac{\theta}{2}e^{i(\varphi_{f}+2\varphi_{d})}-2r_{d}\sin\frac{\theta}{2}e^{i(\varphi_{f}+\varphi_{d}+\frac{\pi}{2})}-\cos^{2}\frac{\theta}{2}e^{i\varphi_{f}}\right\\}.$

In order to make the formula shorten certain symbols are introduced

$r_{g}=\frac{\sigma^{2}}{t_{4}}\left[\left\\{64+20T^{2}\sigma^{4}+T^{4}\right\\}+\frac{T^{2}}{\sigma^{4}}\left\\{16(1-4\sigma^{4})+T^{2}(\sigma^{4}-4)\right\\}^{2}\right]^{1/2},$
$\varphi_{g}=\arctan\frac{g_{1}}{g_{0}}=\arctan\left\\{\frac{T[16(1-4\sigma^{4})+T^{2}(\sigma^{4}-4)]}{\sigma^{2}(64+20T^{2}\sigma^{4}+T^{4})}\right\\},$
$r_{d}=\frac{4T}{t_{4}}\left[\sigma^{4}(16-T^{2})(1-\sigma^{4})^{2}+16T^{2}(1-\sigma^{8})^{2}\right]^{1/2},$
$\varphi_{d}=\arctan\left(\frac{16T^{2}(1-\sigma^{8})}{4T\sigma^{2}(1-\sigma^{4})(16-T^{2})}\right)=\arctan\left(4T\frac{(1-\sigma^{8})}{\sigma^{2}(1-\sigma^{4})(16-T^{2})}\right).$

Thus, suggested restoration scheme (63) in limiting cases precisely renew
Gaussian signal (25)

$f_{(t=1)}(x)=f_{(t=0)}(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp{\left(-\frac{x^{2}}{2\sigma^{2}}\right)}.$
(64)

Thereby, side by side with known limiting cases arises possibility of signal
intensity distribution restoration in the region $t=[0..1]$. In our opinion
investigation of intensity distribution in the mentioned region is an urgent
problem, however it passes the limits of the present investigation.

##  5 Conclusions

In this work we propose a mutual space-frequency distribution as
generalization of Weyl distribution (16). The mutual distributions is
characterized by a certain parameter $t$. It is a generalization of space-
frequency representations suggested by Wigner and Weyl and comprises them as
limiting cases. In the course of the investigation it has been determined that
transition from Weyl distribution into Wigner one occurs by means of mutual
region as a curve at the information diagram of mutual coordinates $(x,p)$
Fig.2. When modificating the mutual parameter $t$ the distribution changes and
simultaneously suffers deformation for the angle propotional to parameter $t$.
Fig.3 depicts the region of mutual distribution close to Weyl distribution
(quasi-Weyl region of distribution). Fig. 3(a) illustrates Weyl distribution
of Gaussian signal that is formed from the mutual distribution (51) when the
value of mutual parameter $t=0$. Increase of value of mutual parameter till
$t=0,1$ leads to the curve of the mutual distribution at informational diagram
(Fig.3(b)). When $t=0,25$ beside curve also starts the pocess of deformation
that leads to the transformation of Weyl distribution into Wigner distribution
(Fig.3(c)). The peculiar is the value of parameter $t=0,5$. mutual
distribution in this case is placed precisely in the middle between limiting
cases of Weyl and Wigner distributions (Fig.5). In the process of increase of
the mutual parameter $t$ a transformation of Weyl distribution into Wigner
distribution takes place by means of change of the curve counterclockwise to
the mutual space-frequency distribution. In the limiting case $t=1$ from the
mutual distribution Wigner distribution is formed (Fig.4(a)). When the mutual
parameter $t$ decreases in the region of Wigner distribution the curve at the
informational diagram of joined coordinates $(x,p)$ is changed (Fig.4(b,c)).
In the region of Wigner distribution when the parameter $t$ decreases the
curve of mutual distribution is changed clockwise. When the value $t=0,5$ the
mutual distribution is formed what can be observed at Fig.5. Thereby, we come
to conclusion that in the process of changing the parameter $t$ of the mutual
distribution the Weyl and Wigner distributions move towards one another at
informational diagram and are put in equilibrium in the point $t=0,5$ (Fig.5).
It can be stated that Wigner distribution is formed as a change of curve of
Weyl distribution at informational diagram for the angle $\theta<90^{0}$.
Similarly, Weyl distribution is formed as a change of curve of Wigner
distribution in the contrary way. It should be noticed that in the mutual
region the distribution becomes complex one (Fig.6(b-d)), (Fig.7(b-d)).
However, in the known liniting cases only the real part has a contribution:
$t=0$ (Fig.6(a)),(Fig.7(a)) and $t=1$ (Fig.6(e)),(Fig.7(e)).

In this paper we propose new space-frequency distribution (16), which could
play an important role in the optical information processing schemes
describing.

## References

  * [1] E.P. Wigner, Phys. Rev. 40, 749(1933).
  * [2] D. Gabor, Journal of the IEE 93, 429(1946).
  * [3] J. Ville, Cables et Transmissions 2A, 61(1948).
  * [4] T. Claasen W. Mecklenbrauker, Philips Jour. of Research 35, 217(1980).
  * [5] T. Claasen W. Mecklenbrauker, Philips Jour. of Research 35, 276(1980).
  * [6] T. Claasen W. Mecklenbrauker, Philips Jour. of Research 35, 372(1980).
  * [7] L.Cohen, 77, 941(1989).
  * [8] B. Boashash, Time-frequency Signal Analysis and Processing (Elsevier, 2003, 743 p.).
  * [9] G. Cristobal, C. Gonzalo, J. Bescos, Advances in Electronics and Electron Physics Series 80, 309(1991).
  * [10] L. Cohen Time-frequency Analysis (Prentice Hall, 1995, 299 p.).
  * [11] A. Mertings Signal Analysis (WileySons, 1999, 310 p.).
  * [12] Ch. Kurtsiefer, T. Pfau, J. Mlynek, Nature 386, 150(1997).
  * [13] G. Breitenbach, S. Schiller, J. Mlynek, Nature 387, 471(1997).
  * [14] D. Smithey, M. Beck, M. Raymer, Phys. Rev. Lett. 70(9), 1244(1996).
  * [15] L. Cohen, J. Math. Phys. 7, 781(1966).
  * [16] N. De Bruijn, Nieuw Archief voor Wiskunde (3)(XXI), 205(1973).
  * [17] L. Durak, PhD thesis, Institute of engineering and science of Bilken University, Turkey, 2003, 139 p.
  * [18] M. Shovgenyuk, Yu. Kozlovskii, Dep. NAS Ukraine, 6, 92(2000).
  * [19] S. Chountasis, A. Vourdas and C. Bendjaballah, Phys. Rev. A. 60(5), 3467(1999).
  * [20] S. Chountasis, A. Vourdas, Phys. Rev. A. 58(2), 848(1998).
  * [21] S. Chountasis, A. Vourdas, Phys. Rev. A. 58(3), 1794(1999).
  * [22] M. Shovgenyuk, Preprint ICMP-92-25U, Lviv, (1992). (www.icmp.lviv.ua).
  * [23] Yu.M. Kozlovskii, Ukr.J.Phys.Opt. 3, 124(2003).
  * [24] Yu.M. Kozlovskii, Preprint ICMP-06-27U, Lviv, (2006). (www.icmp.lviv.ua).
  * [25] Yu.M. Kozlovskii, Ukr.Phys.Journ. 3, 124(2008).

Figure 1: Real and Image parts of the Gaussian mutual distribution constant
$C_{t}^{\sigma}$.

Figure 2: Redistribution of the Gaussian mutual distribution from Weyl to
Wigner domain.

Figure 3: Rotational displaysment of the Gaussian mutual distribution in the
Weyl domain at different values of mutual parameter $t$ (Real part).

Figure 4: Rotational displaysment of the Gaussian mutual distribution in the
Wigner domain at different values of mutual parameter $t$ (Real part).

Figure 5: Mutual distribution of the Gaussian signal at the value of mutual
parameter $t=0.5$ (Real part).

Figure 6: Mutual space-frequency distribution in direction $x=0$ of Gaussian
signal at different values of mutual parameter $t$, solid line real part and
dot line image part.

Figure 7: Mutual space-frequency distribution in direction $p=0$ of Gaussian
signal at different values of mutual parameter $t$, solid line real part and
dot line image part.


# Consensus on Lie groups for
the Riemannian Center of Mass

Spencer Kraisler, Student Member, IEEE, Shahriar Talebi, Student Member, IEEE,
and Mehran Mesbahi, Fellow, IEEE The authors are with the William E. Boeing
Department of Aeronautics and Astronautics, University of Washington, Seattle,
WA, USA. Emails<EMAIL_ADDRESS>The research of the
authors has been supported by NSF grant ECCS-2149470 and AFOSR grant
FA9550-20-1-0053.

###### Abstract

In this paper, we develop a consensus algorithm for distributed computation of
the Riemannian Center of Mass (RCM) on Lie Groups. The algorithm is built upon
a distributed optimization reformulation that allows developing an intrinsic,
distributed (without relying on a consensus subroutine), and a computationally
efficient protocol for the RCM computation. The novel idea for developing this
fast distributed algorithm is to utilize a Riemannian version of distributed
gradient flow combined with a gradient tracking technique. We first guarantee
that, under certain conditions, the limit point of our algorithm is the RCM
point of interest. We then provide a proof of global convergence in the
Euclidean setting, that can be viewed as a “geometric” dynamic consensus that
converges to the average from arbitrary initial points. Finally, we proceed to
showcase the superior convergence properties of the proposed approach as
compared with other classes of consensus optimization-based algorithms for the
RCM computation.

###### Index Terms:

Consensus on Lie Groups, Riemannian Center of Mass, Karcher Mean, Multiagent
Systems

##  I Introduction

Consensus algorithms are a ubiquitous class of protocols with pertinent
applications to fields such as distributed estimation, optimization, and
control of multi-agent systems. At a foundational level, consensus algorithms
steer a set of dynamical agents towards a single point. Average consensus
necessitates this point to be the average of the agents’ initial states.
Statistical and computational advantages of average consensus algorithms have
found a wide range of applications in distributed resource allocation,
formation control, and distributed estimation [1, 2, 3].

While consensus algorithms have predominantly been studied for Euclidean
spaces, there has been a number of efforts in the literature to generalize
this protocol to Riemannian manifolds [4, 5, 6, 7, 8, 9, 10]. One notable
application is 3D localization of camera sensors [11]. In this direction,
useful notions from Euclidean geometry pertinent to the notion of average
consensus generalize to Riemannian manifolds. This includes the Riemannian
Center of Mass (RCM), replicating the notion of an average [12].

The Riemannian consensus algorithms discussed so far achieve consensus but not
consensus to the RCM (or the so-called “RCM consensus”). To the best of our
knowledge, [13, §3.2] is the earliest work on the topic of RCM consensus. The
method was later extended to a larger regime of Riemannian manifolds in [14].
In order to guarantee RCM consensus, this approach utilizes a consensus
subroutine within each iteration. However, practical implementation of this
subroutine may not be always favorable and time-efficient. This calls for a
fast distributed RCM consensus algorithm that is,  1. intrinsic, and as such,
parameterization of the manifold does not affect its properties;  2.
completely distributed without relying a consensus subroutine; and  3. has at
least a linear convergence rate.

In this paper, we first reformulate the RCM for a set of points on a
Riemannian manifold as the unique optimal solution to a consensus optimization
problem. This perspective calls to explore an array of techniques including
distributed optimization on Riemannian manifolds and dynamic average
consensus. In this direction, we also provide an optimizer for Lie groups
equipped with a bi-invariant metric with the aforementioned properties, thus
achieving a distributed RCM consensus algorithm. We prove that, under certain
conditions, the limit point of our algorithm is the RCM of the original
points. We provide global convergence guarantees in the Euclidean case and
compare the performance of our proposed algorithm with distributed constrained
optimization approaches such as penalty-based and Lagrangian-based methods.

The rest of the paper is organized as follows. §II offers the problem
formulation and a technical background on Riemannian geometry, RCM and Lie
groups. In §III, we present the RCM optimization reformulation utilized to
develop fast and simple distributed RCM consensus. We then proceed with limit
point analysis and global convergence guarantees of the proposed algorithm for
the Euclidean space in §IV. Finally, several simulation scenarios are
presented in §V, followed by concluding remarks in §VI.

##  II Problem Statement and Background

Before delving deeper into necessary mathematical background, let us first
present the problem statement. Consider a Lie group $\mathcal{G}$ equipped
with a bi-invariant Riemannian metric and a communication network of agents
$[N]=\\{1,2,\ldots,N\\}$ represented by an undirected connected graph
$G=([N],E)$, were $E$ represents the edge set. Given $N$ points
$z_{i}\in\mathcal{G}$, the goal is to design a distributed (dynamic) protocol
(with respect to $G$) for each agent as
$\dot{x}_{i}(t)=\mathbf{F}_{i}(t,x_{i},x_{j}:j\sim i)$, steering the agents’
states to the RCM of $\\{z_{1},...,z_{N}\\}$.

###  II-A Riemannian Geometry

For geometric concepts, we follow the standard notation as in [15]. Let
$(\mathcal{M},\langle.,.\rangle)$ be a Riemannian manifold and
$\mathrm{d}(.,.)$ the corresponding induced geodesic distance. In this paper,
we assume that $\mathcal{M}$ has a bounded sectional curvature. We denote its
N-fold Cartesian product as $\mathcal{M}^{N}$ and use boldface letter
${\bm{x}}=(x_{1},\cdots,x_{N})\in\mathcal{M}^{N}$ to denote a point on the
product space.

We denote the open geodesic ball centered at $x\in\mathcal{M}$ with radius
$r>0$ as,

$\mathcal{B}(x,r):=\\{y\in\mathcal{M}:\mathrm{d}(x,y)<r\\}.$

A subset $A\subset\mathcal{M}$ is called geodesically convex (g-convex) if for
any $x,y\in A$, there exists a unique minimizing geodesic
$\gamma:[0,T]\to\mathcal{M}$ connecting $x$ to $y$ such that $\gamma$ is
contained in $A$ and there exists no other minimizing geodesic connecting the
pair.111Some authors call this strong geodesic convexity. A function
$f:A\to\mathbb{R}$ is called g-convex if for any geodesic $\gamma$ contained
in $A$, $f\circ\gamma$ is convex in the Euclidean sense [16]. The convexity
radius of $\mathcal{M}$ is defined as

$r^{*}:=\frac{1}{2}\min(\mathrm{inj}(\mathcal{M}),\frac{\pi}{\sqrt{\Delta}}),$

where $\mathrm{inj}(\mathcal{M})$ is the injectivity radius and $\Delta$ is
the upper bound on the sectional curvature of $\mathcal{M}$. If $r\leq r^{*}$,
then $\mathcal{B}(x,r)$ is g-convex [12, 4].

Next, for ${\bm{z}}=(z_{1},\dots,z_{N})\in\mathcal{M}^{N}$, define its
Riemannian Center of Mass (RCM) as any minimizer of,

$\min_{y\in\mathcal{M}}\;\sum_{i=1}^{N}\mathrm{d}(y,z_{i})^{2};$ (1)

the RCM may not exist nor be unique. We denote the RCM of ${\bm{z}}$ as
$\text{RCM}({\bm{z}})$ or $\bar{z}$ when it exists and is unique. One can
induce existence and uniqueness guarantees as follows. Define the convexity
submanifold:

$\mathcal{C}:=\\{{\bm{z}}\in\mathcal{M}^{N}:\exists y\in\mathcal{M},r<r^{*}\
\text{s.t. }\ z_{i}\in\mathcal{B}(y,r)\ \forall i\\}.$ (2)

If ${\bm{z}}\in\mathcal{C}$, then $\text{RCM}({\bm{z}})$ exists and is unique
[12, Thm. 2.1]. Now, let $\mathcal{B}$ be any geodesic ball satisfying the
condition in equation 2 for ${\bm{z}}$. If some $z\in\mathcal{B}$ satisfies
the Karcher equation,

$\sum_{i=1}^{N}\log_{z}(z_{i})=0,$ (3)

then we necessarily have $z=\text{RCM}({\bm{z}})$.

###  II-B Lie Groups

A Lie group $\mathcal{G}$ is a smooth manifold with a group structure where
the group and inverse operators are smooth mappings. Every (Lie) group admits
the identity element $e$ such that $xe=ex=x$ and $xx^{-1}=x^{-1}x=e$ for all
$x\in\mathcal{G}$. The tangent space of $\mathcal{G}$ at the identity is
called the Lie algebra $\mathfrak{g}:=T_{e}\mathcal{G}$. The Lie algebra is a
vector space equipped with a vector multiplication operator called the
commutator. We suggest [17, 18] for introductions to Lie groups with emphasis
on Riemannian geometry and control theory.

A Riemannian metric $\langle.,.\rangle$ on $\mathcal{G}$ is called left-
invariant if

$\langle\xi,\eta\rangle_{e}=\langle dL_{x}\xi,dL_{x}\eta\rangle_{x},$

for any $x\in\mathcal{G}$, $\xi,\eta\in\mathfrak{g}$. Here, $dL_{x}$ is the
differential of the left-translation map $L_{x}:y\mapsto xy$. Right-invariant
metrics are defined similarly with the right-translation map $R_{x}:y\mapsto
yx$. Bi-invariant metrics are both left- and right-invariant. If
$\mathcal{G}\subset GL(n)$ is a matrix Lie group, then $dL_{x}\xi=x\xi$.

All Lie groups admit left- and right-invariant metrics [15, Lemma 3.10]. A Lie
group admits a bi-invariant metric if and only if it is isomorphic to the
Cartesian product of a compact Lie group and a vector space [19, Lemma 7.5].
An example of a Lie group with a bi-invariant metric is the unit quaternions
equipped with the dot product.

The manifold and Lie group exponential coincide at the identity when the
metric is bi-invariant [17, Prop. 3.10.]. Since left-translation is an
isometry, the following identities hold for these Lie groups [15, Prop. 5.9]:
for any $x,y\in\mathcal{G}$ and $\xi\in\mathfrak{g}$,

$\displaystyle\exp_{x}(dL_{x}\xi)$ $\displaystyle=L_{x}(\mathrm{Exp}(\xi)),$
$\displaystyle\log_{x}(y)$ $\displaystyle=dL_{x}\mathrm{Log}(x^{-1}y),$

where $\mathrm{Exp}$ and $\exp_{x}$ are, respectively, the Lie and manifold
exponential at $x$; $\mathrm{Log}$ and $\log_{x}$, on the other hand, are
their corresponding inverse mappings (wherever they are well-defined). This
also implies that on $T_{x}\mathcal{G}$, we have

$\nabla_{x}\left[\frac{1}{2}\mathrm{d}(x,y)^{2}\right]=-\log_{x}(y)=-dL_{x}\mathrm{Log}(x^{-1}y).$

##  III The Algorithm and its Derivation

In this section, we present a reformulation of equation 1 and propose a
solution algorithm that is intrinsic, distributed, and has an (empirically)
linear convergence rate.

###  III-A RCM Consensus Reformulation

Using a consensus reformulation, we can reformulate equation 1 as a g-convex
consensus optimization problem,

$\displaystyle\min_{{\bm{x}}\in\mathcal{M}^{N}}$
$\displaystyle\mathbf{f}({\bm{x}}):=\sum_{i=1}^{N}\mathrm{d}(x_{i},z_{i})^{2}$
(5a) s.t. $\displaystyle\;{\bm{x}}\in\mathcal{A},$ (5b)

where $\mathcal{A}=\\{(x,...,x):x\in\mathcal{M}\\}$ is the agreement
submanifold. Also, note that on a connected graph, equation 5b is equivalent
to zero consensus error $\varphi({\bm{x}})=0$, defined as

$\varphi({\bm{x}}):=\frac{1}{2}\sum_{\\{i,j\\}\in
E}\mathrm{d}(x_{i},x_{j})^{2}=\frac{1}{4}\sum_{i=1}^{N}\sum_{j\sim
i}\mathrm{d}(x_{i},x_{j})^{2}.$ (6)

Thus, the unique solution to equation 5 on $\mathcal{C}$ (where its existence
and uniqueness is guaranteed) corresponds to
${\bm{x}}^{*}:=(\bar{z},...,\bar{z})$.

((a))

((b))

Figure 1:  (a) Block diagram of Algorithm 1. (b) Illustration of the
correspondence in Algorithm 1 between the state consensus on the Lie group
(left) and the gradient dynamic consensus on its Lie algebra (right),
interconnected by two mappings.

This optimization reformulation enables us to design a distributed algorithm
for RCM consensus using local first order information. That is, for agent $i$
we use,222These closed-form expressions only hold for Lie groups with bi-
invariant metrics. Nonetheless, equation 5 is still distributed for arbitrary
Riemannian manifolds but equation 7 will be different.

$\displaystyle\nabla_{x_{i}}\mathbf{f}({\bm{x}})$
$\displaystyle=-dL_{x_{i}}\mathrm{Log}(x_{i}^{-1}z_{i}),$ (7a)
$\displaystyle\nabla_{x_{i}}\varphi({\bm{x}})$ $\displaystyle=\textstyle-
dL_{x_{i}}\sum_{j\sim i}\mathrm{Log}(x_{i}^{-1}x_{j}),$ (7b)

which depend only on the variables $\\{x_{i},x_{j}:j\sim i\\}$ that are
locally available. Therefore, a fast first order optimizer will be a solution
to our problem.

###  III-B Our Algorithm

######  Algorithm 1.

Given ${\bm{z}}=(z_{1},\cdots,z_{N})\in\mathcal{C}\subset\mathcal{G}^{N}$, let
$\mathbf{L}:\mathfrak{g}^{N}\to\mathfrak{g}^{N}$ be the graph Laplacian of
$G$, represented as a linear operator on the product Lie algebra
$\mathfrak{g}^{N}$. For each agent $i$, set $x_{i}(0):=z_{i}\in\mathcal{G}$
and $w_{i}(0):=0\in\mathfrak{g}$. The proposed algorithm assumes the following
dynamics:

$\displaystyle\dot{{\bm{x}}}$
$\displaystyle=-\nabla\varphi({\bm{x}})-dL_{{\bm{x}}}{\bm{v}}$ (8a)
$\displaystyle\dot{{\bm{w}}}$ $\displaystyle=\mathbf{L}{\bm{v}},$ (8b)

where ${\bm{v}}=-{\bm{w}}+dL_{{\bm{x}}^{-1}}\nabla\mathbf{f}({\bm{x}})$.

In order to show that Algorithm 1 is fully distributed, notice that the
dynamics for each agent $i$ reduces to

$\displaystyle\dot{x}_{i}$
$\displaystyle=dL_{x_{i}}\big{(}\textstyle\sum_{j\sim
i}\mathrm{Log}(x_{i}^{-1}x_{j})-v_{i}\big{)},$ $\displaystyle\dot{w}_{i}$
$\displaystyle=\textstyle\sum_{j\sim i}(v_{i}-v_{j}),$

with $v_{i}=-w_{i}+\mathrm{Log}(z_{i}^{-1}x_{i})$. Later, we show that each
$v_{i}\in\mathfrak{g}$ is tracking a global information about the distributed
cost $\mathbf{f}$. Also, $w_{i}$ is a latent state for a dynamic consensus
algorithm on the Lie algebra equation 10.

###  III-C Derivation of our Algorithm

Similar to the gradient flow as the simplest (continuous) optimizer for an
optimization problem, the Distributed Gradient Flow (DGF) can be viewed as the
simplest optimizer for a distributed optimization problem. The state dynamics
for DGF assumes the form,

$\dot{{\bm{x}}}=-\nabla\varphi({\bm{x}})-\nabla\mathbf{f}({\bm{x}}).$ (9)

Locally, each agent’s state follows the dynamics

$\textstyle\dot{x}_{i}=dL_{x_{i}}\sum_{j\sim
i}\mathrm{Log}(x_{i}^{-1}x_{j})+dL_{x_{i}}\mathrm{Log}(x_{i}^{-1}z_{i}).$

The first forcing term above drives the consensus error $\varphi$ to zero,
whereas the second term attempts to minimize the local cost
$x\mapsto\frac{1}{2}\mathrm{d}(x,z_{i})^{2}$. This is a first attempt at
minimizing the cost
$x\mapsto\mathbf{f}(x,...,x)=\sum_{i=1}^{N}\frac{1}{2}\mathrm{d}(x,z_{i})^{2}$
in a distributed manner, whose global minimizer is $x=\bar{z}$. However, this
method easily fails due to the excess of points that satisfy
$\nabla\varphi({\bm{x}})=-\nabla\mathbf{f}({\bm{x}})$ with
${\bm{x}}\not\in\mathcal{A}$.

One can think of DGF as an “open loop” distributed optimizer for Equation 5,
i.e., while Equation 9 includes global state feedback with the consensus term,
there is no global cost gradient feedback. To resolve this issue, each agent
needs access to the global information about the average gradients of local
costs of all agents. The simplest way to introduce this cost gradient feedback
is through Gradient Tracking (GT) [20]. This involves implementing a dynamic
consensus algorithm on the Lie algebra for tracking this average gradient:

$\displaystyle\dot{\mathbf{w}}$
$\displaystyle=-\mathbf{L}\mathbf{w}+\mathbf{L}\mathbf{u},$ (10a)
$\displaystyle\mathbf{v}$ $\displaystyle=-\mathbf{w}+\mathbf{u}.$ (10b)

Here, $w_{i}(0)=0$ and we set
$u_{i}=dL_{x_{i}^{-1}}\nabla_{x_{i}}\mathbf{f}({\bm{x}})$. Under these
dynamics, each $v_{i}$ will track the same global information
$\frac{1}{N}\sum_{j=1}^{N}dL_{x_{j}^{-1}}\nabla_{x_{j}}\mathbf{f}({\bm{x}})$.

Finally, combining DGF with GT, which is achieved by replacing the term
$-\nabla\mathbf{f}({\bm{x}})$ in equation 9 with $-dL_{{\bm{x}}}{\bm{v}}$,
gives us Equation 8. In other words, we are “closing the loop” on the DGF by
introducing a cost gradient feedback for each agent using GT. A schematic
diagram of this design is also depicted in Figure 1(a) with its geometric
interpretation in Figure 1(b).

######  Remark 1.

We do no require $z_{i}$’s to be constant. Suppose each $z_{i}$ is time-
varying and we wish for each $x_{i}$ to track the time-varying RCM in a
distributed manner. Indeed, we have observed that our algorithm works
empirically in this dynamic consensus setup.

######  Remark 2.

Note that we have initialized ${\bm{x}}(0)=z_{i}$. This can be relaxed to only
requiring ${\bm{x}}(0)\in\mathcal{C}$. The “closer” the points are initialized
to $\bar{z}$, the faster the convergence. However, since $\bar{z}$ is unknown
a priori, we simply chose ${\bm{x}}_{i}(0)=z_{i}$; another candidate would be
${\bm{x}}_{i}(0)=e$.

##  IV Main Results

In this section, we first provide the limit point analysis of Algorithm 1 in
the general case and then guarantee convergence for the Euclidean case.

######  Proposition 1 (Limit point).

Let $\mathcal{B}\subset\mathcal{G}$ be a geodesic ball with radius $r<r^{*}$
containing $\\{z_{1},...,z_{N}\\}$, and consider the trajectory
$({\bm{x}}(.),{\bm{w}}(.))$ generated by Algorithm 1. Suppose
$({\bm{x}}(.),{\bm{w}}(.))$ converges to some $({\bm{x}}^{*},{\bm{w}}^{*})$
such that ${\bm{x}}^{*}\in\mathcal{B}^{N}$. Then $x_{i}^{*}=\bar{z}$ and
$w_{i}^{*}=\mathrm{Log}(z_{i}^{-1}\bar{z})$ for each $i$.

###### Proof.

We can compactly express Algorithm 1 as

$\begin{bmatrix}\dot{{\bm{x}}}\\\
\dot{{\bm{w}}}\end{bmatrix}=\begin{bmatrix}-\nabla\varphi({\bm{x}})-dL_{\bm{x}}[-{\bm{w}}+\mathrm{Log}({\bm{z}}^{-1}{\bm{x}})\\\
\mathbf{L}[-{\bm{w}}+\mathrm{Log}({\bm{z}}^{-1}{\bm{x}})]\end{bmatrix}.$

Since we assume $({\bm{x}}^{*},{\bm{w}}^{*})$ is a fixed point,

$\mathbf{L}\left[-{\bm{w}}^{*}+\mathrm{Log}({\bm{z}}^{-1}{\bm{x}}^{*})\right]=0,$

and hence,

$-{\bm{w}}^{*}+\mathrm{Log}({\bm{z}}^{-1}{\bm{x}}^{*})=(\xi,...,\xi)=:\xi{\bm{1}}\in\mathfrak{g}^{N},$

for some $\xi\in\mathfrak{g}$. Hence,
$dL_{({\bm{x}}^{*})^{-1}}\nabla\varphi({\bm{x}}^{*})=\xi{\bm{1}}$. Then

$N\|\xi\|^{2}=\langle\xi{\bm{1}},\xi{\bm{1}}\rangle=\langle
dL_{({\bm{x}}^{*})^{-1}}\nabla\varphi({\bm{x}}^{*}),\xi{\bm{1}}\rangle.$

Since the metric is bi-invariant, we have

$\displaystyle
dL_{(x_{i}^{*})^{-1}}\mathrm{Log}_{x_{i}^{*}}(x_{j}^{*})=\mathrm{Log}((x_{i}^{*})^{-1}x_{j}^{*})=-dL_{(x_{j}^{*})^{-1}}\mathrm{Log}_{x_{j}^{*}}(x_{i}^{*})$

implying $\textstyle\sum_{i=1}^{N}dL_{(x_{i}^{*})^{-1}}\sum_{j\sim
i}\log_{x_{i}^{*}}(x_{j}^{*})=0$. Therefore

$\displaystyle\langle
dL_{({\bm{x}}^{*})^{-1}}\nabla\varphi({\bm{x}}^{*}),\xi{\bm{1}}\rangle=\left\langle\sum_{i=1}^{N}\sum_{j\sim
i}\mathrm{Log}((x_{i}^{*})^{-1}x_{j}^{*}),\xi\right\rangle=0.$

Thus we must have $\xi=0$, and therefore $\nabla\varphi({\bm{x}}^{*})=0$.

Next, it is shown in [4, Theorem 5] that the critical points of $\varphi$
restricted to $\mathcal{C}$ coincide with $\mathcal{A}$. So,
${\bm{x}}^{*}\in\mathcal{C}$ and $\nabla\varphi({\bm{x}}^{*})=0$ together
imply ${\bm{x}}^{*}\in\mathcal{A}$. Next, note that

$\textstyle\frac{d}{dt}\sum_{i=1}^{N}\dot{w}_{i}(t)=0.$

As such, since ${\bm{w}}(0)=0$, we must have $\sum_{i=1}^{N}w_{i}^{*}=0$.
Therefore, by the fact that
$-{\bm{w}}^{*}+\mathrm{Log}({\bm{z}}^{-1}{\bm{x}}^{*})=0$, we arrive at

$\textstyle\sum_{i=1}^{N}\mathrm{Log}(z_{i}^{-1}x^{*})=0.$

Since $x^{*},z_{1},...,z_{N}\in\mathcal{B}$ and $x^{*}$ satisfies Equation 3,
$x^{*}=\bar{z}$, and thus $x_{i}^{*}=\bar{z}$ and
$w_{i}^{*}=\mathrm{Log}(z_{i}^{-1}\bar{z})$. ∎

Next, we provide the convergence guarantees of our algorithm for
$\mathcal{G}=\mathbb{R}^{n}$; the corresponding analysis for arbitrary Lie
groups is the subject of our future work. Herein, we write
${\bm{x}}=\mathrm{col}(x_{1},...,x_{N})\in\mathbb{R}^{Nn}$ to denote vertical
concatenation and use similar notation for ${\bm{z}}$, ${\bm{x}}$, and
${\bm{u}}$. Then the dynamics of Algorithm 1 reduces to the following matrix
form:

$\displaystyle\dot{{\bm{x}}}$ $\displaystyle=-(\mathbf{L}\otimes
I_{n}){\bm{x}}+{\bm{w}}-({\bm{x}}-{\bm{z}})$ (11a)
$\displaystyle\dot{{\bm{w}}}$ $\displaystyle=-(\mathbf{L}\otimes
I_{n}){\bm{w}}+(\mathbf{L}\otimes I_{n})({\bm{x}}-{\bm{z}}),$ (11b)

where $\otimes$ denote the Kronecker product. The RCM of initial states
${\bm{z}}$ reduces to the Euclidean average

$\textstyle\bar{z}=\text{RCM}({\bm{z}})=\frac{1}{N}\sum_{i=1}^{N}z_{i}\in\mathbb{R}^{n},$

and so
${\bm{x}}^{*}=\mathrm{col}(\bar{z},...,\bar{z})={\bm{1}}_{N}\otimes\bar{z}$.
The next result establishes guaranteed average consensus starting from
arbitrary initial points.

######  Theorem 1 (Convergence in $\mathbb{R}^{n}$).

Suppose ${\bm{x}}(t),{\bm{w}}(t)$ is the trajectory generated by equation 11
over a connected graph $G$ with ${\bm{w}}(0)=0$. Then, starting from any
arbitrary starting point ${\bm{x}}(0)$, we have an exponentially convergent
trajectory with $\lim_{t\to\infty}x_{i}(t)=\bar{z}$ and
$\lim_{t\to\infty}w_{i}(t)=\bar{z}-z_{i}$ for all $i$.

###### Proof.

First, by noting that ${\bm{v}}=-{\bm{w}}+{\bm{x}}-{\bm{z}}$ we reformulate
the system in equation 11 as

$\begin{bmatrix}\dot{{\bm{v}}}\\\ \dot{{\bm{x}}}\end{bmatrix}=(A\otimes
I_{n})\begin{bmatrix}{\bm{v}}\\\ {\bm{x}}\end{bmatrix},$ (12)

where

$A\coloneqq\begin{bmatrix}-\mathbf{L}-I_{N}&-\mathbf{L}\\\
-I_{N}&-\mathbf{L}\end{bmatrix}.$ (13)

Before proceeding, we need the following results on characterizing the
spectrum of $A$ with its proof deferred to the end of this section.

######  Lemma 1.

Suppose $G$ is connected. Then zero is a simple eigenvalue of the matrix $A$
in (13). Furthermore, all non-zero eigenvalues of this matrix have negative
real parts.

The spectral properties of $A$ established above then implies that the
dynamics of equation 13 is marginally stable. Let $p,q$ be, respectively, the
right and left eigenvectors of $A$ associated with $\lambda=0$ such that
$p^{\top}q=1$. So, $p\in\mathcal{N}(A)$. Also, since 0 is a simple eigenvalue
of $A^{\top}$, we get

$\mathcal{N}(A^{\top})=\mathrm{span}\left(\begin{bmatrix}{\bm{1}}_{N}\\\
-{\bm{1}}_{N}\end{bmatrix}\right).$

Thus, we set $p=\begin{bmatrix}0_{N}\\\ {\bm{1}}_{N}\end{bmatrix}$,
$q=\frac{1}{N}\begin{bmatrix}-{\bm{1}}_{N}\\\ {\bm{1}}_{N}\end{bmatrix}$. By
properties of Kronecker product, we obtain that $A\otimes I_{n}$ is also
marginally stable with a zero eigenvalue of algebraic and geometric
multiplicity $n$. This eigenvalue has corresponding right eigenvectors
$p\otimes{\bm{e}}_{i}$ and left eigenvectors $q\otimes{\bm{e}}_{i}$. Here,
${\bm{e}}_{i}\in\mathbb{R}^{n}$ is the $i$-th standard basis for
$i=1,2,\cdots,n$.

Next, by computing $e^{(A\otimes I_{n})t}$ using the Jordan decomposition and
taking the limit as $t\to\infty$ we obtain that (see [21, Proposition 3.11]
for a similar computation)

$\lim_{t\to\infty}e^{(A\otimes
I_{n})t}=\sum_{i=1}^{n}(p\otimes{\bm{e}}_{i})(q\otimes{\bm{e}}_{i})^{\top}=(pq^{\top})\otimes
I_{n}.$

Therefore, Equation 13 converges as follows

$\begin{bmatrix}{\bm{v}}^{*}\\\
{\bm{x}}^{*}\end{bmatrix}\coloneqq\lim_{t\to\infty}\begin{bmatrix}{\bm{v}}(t)\\\
{\bm{x}}(t)\end{bmatrix}=((pq^{\top})\otimes
I_{n})\begin{bmatrix}{\bm{v}}(0)\\\ {\bm{x}}(0)\end{bmatrix},$

where ${\bm{w}}(0)={\bm{0}}_{Nn}$ and ${\bm{x}}(0)$ is the arbitrary starting
point. Thus,
${\bm{v}}(0)=-{\bm{w}}(0)+{\bm{x}}(0)-{\bm{z}}={\bm{x}}(0)-{\bm{z}}$, and

$\displaystyle\begin{bmatrix}{\bm{v}}^{*}\\\ {\bm{x}}^{*}\end{bmatrix}$
$\displaystyle=\frac{1}{N}\left(\begin{bmatrix}{\bm{0}}_{N\times
N}\quad{\bm{0}}_{N\times N}\\\
-{\bm{1}}_{N}{\bm{1}}_{N}^{\top}\quad{\bm{1}}_{N}{\bm{1}}_{N}^{\top}\end{bmatrix}\otimes
I_{n}\right)\begin{bmatrix}{\bm{x}}(0)-{\bm{z}}\\\ {\bm{x}}(0)\end{bmatrix}$
$\displaystyle=\begin{bmatrix}{\bm{0}}_{Nn}\\\
\frac{1}{N}({\bm{1}}_{N}{\bm{1}}_{N}^{\top}\otimes
I_{n}){\bm{z}}\end{bmatrix}=\begin{bmatrix}{\bm{0}}_{Nn}\\\
{\bm{1}}_{N}\otimes\bar{z}\end{bmatrix}.$

Therefore ${\bm{x}}^{*}={\bm{1}}_{N}\otimes\bar{z}$ which completes the proof.
∎

###### Proof of Lemma 1.

Let $\begin{bmatrix}{\bm{v}}\\\ {\bm{x}}\end{bmatrix}\in\mathcal{N}(A)$ be
non-zero. Then

$\displaystyle-\mathbf{L}{\bm{v}}-{\bm{v}}-\mathbf{L}{\bm{x}}$
$\displaystyle=0,$ (15a) $\displaystyle-{\bm{v}}-\mathbf{L}{\bm{x}}$
$\displaystyle=0.$ (15b)

By combining the two, we get $\mathbf{L}{\bm{v}}=0$ implying that
${\bm{v}}\in\mathrm{span}({\bm{1}}_{N})$. It follows from equation 15b that
$\mathbf{L}{\bm{x}}=-{\bm{v}}=\alpha{\bm{1}}_{N}$. Then
${\bm{1}}^{\top}\mathbf{L}{\bm{x}}=0=\alpha{\bm{1}}_{N}^{\top}{\bm{1}}_{N}$,
and thus $\alpha=0$. But this implies that ${\bm{v}}=0$. Thus
$\mathbf{L}{\bm{x}}=0$, which implies that
${\bm{x}}\in\mathrm{span}({\bm{1}}_{N})$. Therefore, we can conclude that

$\mathcal{N}(A)=\left\\{\begin{bmatrix}{\bm{0}}_{N}\\\
{\bm{x}}\end{bmatrix}:{\bm{x}}\in\mathrm{span}({\bm{1}}_{N})\right\\}.$

Since the nullspace has dimension 1, it follows that 0 is a simple eigenvalue
of $A$. Next, let $\begin{bmatrix}{\bm{v}}\\\ {\bm{x}}\end{bmatrix}$ be an
eigenvector of $A$ associated with a non-zero eigenvalue $\lambda$. For the
sake of contradiction, suppose $\Re[\lambda]\geq 0$. Then

$\displaystyle\begin{cases}\lambda{\bm{v}}&=-\mathbf{L}{\bm{v}}-{\bm{v}}-\mathbf{L}{\bm{x}}\\\
\lambda{\bm{x}}&=-{\bm{v}}-\mathbf{L}{\bm{x}}\end{cases}$

Note, if ${\bm{x}}=0$, then ${\bm{v}}=0$. So, in order for
$\begin{bmatrix}{\bm{v}}\\\ {\bm{x}}\end{bmatrix}$ to be an eigenvector, we
need ${\bm{x}}\not=0$. By solving for ${\bm{x}}$, we obtain

$-\lambda{\bm{x}}=(\lambda I+\mathbf{L})^{2}{\bm{x}}.$ (16)

So, $(-\lambda,{\bm{x}})$ is an eigenpair of the matrix $(\lambda
I+\mathbf{L})^{2}$. The eigenvalues of $(\lambda I+\mathbf{L})^{2}$ are
$\\{(\lambda+\rho)^{2}:\rho\in\sigma(\mathbf{L})\\}$. So, there exists some
$\rho\in\sigma(\mathbf{L})$ such that $-\lambda=(\lambda+\rho)^{2}$. Let
$\lambda=a+bi$ and note the hypothesis of contradiction is $a\geq 0$ and
$a^{2}+b^{2}\neq 0$. By separating the real and imaginary parts, we obtain

$\displaystyle-a$ $\displaystyle=(a+\rho)^{2}-b^{2},$ $\displaystyle-b$
$\displaystyle=2(a+\rho)b.$

Suppose $b=0$, then the first equation and the hypothesis $a\geq 0$ implies
that $a=0$. But, this implies $\lambda=0$ which is a contradiction. Next,
suppose $b\not=0$, then the second equation implies $-1=2a+2\rho$. Thus
$a=\frac{-1-2\rho}{2}<0$, which is a contradiction with $a\geq 0$ as $\rho\geq
0$. Therefore, if $\lambda\neq 0$ then we must have $a<0$, and so all non-zero
eigenvalues of $A$ have negative real part. This completes the proof. ∎

##  V Simulations

In this section and for illustration purposes, we consider the special
orthogonal Lie group $\mathcal{G}=SO(3)$, with its parameterization of
$3\times 3$ rotation matrices. The associated Lie algebra, denoted
$\mathfrak{so}(3)$, is the space of $3\times 3$ skew-symmetric matrices. We
equip $SO(3)$ with the following bi-invariant Riemannian metric:

$\langle
u,v\rangle_{x}:=\frac{1}{2}\mathrm{tr}(u^{\top}v)=\frac{1}{2}\mathrm{tr}(\xi^{\top}\eta),$

where $x\in SO(3)$, $u,v\in T_{x}SO(3)$, and $\xi,\eta\in\mathfrak{so}(3)$
with $\xi=x^{-1}u,\eta=x^{-1}v$. Then, the corresponding geodesic distance
reduces to

$d(x,y)^{2}=\|\mathrm{Log}(x^{\top}y)\|^{2}=-\frac{1}{2}\mathrm{tr}\left[\mathrm{Log}(x^{\top}y)^{2}\right].$

In the following, two scenarios are in order where we randomly initialize
agents on $SO(3)$.

###  V-A Scenario 1 (Comparison of Consensus Algorithms)

In this section, we run simulations comparing our solver to three other
consensus algorithms. We randomly initialized $N=10$ agents
${\bm{x}}\in\mathcal{C}$ and run each algorithm with the same initial
conditions in order to compare them. We choose the following two metrics for
comparison: the consensus error equation 6, and the RCM Error $\textstyle
E_{\mathrm{RCM}}({\bm{x}})=\sum_{i=1}^{N}\mathrm{d}(x_{i},\bar{z})^{2}$ which
aggregates the error of each point $x_{1},\ldots,x_{N}$ from the RCM
$\bar{z}$. Using these two metrics, we compare our solver against three other
algorithms described in the following and illustrate the results in Figure 2.
For the technical details of these algorithm and their implementation see the
extended version of this work [Kraisler???] and the code [22].

####  V-A1 Algorithm 1

To implement the continuous dynamics of our algorithm, we consider their
forward-Euler discretization on each tangent space with step size
$\epsilon=0.1$ and use the Lie group exponential mapping as our choice of
retraction.

####  V-A2 Riemannian Consensus Algorithm

As stated before, [4] was one of the first works to introduce a consensus
algorithm for arbitrary Riemannian manifolds, referred to as “R. Tron et. al”
in Figure 2.

####  V-A3 Penalty method

The penalty approach is a commonly used solver for constrained optimization
problems. We implemented [23, Algorithm 14.3.1], while introducing a
distributed implementation for it. In particular, we chose penalty parameter
$\mu(s)=\frac{1}{\sqrt{s}}$, number of gradient descent iterations $S=50$, and
step size $\epsilon=0.1$.

####  V-A4 Lagrangian method

A Lagrangian method is an approach to solving constrained optimization
problems that involves finding saddle points of the Lagrangian function. We
implemented the solver described in [24, 4.4.1].

It can be seen in Figure 2 that both error metrics for our algorithm are
vanishing with linear rate. Also, note in Figure 2 that Riemannian Consensus
Algorithm has a seemingly linear rate of decrease for the consensus error,
similar to our algorithm but with faster rate. However, we emphasize that this
algorithm does not converge to the RCM and it only synchronizes the agents to
a point. This is evident in Figure 2, where the RCM Error for this algorithm
stops decreasing. Finally, as illustrated in Figure 2, the Penalty and
Lagrangian method have a sub-linear rate of decrease for both the consensus
error and the RCM error.

Figure 2: The consensus error (left) and the RCM Error (right) at each
iteration, comparing our algorithm (in blue) with three other algorithms.

###  V-B Scenario 2 (Linear rate of convergence)

In this scenario, we generate $100$ different problem instances with $N=10$
agents randomly initialized on $\mathcal{C}\subset SO(3)^{10}$. We run the
algorithm on each instance for $200$ iterations for visualization of the
converge rate. We illustrate the result in Figure 3 showing the statistics of
the RCM Error at each time step, confirming a linear rate. We note that
although convergences analysis of the algorithm has only been provided for the
Euclidean case, these numerical experiments–using different problem
parameters–demonstrate the effectiveness of the proposed approach for
distributed RCM consensus on Lie groups.

Figure 3: The RCM error of our proposed algorithm over $100$ randomly sampled
problem instances.

##  VI Conclusions

In this work, we reformulate the Riemannian Center of Mass (RCM) as a
distributed optimization problem and propose a novel and fast distributed
solver on bi-invariant Lie groups. This in turn provides the first
(completely) distributed solver for the RCM consensus problem with no
subroutines. The key idea behind the proposed algorithm is a correspondence
between a consensus protocol on the Lie group and a dynamic consensus protocol
on its Lie algebra. We have shown the properties of limit points of our
algorithm in the general setting, and guaranteed its convergence in the
Euclidean setting; the convergence analysis for the general case is the
subject of our ongoing work.

## References

  * [1] R. Olfati-Saber, J. A. Fax, and R. M. Murray, “Consensus and cooperation in networked multi-agent systems,” Proceedings of the IEEE, vol. 95, no. 1, pp. 215–233, 2007.
  * [2] S. H. Strogatz, Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Life.  Hachette UK, 2012.
  * [3] S. S. Kia, B. Van Scoy, J. Cortes, R. A. Freeman, K. M. Lynch, and S. Martinez, “Tutorial on dynamic average consensus: The problem, its applications, and the algorithms,” IEEE Control Systems Magazine, vol. 39, no. 3, pp. 40–72, 2019.
  * [4] R. Tron, B. Afsari, and R. Vidal, “Riemannian consensus for manifolds with bounded curvature,” IEEE Transactions on Automatic Control, vol. 58, no. 4, pp. 921–934, 2012.
  * [5] R. Olfati-Saber, “Swarms on sphere: A programmable swarm with synchronous behaviors like oscillator networks,” in Proceedings of the 45th IEEE Conference on Decision and Control, pp. 5060–5066, IEEE, 2006.
  * [6] A. Sarlette, S. Bonnabel, and R. Sepulchre, “Coordinated motion design on Lie groups,” IEEE Transactions on Automatic Control, vol. 55, no. 5, pp. 1047–1058, 2010.
  * [7] R. Sepulchre, “Consensus on nonlinear spaces,” IFAC Proceedings Volumes, vol. 43, no. 14, pp. 1029–1039, 2010.
  * [8] S. Chen, L. Zhao, W. Zhang, and P. Shi, “Consensus on compact Riemannian manifolds,” Inf. Sci., vol. 268, pp. 220–230, June 2014.
  * [9] S. Chen, P. Shi, W. Zhang, and L. Zhao, “Finite-time consensus on strongly convex balls of Riemannian manifolds with switching directed communication topologies,” J. Math. Anal. Appl., vol. 409, pp. 663–675, Jan. 2014.
  * [10] S. Kraisler, S. Talebi, and M. Mesbahi, “Distributed consensus on manifolds using the Riemannian center of mass,” IEEE Conference on Control Technology and Applications, 2023.
  * [11] R. Tron and R. Vidal, “Distributed 3-d localization of camera sensor networks from 2-d image measurements,” IEEE Transactions on Automatic Control, vol. 59, no. 12, pp. 3325–3340, 2014.
  * [12] B. Afsari, “Riemannian ${L}^{p}$ center of mass: Existence, uniqueness, and convexity,” Proceedings of the American Mathematical Society, vol. 139, no. 2, pp. 655–673, 2011.
  * [13] R. Tron, R. Vidal, and A. Terzis, “Distributed pose averaging in camera networks via consensus on SE(3),” in 2008 Second ACM/IEEE International Conference on Distributed Smart Cameras, pp. 1–10, IEEE, 2008\.
  * [14] R. Tron, B. Afsari, and R. Vidal, “Average consensus on riemannian manifolds with bounded curvature,” in 2011 50th IEEE Conference on Decision and Control and European Control Conference, pp. 7855–7862, IEEE, 2011.
  * [15] J. M. Lee, Introduction to Riemannian Manifolds, vol. 176.  Springer, 2nd ed., 2018.
  * [16] N. Boumal, An Introduction to Optimization on Smooth Manifolds.  Cambridge University Press, 2023.
  * [17] A. Arvanitogeōrgos, An Introduction to Lie Groups and the Geometry of Homogeneous Spaces, vol. 22.  American Mathematical Soc., 2003.
  * [18] Y. L. Sachkov, “Control theory on Lie groups,” Journal of Mathematical Sciences, vol. 156, no. 3, pp. 381–439, 2009.
  * [19] J. Milnor, “Curvatures of left invariant metrics on Lie groups,” Advances in mathematics, vol. 21, no. 3, pp. 293–329, 1976.
  * [20] A. Nedic, A. Olshevsky, and W. Shi, “Achieving geometric convergence for distributed optimization over time-varying graphs,” SIAM Journal on Optimization, vol. 27, no. 4, pp. 2597–2633, 2017.
  * [21] M. Mesbahi and M. Egerstedt, Graph Theoretic Methods in Multiagent Networks.  Princeton University Press, 2010.
  * [22] S. Kraisler, S. Talebi, and M. Mesbahi, “RCM-Consensus-Lie-Group,” GitHub repository, 2023.  Available online at https://github.com/spencerkraisler/RCM-Consensus-Lie-Group.
  * [23] A. R. Conn, N. I. Gould, and P. L. Toint, Trust Region Methods.  SIAM, 2000.
  * [24] D. P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods.  Academic press, 2014.

  *[RCM]: Riemannian Center of Mass
  *[g-convex]: geodesically convex