mydate [],

# All finite-mass Dirac monopoles

Filip Blaschke<EMAIL_ADDRESS>Research Centre for Theoretical
Physics and Astrophysics, Institute of Physics,
Silesian University in Opava, Bezručovo náměstí 1150/13, 746 01 Opava, Czech
Republic  Institute of Experimental and Applied Physics,
Czech Technical University in Prague, Husova 240/5, 110 00 Prague 1, Czech
Republic  Petr Beneš<EMAIL_ADDRESS>Institute of Experimental and
Applied Physics,
Czech Technical University in Prague, Husova 240/5, 110 00 Prague 1, Czech
Republic

###### Abstract

We present a “primitive” way of realizing finite-mass Dirac monopoles in
$U(1)$ gauge theories involving a single non-minimally interacting scalar
field. Typically, the energy density of this type of monopole is not
concentrated at its core, but it is distributed in a spherical shell, as we
illustrate on several exact solutions in the BPS limit. We show that our
construction can be interpreted as a limit of infinitely massive $W$ bosons
coupled to electromagnetic field-strength via a dipole moment. Combining our
approach with ideas of Weinberg and Lee, we present a general landscape of
$U(1)$ gauge models that support a finite-mass Dirac monopole. In fact, all
classical monopoles, i.e., Wu–Yang, ’t Hooft–Polyakov, Cho–Maison, etc., are
special points on this landscape.

Magnetic monopole; exact solutions; BPS limit

##  I Bare monopoles & Dressed monopoles

In the classical Maxwell theory ($e$ denotes the electric charge)

$\mathcal{L}=-\frac{1}{2e^{2}}F_{\mu\nu}F^{\mu\nu}\,,$ (1)

a magnetic monopole is a singular configuration of $U(1)$ gauge fields, i.e.,
the solution of $\vec{\nabla}\times\vec{A}=q\,\vec{r}/r^{3}$, with $q$ being a
dimensionless constant. As is well known, the vector potential can be
described everywhere except on a line going from the monopole to infinity, the
so-called Dirac string Dirac:1931kp . For instance, a static monopole at the
origin with a Dirac string lying on negative $z$-axis is given as

$A_{i}^{\rm D}=q\frac{\varepsilon_{i3j}x_{j}}{r(r+z)}\,.$ (2)

The shape of the string can be changed via gauge transformation.

The above description, which we dub a _bare monopole_ , is singular in three
logically distinct ways. Namely, 1) there is an unphysical singularity along
the Dirac string, 2) there is a physical singularity at the origin in both
gauge fields and in the magnetic field $\vec{B}=q\,\vec{r}/r^{3}$ (we define
$B_{i}\equiv\frac{1}{2}\varepsilon_{ijk}F_{jk}$) and 3) the energy density
$\mathcal{E}=q^{2}/(2e^{2}r^{4})$ diverges at the origin in such a way that
the total energy (i.e., the classical mass of the monopole) is infinite.

These singularities are invariably clues that point to the incompleteness of
the theoretical model. However, in so far as physics is concerned, they are
not equivalently serious. The presence of Dirac string is simply a failure of
global description of the $U(1)$ fibre bundle. The singularity in $\vec{B}$ is
a result of taking the bird’s eye point of view and describing the monopole as
a point particle. By itself, this is not a problem as we can reasonably expect
that in a more fundamental theory this singularity will be smoothed out by
monopole’s microscopic degrees of freedom. However, the real red flag is the
absence of finite energy solutions.

All in all, we are compelled to depart from the pure Maxwell theory.
Fortunately, unlike the case of _electric_ monopoles, we do not have to
abandon classical field theory to establish the existence of _finite-mass
Dirac monopoles_.

Let us now follow the (somewhat chronological) path that leads to finite-mass
field-theoretical descriptions. First, there was an observation of Wu and Yang
Wu:1976ge  that one can get rid of Dirac string by embedding (2) into an
$SU(2)$ gauge field ($q=1$):

$A_{i}^{\rm WY}=\frac{\varepsilon_{iaj}\sigma_{a}x_{j}}{2r^{2}}=A_{i}^{\rm
D}U\sigma_{3}U^{\dagger}-\mathrm{i}\,U\partial_{i}U^{\dagger}\,,$ (3)

where $\sigma_{a}$ are the Pauli matrices and where

$U=\left(\begin{array}[]{ll}\cos\big{(}\theta/2\big{)}&-\sin\big{(}\theta/2\big{)}\,\mathrm{e}^{-\mathrm{i}\varphi}\\\\[2.0pt]
\sin\big{(}\theta/2\big{)}\,\mathrm{e}^{\mathrm{i}\varphi}&\phantom{-}\cos\big{(}\theta/2\big{)}\end{array}\right)$
(4)

is a (singular) $SU(2)$ gauge transformation that relates the Wu–Yang monopole
$A_{i}^{\rm WY}$ and the Dirac monopole $A_{i}^{\rm D}$. By itself, however,
this enhancement of symmetry is not sufficient to get rid of other
singularities, neither in the $SU(2)$ “magnetic” fields, nor in the energy
density.

But, as was realized independently by ’t Hooft tHooft:1974kcl  and Polyakov
Polyakov:1974ek , the Wu–Yang monopole can be made completely regular via
spontaneous symmetry breaking instigated by adjoint scalars. Loosely speaking,
these additional fields condense at monopole’s core and smooth out the $r=0$
singularity – the monopole becomes _dressed_. Due to this field-dressing, the
’t Hooft–Polyakov monopole has none of the three singularities and represents
the best classical description of a magnetic monopole.

However, as E. Weinberg and K. Lee pointed out Lee:1994sk  – and what we want
to propagate in this paper as well – it is also possible to come up with
finite-mass monopoles without resorting to non-Abelian gauge fields (2).

The key idea of Lee:1994sk  is to couple an $U(1)$ gauge field to a complex
vector field $W_{\mu}$ with non-zero dipole moment tensor
$d_{\mu\nu}=\mathrm{i}(W_{\mu}^{\dagger}W_{\nu}-W_{\mu}W_{\nu}^{\dagger})$
that screens the bare monopole’s charge via dipole moment interactions. In
other words:

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle-\frac{1}{4g^{2}}F_{\mu\nu}^{2}-\frac{\eta}{2}d_{\mu\nu}F^{\mu\nu}-\frac{\chi
g^{2}}{4}d_{\mu\nu}^{2}$ (5)
$\displaystyle{}-\frac{1}{2}\big{|}D_{\mu}W_{\nu}-D_{\nu}W_{\mu}\big{|}^{2}+m^{2}\big{|}W_{\mu}\big{|}^{2}\,,$

where $\eta$ measures the strength of the dipole interactions with external
el. mag. field, while $\chi$ measures the dipole-dipole self-interactions. For
a static configuration, the energy density can be rewritten as

$\displaystyle\mathcal{E}$ $\displaystyle=$
$\displaystyle\frac{\chi}{4}\bigg{(}g\,d_{ij}+\frac{\eta}{\chi\,g}F_{ij}\bigg{)}^{2}+\frac{1}{4g^{2}}\bigg{(}1-\frac{\eta^{2}}{\chi}\bigg{)}F_{ij}^{2}$
(6)
$\displaystyle{}+\frac{1}{2}\big{|}D_{i}W_{j}-D_{j}W_{i}\big{|}^{2}+m^{2}\big{|}W_{i}\big{|}^{2}\,.$

In order to cancel the singularity in $F_{ij}$ we need $d_{ij}$ to have the
same type of singularity and – as we see – we must require $\eta^{2}=\chi$.
However, by itself this is not enough to make $\mathcal{E}$ regular Lee:1994sk
. Critically, the mass term would still give us $\sim 1/r^{2}$. Furthermore,
simply setting $m=0$ from the beginning would not help either. This can be
seen in the limit $\eta=\lambda=1$, where the above model can be rewritten as
pure $SU(2)$ model upon identifications $A_{\mu}=A_{\mu}^{3}$ and
$W_{\mu}=(A_{\mu}^{1}+\mathrm{i}A_{\mu}^{2})/(\sqrt{2}g)$. The static solution
is gauge-equivalent to Wu–Yang monopole which we know has divergent energy.

Thus a further field must be added, namely a real scalar field $\phi$ that
develops a non-zero expectation value at the vacuum, but which vanishes at the
monopole’s core. The idea is to make mass of $W_{\mu}$ dependent on $\phi$ in
such a way that $m^{2}(\phi)\to m^{2}(0)=0$ as $r\to 0$.

All in all, the Weinberg–Lee’s procedure for regularising magnetic monopole
can be encompassed by the following Lagrangian:

$\displaystyle\mathcal{L}_{\rm WL}$ $\displaystyle=$
$\displaystyle-\frac{1}{4g^{2}}F_{\mu\nu}^{2}-\frac{\eta}{2}d_{\mu\nu}F^{\mu\nu}-\frac{\chi
g^{2}}{4}d_{\mu\nu}^{2}$ (7)
$\displaystyle{}-\frac{1}{2}\big{|}D_{\mu}W_{\nu}-D_{\nu}W_{\mu}\big{|}^{2}+m^{2}(\phi)\big{|}W_{\mu}\big{|}^{2}$
$\displaystyle{}+\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{\lambda}{2}\big{(}\phi^{2}-v^{2}\big{)}^{2}\,.$

The total energy of a monopole solution is finite and regular if we take
$\eta^{2}=\chi$ and make an appropriate choice of $m^{2}(\phi)$. In Lee:1994sk
, these solutions were dubbed _non-topological magnetic monopoles_ , as the
mechanism that renders monopoles in these models to have finite classical
masses has nothing to do with topological reasoning á la ’t Hooft–Polyakov.
Nevertheless, the ’t Hooft–Polyakov monopole is, in fact, a special point in
the above family of models, namely $\eta=\chi=1$ and
$m^{2}(\phi)=g^{2}\phi^{2}$. Indeed, with these assignments (7) becomes a
unitary gauge of $SU(2)$ theory with an adjoint scalar
$\Phi=\frac{1}{2}\sigma_{3}\phi$.

To summarise, Eq. (7) gives us finite-mass magnetic monopoles that, compared
to bare monopoles of Eq. (1), are dressed in complex vector fields $W_{\mu}$
with explicit dipole moment interactions and a scalar field $\phi$ that
controls the mass of $W_{\mu}$’s, making it effectively zero at the monopole’s
core. This idea can be depicted as in Fig. 1.

Figure 1: Transition from a bare monopole of (1) to a “dressed” monopole of
(7).

In the next section, we introduce a “primitive” way of regularizing magnetic
monopole’s energy that is seemingly unrelated to these arguments. Our
construction will lead us to novel solutions that can be described exactly in
the BPS limit $\lambda\to 0$. However, in section III we will unify our and
Weinberg–Lee’s points of view into a much broader landscape of theories that
support finite-mass Dirac monopoles.

##  II A primitive way

The core idea of this paper is to have a scalar field $\phi$ interact directly
with the el. mag. field via non-canonical kinetic term:

$\mathcal{L}=-\frac{h^{\prime\,2}(\phi/v)}{4g^{2}}F_{\mu\nu}^{2}+\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{\lambda}{2}\big{(}\phi^{2}-v^{2}\big{)}^{2}\,.$
(8)

From the outset, there seems to be no need for complex vector field $W_{\mu}$
and dipole moment interactions (but see Sec. III). Indeed, this model can
yield finite-energy magnetic monopole solutions for appropriate choices of so-
far arbitrary function $h^{\prime}(\phi/v)$ (this notation with the prime
denoting a derivative becomes obvious when we discuss the BPS limit). Without
loss of generality we take the normalization $h^{\prime}(1)=1$.

We can depict this primitive way of a field-dressing by means of just a scalar
field as in Fig. 2.

Figure 2: Transition from a bare monopole of (1) to a “dressed” monopole of
(8).

Let us analyze the full equations of motion for a static configuration,

$\displaystyle\vec{\nabla}\times\Big{(}h^{\prime\,2}(\phi/v)\,\vec{B}\Big{)}=0\,,$
(9a)
$\displaystyle\vec{\nabla}^{2}\phi=\frac{h^{\prime}(\phi/v)\,h^{\prime\prime}(\phi/v)}{g^{2}v}\vec{B}^{2}+2\lambda\phi\big{(}\phi^{2}-v^{2}\big{)}\,.$
(9b)

We see that in the monopole background, $\vec{B}=q\,\vec{r}/r^{3}$, the first
equation is solved identically if the scalar field depend only on the
(dimensionless) radius, i.e., $\phi\equiv v\sigma(\rho)$, where $\rho=gvr$.
The second equation (subject to the boundary conditions $\sigma(0)=0$ and
$\sigma(\infty)=1$) then reads

$\frac{1}{\rho^{2}}\partial_{\rho}\big{(}\rho^{2}\partial_{\rho}\sigma\big{)}=q^{2}\frac{h^{\prime}(\sigma)\,h^{\prime\prime}(\sigma)}{\rho^{4}}+\frac{2\lambda}{g^{2}}\sigma\big{(}\sigma^{2}-1\big{)}\,.$
(10)

From the expression for the energy density

$\frac{\mathcal{E}}{g^{2}v^{4}}=q^{2}\frac{h^{\prime\,2}(\sigma)}{2\rho^{4}}+\frac{1}{2}(\partial_{\rho}\sigma)^{2}+\frac{\lambda}{2g^{2}}\big{(}\sigma^{2}-1\big{)}^{2}$
(11)

it is clear that in order for $\mathcal{E}$ to be regular at origin, we must
demand $h^{\prime}(0)=0$. In fact, if we assume a polynomial contact
$h^{\prime}\sim h_{0}\sigma^{\alpha}$ as $\sigma\to 0$, it follows that
$\alpha\in[1,2]$ and the asymptotic behaviour of $\sigma$ and $\mathcal{E}$
near the origin falls into two categories:

$\displaystyle 1<\alpha\leq 2$ $\displaystyle\ :$ $\displaystyle\hskip
17.07164pt\left\\{\begin{array}[]{ccl}\sigma&\sim&\displaystyle\bigg{[}\frac{\rho}{qh_{0}(\alpha-1)}\bigg{]}^{\frac{1}{\alpha-1}}\,,\\\\[10.0pt]
\displaystyle\frac{{\mathcal{E}}}{g^{2}v^{4}}&\sim&\displaystyle\bigg{[}\frac{\rho^{2-\alpha}}{qh_{0}(\alpha-1)^{\alpha}}\bigg{]}^{\frac{2}{\alpha-1}}\,,\end{array}\right.\hskip
28.45274pt$ (12c) and $\displaystyle\alpha=1$ $\displaystyle\ :$
$\displaystyle\hskip
17.07164pt\left\\{\begin{array}[]{ccl}\sigma&\sim&\displaystyle\sigma_{0}\exp\\!\bigg{(}-\frac{|qh_{0}|}{\rho}\bigg{)}\,,\\\\[10.0pt]
\displaystyle\frac{\mathcal{E}}{g^{2}v^{4}}&\sim&\displaystyle\frac{q^{2}h_{0}^{2}\sigma_{0}^{2}}{\rho^{4}}\exp\\!\bigg{(}-2\frac{|qh_{0}|}{\rho}\bigg{)}\,.\end{array}\right.\hskip
22.76219pt$ (12f)

In the first category ($1<\alpha\leq 2$) both $\sigma$ and $\mathcal{E}$ fall
off in the origin as a power law with a positive exponent. In the second
category ($\alpha=1$), which is just a limit case $\alpha\to 1$, both $\sigma$
and $\mathcal{E}$ develop non-analytic profiles $\sim\mathrm{e}^{-1/\rho}$. As
we shall see, this behaviour is indicative of a maximal region around
monopole’s core where almost no energy is stored. Thus, this kind of monopole
is “hollow” (a term which we adopt from Bazeia:2018fhg ). It should be also
pointed out, that only in the case $\alpha=2$, the energy density does not
vanish at the origin. In all other cases, the energy is stored in a spherical
shell.

###  II.1 The BPS limit

In the BPS limit $\lambda\to 0$, the energy density can be completed into a
total square as

$\mathcal{E}=\frac{1}{2}\bigg{(}\vec{\nabla}\phi-\frac{h^{\prime}(\phi/v)}{g}\vec{B}\bigg{)}^{2}+\frac{v}{g}\vec{\nabla}\big{(}h(\phi/v)\,\vec{B}\big{)}\,,$
(13)

where we used
$h(\phi/v)\,\vec{\nabla}\cdot\vec{B}=h(0)q\delta^{3}(\vec{r})=0$. Thus, we
demand that the primitive function of $h^{\prime}$ is zero at each zero of
$\phi$, namely $h(0)=0$. Notice that this fixes the constant of integration.
With this additional condition, one can easily derive the formula for mass of
the BPS monopole:

$M_{\rm BPS}\ =\ \frac{4\pi vq}{g}h(1)\ =\ \frac{4\pi
vq}{g}\int_{0}^{1}\\!h^{\prime}(\sigma)\,\mathrm{d}\sigma\,.$ (14)

Let us illustrate the BPS solutions on a particular family of functions

$h^{\prime}(\sigma)=\sigma^{1+1/n}\,,$ (15)

where $n\geq 1$. The BPS equation
$\vec{\nabla}\phi=h^{\prime}(\phi/v)\,\vec{B}/g$ can be easily integrated. For
instance, taking a single monopole at the origin,111Exact multi-monopole
solutions can be readily obtained as well, simply by taking
$\vec{B}=q\sum_{i}(\vec{r}-\vec{r}_{i})/|\vec{r}-\vec{r}_{i}|^{3}$, where the
sum goes over positions $\vec{r}_{i}$ of the monopoles. i.e.,
$\vec{B}=q\,\vec{r}/r^{3}$, the result reads

$\displaystyle\phi=\frac{v}{(1+\frac{q}{n\rho})^{n}}\hskip 14.22636pt$ and
$\displaystyle\hskip
14.22636pt\frac{\mathcal{E}}{g^{2}v^{4}}=\frac{q^{2}}{\rho^{4}(1+\frac{q}{n\rho})^{2n+2}}\,,\hskip
22.76219pt$ (16a) in full accordance with our asymptotic formulas (12c). The
limit $n\to\infty$ takes us directly to the “hollow” monopole case:
$\displaystyle\phi=v\,\mathrm{e}^{-q/\rho}\hskip 14.22636pt$ and
$\displaystyle\hskip
14.22636pt\frac{\mathcal{E}}{g^{2}v^{4}}=\frac{q^{2}}{\rho^{4}}\mathrm{e}^{-2q/\rho}\,,$
(16b)

corresponding to the asymptotic behaviour (12f). We display the energy
densities (16) in Fig. 3. The mass of the monopole follows from (14) as

$M=\frac{4\pi vq}{g}\,\frac{n}{2n+1}\,.$ (17)

Figure 3: Energy density profiles (16) for various $n$. We take $q=1$.

##  III Outlook: Broadening the landscape

At first, it seems that our primitive approach of the previous section runs
orthogonal to the Weinberg–Lee’s ideas of charge screening via dipole moment
interactions. However, there is a natural connection between the two that
results from the following observation. We can rewrite Eq. (8) as

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle-\frac{1}{4g^{2}}F_{\mu\nu}^{2}-\frac{\eta(\phi/v)}{2}d_{\mu\nu}F^{\mu\nu}-g^{2}\frac{\chi(\phi/v)}{4}d_{\mu\nu}^{2}$
(18)
$\displaystyle{}+\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{\lambda}{2}\big{(}\phi^{2}-v^{2}\big{)}^{2}\,.$

Indeed, upon eliminating $d_{\mu\nu}$ via its equation of motion
$\partial_{d_{\mu\nu}}\mathcal{L}=0$, we obtain

$\mathcal{L}=-\frac{1}{4g^{2}}\bigg{[}1-\frac{\eta^{2}(\phi/v)}{\chi(\phi/v)}\bigg{]}F_{\mu\nu}^{2}+\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{\lambda}{2}\big{(}\phi^{2}-v^{2}\big{)}^{2}\,.$
(19)

Hence we return to the original Lagrangian (8) if222One can absorb either
$\eta$ or $\chi$ into the definition of $d_{\mu\nu}$, hence there is no
ambiguity in the correspondence between (18) and (8).

$\eta^{2}=\chi\big{(}1-h^{\prime\,2}\big{)}\,.$ (20)

Furthermore, we can identify $d_{\mu\nu}$ with the dipole moment tensor of the
vector field $W_{\mu}$ in the limit of infinite mass. It is in this sense that
our primitive approach falls into the same purview as the Weinberg–Lee’s
procedure, namely that the finiteness of the monopole’s mass is realized by
screening of the magnetic charge via dipole moment interactions.

There is, however, one important difference. Notice that our approach calls
for field-dependent $\eta$ and $\chi$ couplings that were, however, kept
constant in (7). From the formula (20), we see that at the monopole’s position
$\eta^{2}(0)=\chi(0)$, namely the same condition required by the Weinberg–Lee
scheme for eliminating the singularity in energy density. On the other hand,
in the vacuum $\eta(1)=0$, i.e., the dipole moment interactions are shut down.
This is physically different process then having field-dependent mass of
$W_{\mu}$ fields.

Thus, we see that ours and Weinberg–Lee’s approach are not equivalent, but
they are a special cases of a yet more general procedure that can be subsumed
by the following landscape of $U(1)$ gauge theories:

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle-\frac{f_{1}^{2}(\phi/v)}{4g^{2}}F_{\mu\nu}^{2}-\frac{\eta(\phi/v)}{2}d_{\mu\nu}F^{\mu\nu}-\frac{\chi(\phi/v)g^{2}}{4}d_{\mu\nu}^{2}$
(21)
$\displaystyle{}-\frac{f_{2}^{2}(\phi/v)}{2}\big{|}D_{\mu}W_{\nu}-D_{\nu}W_{\mu}\big{|}^{2}+m^{2}(\phi/v)\big{|}W_{\mu}\big{|}^{2}$
$\displaystyle{}+\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{\lambda}{2}\big{(}\phi^{2}-v^{2}\big{)}^{2}\,.$

In this landscape both approaches are combined: Both mass and dipole couplings
are taken as functions of $\phi/v$ (again, some redundancy exists as one of
those functions can be eliminated by redefinition of $W_{\mu}$ fields).
Furthermore, the kinetic terms of both $A_{\mu}$ and $W_{\mu}$ have field-
dependent pre-factors.

Eq. (21) represents a vast territory. It contains Wu–Yang, ’t Hooft–Polyakov,
Weinberg–Lee’s and our primitive monopoles as a special cases. As already
noted in Lee:1994sk , by adding additional real vector field $Z_{\mu}$
together with appropriate interaction terms, the above theory can be easily
mapped to the bosonic sector of the Electroweak model (or some of its
modification). In this way, both the so-called Cho–Maison monopole Cho:1996qd
and all its finite-mass incarnations Cho:2013vba ; Ellis:2016glu ;
Blaschke:2017pym  are naturally contained in this conceptual scheme.

Further, if we allow for more scalar fields, we can incorporate various non-
linear electrodynamics theories. For instance, to obtain a Born–Infeld
electromagnetism, we simply use a non-dynamical scalar field, say $\varphi$,
as a Lagrange multiplier:

$\displaystyle\mathcal{L}$ $\displaystyle\supset$
$\displaystyle-\frac{\varphi}{4g^{2}}F_{\mu\nu}^{2}-\frac{\beta^{2}}{2}\bigg{(}\varphi+\frac{1}{\varphi}\bigg{)}+2\beta^{2}$
(22) $\displaystyle\xrightarrow[\partial_{\varphi}\mathcal{L}=0]{}$
$\displaystyle\beta^{2}\Bigg{(}1-\sqrt{1+\frac{1}{2g^{2}\beta^{2}}F_{\mu\nu}^{2}}\Bigg{)}\,.$
(23)

Notice that (21) also contains non-Abelian theories with non-minimally coupled
adjoint scalars. This can be seen by writing down such models in unitary gauge
and figuring our relations between field-dependent functions. In this way, the
following sub-landscape can be derived:

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\frac{v^{2}}{2}f_{1}^{2}(\phi/v)\big{(}D_{\mu}\boldsymbol{n}\big{)}^{2}$
(24)
$\displaystyle{}-\frac{1}{4g^{2}}\bigg{[}f_{2}^{2}(\phi/v)\boldsymbol{F}_{\mu\nu}^{2}+f_{3}(\phi/v)\big{(}\boldsymbol{n}\cdot\boldsymbol{F}_{\mu\nu}\big{)}^{2}\bigg{]}$
$\displaystyle{}+\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{\lambda}{2}\big{(}\phi^{2}-v^{2}\big{)}^{2}\,,\hskip
28.45274pt$

where we decomposed the adjoint triplet
$\boldsymbol{\phi}=\phi\,\boldsymbol{n}$ into a scalar field $\phi$ and a unit
iso-vector $\boldsymbol{n}^{2}=1$. Here, we employ $SO(3)$ notation.

We plan to investigate the monopole solutions of these models in a separate
publication.

###### Acknowledgements.

F. B. would like to thank M. Eto and M. Arai for valuable discussions and
express gratitude for the institutional support of the Research Centre for
Theoretical Physics and Astrophysics, Institute of Physics, Silesian
University in Opava and to the Institute of Experimental and Applied Physics,
Czech Technical University in Prague. P. B. is grateful to FSM and Aštar Šeran
for support. This work was supported by the grant of the Ministry of
Education, Youth and Sports of the Czech Republic no. LTT17018.

## References

  * (1) P. A. M. Dirac, _Quantised singularities in the electromagnetic field,_ , _Proc. Roy. Soc. Lond. A_ 133 (1931) 60–72.
  * (2) T. T. Wu and C. N. Yang, _Dirac Monopole Without Strings: Monopole Harmonics_ , _Nucl. Phys. B_ 107 (1976) 365.
  * (3) G. ’t Hooft, _Magnetic Monopoles in Unified Gauge Theories_ , _Nucl. Phys. B_ 79 (1974) 276–284.
  * (4) A. M. Polyakov, _Particle Spectrum in Quantum Field Theory_ , _JETP Lett._ 20 (1974) 194–195.
  * (5) K.-M. Lee and E. J. Weinberg, _Nontopological magnetic monopoles and new magnetically charged black holes_ , _Phys. Rev. Lett._ 73 (1994) 1203–1206, [hep-th/9406021].
  * (6) D. Bazeia, M. A. Marques and G. J. Olmo, _Small and hollow magnetic monopoles_ , _Phys. Rev. D_ 98 (2018) 025017, [1807.01299].
  * (7) Y. M. Cho and D. Maison, _Monopoles in Weinberg-Salam model_ , _Phys. Lett. B_ 391 (1997) 360–365, [hep-th/9601028].
  * (8) Y. M. Cho, K. Kim and J. H. Yoon, _Finite Energy Electroweak Dyon_ , _Eur. Phys. J. C_ 75 (2015) 67, [1305.1699].
  * (9) J. Ellis, N. E. Mavromatos and T. You, _The Price of an Electroweak Monopole_ , _Phys. Lett. B_ 756 (2016) 29–35, [1602.01745].
  * (10) F. Blaschke and P. Beneš, _BPS Cho–Maison monopole_ , _PTEP_ 2018 (2018) 073B03, [1711.04842].


# Eigen-Factors an Alternating Optimization for Back-end Plane SLAM of 3D
Point Clouds

Gonzalo Ferrer1, Dmitrii Iarosh2, Anastasiia Kornilova1 1The authors with the
Skolkovo Institute of Science and Technology (Skoltech), Center for AI
Technology (CAIT). g.ferrer<EMAIL_ADDRESS>author is
with the Software Engineering Department, Saint Petersburg State University.

###### Abstract

Modern depth sensors can generate a huge number of 3D points in few seconds to
be latter processed by Localization and Mapping algorithms. Ideally, these
algorithms should handle efficiently large sizes of Point Clouds under the
assumption that using more points implies more information available. The
Eigen Factors (EF) is a new algorithm that solves SLAM by using planes as the
main geometric primitive. To do so, EF exhaustively calculates the error of
all points at complexity $O(1)$, thanks to the Summation matrix $S$ of
homogeneous points.

The solution of EF is highly efficient: i) the state variables are only the
sensor poses – trajectory, while the plane parameters are estimated previously
in closed from and ii) EF alternating optimization uses a Newton-Raphson
method by a direct analytical calculation of the gradient and the Hessian,
which turns out to be a block diagonal matrix. Since we require to
differentiate over eigenvalues and matrix elements, we have developed an
intuitive methodology to calculate partial derivatives in the manifold of
rigid body transformations $\operatorname*{\mathnormal{SE}(3)}3$, which could
be applied to unrelated problems that require analytical derivatives of
certain complexity.

We evaluate EF and other state-of-the-art plane SLAM back-end algorithms in a
synthetic environment. The evaluation is extended to ICL dataset (RGBD) and
LiDAR KITTI dataset. Code is publicly available at https://github.com/prime-
slam/EF-plane-SLAM.

###### Index Terms:

SLAM, Plane SLAM, Planar Bundle Adjustment, Differentiation in
$\operatorname*{\mathnormal{SE}(3)}3$.

##  I Introduction

The advent of modern sensing capabilities is offering the possibility to
perceive 3D depth, i.e. Point Clouds (PC), at a high streaming rate from
diverse sensors, such as 3D LiDARs, Time of Flight cameras, RGB-D, etc. The
more points are used, the more information is gathered from the environment.
However, is it possible to evaluate all observed points and improve accuracy
without a considerable increase in complexity? This is the main question this
paper tries to address.

Aligning a sequence of observed PCs can be formulated as a Simultaneous
Localization and Mapping (SLAM) [1] problem, where full trajectory and map
geometry are purely estimated from the sensor’s observed 3D points. Similarly,
this problem definition fits perfectly the Bundle Adjustment (BA) [2] problem.
In this paper, we are interested in the design and study of a new observation
function that allows to formulate the SLAM (or BA) problem by directly
considering the sensor’s observed points in an efficient manner.

An important hindrance is that the sensor’s points from PCs are sampled from
the scene geometry, consequently it is highly unlikely to observe the same
exact point on later observations. Therefore, one needs to define some
representation of this inherent geometry that is capable of providing the
error-point in an efficient manner. We will choose planar surfaces as
landmarks due to their regularity and simple constrain equation. These
geometrical surfaces are commonly present in indoors and urban areas. Plane
SLAM [3, 4] approaches use planes as landmarks, but the observations are pre-
processed into planar features, hence, 3D points are not directly used with
their consequent performance deterioration.

We propose a plane SLAM technique for back-end optimization, the Eigen-Factors
(EF), tailored for 3D Point Clouds (PC) from RGBD and LiDAR, while being
capable of processing a large amount of points with constant complexity. This
paper is an improved version of the former EF [5], overcoming its main
limitation (first order optimization). The EF uses planar geometric
constraints for the alignment of a trajectory described by a moving sensor. To
do so, EF aggregates all observed points into the Summation matrix $S$, a
summation of outer product of homogeneous points and directly calculates the
exact point error at $O(1)$ complexity.

In addition, EF optimizes the point to plane error w.r.t. the trajectory. The
plane variables are calculated in closed form in an early step and the
resultant problem is an alternating optimization.

The gradient and the Hessian are directly calculated by reformulating the
problem in terms of manifold optimization on
$\operatorname*{\mathnormal{SE}(3)}3$ and retractions. Unfortunately, EF does
not allow for the Gauss-Newton optimization approximation, since there are not
residuals but eigenvalues, so instead we use the Newton-Raphson optimization.
In order to differentiate over $\operatorname*{\mathnormal{SE}(3)}3$, we have
proposed a methodology which greatly simplifies the calculation of derivatives
of any order and it could be used on other problems where an analytical
expression is hard to derive, as in our case through eigenvalues or matrix
elements. For our particular configuration, the resulting Hessian is block
diagonal which has an impact on the optimization process, since its inversion
is extremely efficient.

The contributions of the Eigen Factors are the following:

  * •

Matrix summation $S$ to calculate all point to plane errors at $O(1)$
complexity;

  * •

Alternating optimization where first the plane parameters are calculated in
closed form and then the reduced dimensionality problem (only trajectory) is
estimated by Newton-Raphson;

  * •

Efficient second order optimization of EF with analytical calculation of the
gradient and the Hessian (block diagonal).

  * •

Methodology to calculate derivatives of any function involving the manifold of
rigid body transformations $\operatorname*{\mathnormal{SE}(3)}3$;

The evaluation compares EF with other state-of-the-art plane SLAM back-end
algorithms in synthetic, RGBD and LiDAR environments, since we believe these
are the most useful settings for plane SLAM. In order to clearly understand
all the details of each method, we have decided to compare only the back-end
and provide to all methods the same exact input. The front-end has an impact
to the final performance of the methods so we have decided to isolate this
effect.

The ideal application of EF is map refinement of 3D PC by optimizing the plane
error over trajectories.

##  II Related Work

The initial motivation for EF is to align PC minimizing the point error. The
Iterative Closest Point (ICP) algorithm [6, 7, 8, 9] fits a pair of PCs by
iteratively finding correspondences and solving the optimization problem [10,
11, 12]. Existing solutions include matching point-to-point, point-to-line
[13], point-to-plane, and other variants emulating planes [14, 15] or assuming
a general geometry. These techniques require to process all points every
iteration and can achieve high accurate results.

One can bring this idea to a sequence of PC or Multi-view registration. An
early work [16] proposed to match all points in every view, reducing the
overall point error, but requiring huge memory resources and intense
computations. Improvements on Multi-view registration targeted efficiency, for
instance considering pairs of PC as constraints [17].

Although computing the overall point error is a natural choice for mapping and
SLAM, researchers have been trying to reduce the a priory high computational
requirements of these techniques. On this regard, processing plane landmarks
for each PC observation could be seen as a proxy for computing the point error
but at a reduced computational cost. Therefore, plane Landmark-based SLAM
optimizes trajectory and a map of plane-landmarks [3, 4]. There have been
proposed multitude of variants, some of them focusing on the plane
representation [18, 19], adding planar constraints [20] or a mixture of points
and planes [21, 22]. The key problem of plane landmarks is that after
extracting planar features from noisy points, there is an inevitable loss of
information compared to considering jointly all points. EF on the other hand,
keeps all points’ information.

Modern mapping approaches, in order to ensure real-time (or close) conditions,
aggregate observations to the map after an accurate scan-to-map, effectively
marginalizing poses. Examples of these approaches include LOAM [23, 24] which
combines PC alignment with mapping or LIO-mapping [25]. Surfel-based
approaches [26, 27] follow the same principle dividing the scene in planar
disks, obtaining impressive results. When considering loop closure, pure
mapping approaches present limitations and require an extra effort to update
their mapping solutions properly. SLAM approaches do not suffer this
limitation, since in general both the trajectory and the map are estimated
simultaneously, so they handle more naturally loop closure situations at the
cost of higher computational resources.

When considering overall point error, efficiency is of the essence. Our
previous paper on EF [5] introduces for the first time the $S$ matrix as a
summation of outer products the homogeneous points from each sensor pose. This
allows to optimize the point to plane error with complexity $O(1)$, which is
the key requirement for our approach to work efficiently. Later, $\pi$-LSAM
[28] makes use of the same matrix $S$ to calculate the overall point to plane
error to solve the Planar Bundle adjustment estimating jointly poses and
plane-landmarks. In a later work, the authors refer to this homogeneous
summation matrix $S$ as the Integrated Cost Matrix (ICM) [22].

The BAREG [29] algorithm updates incrementally a covariance matrix $3\times
3$, without explicit re-evaluation of points, maintaining the same advantages
as the $S$ matrix. BALM2 [30] uses the same concept of the $S$ matrix, which
is referred there as Cluster Points and formalizes some properties of the $S$
matrix such as superposition and transformation.

The idea of using the eigendecomposition from EF as an objective function is
adopted by BALM [31], whose derivation is from a $3\times 3$ covariance
instead of the $4\times 4$ homogeneous $S$ matrix as in EF. The authors
improve the efficiency [30] by calculating analytical derivatives of the
gradient and Hessian, in the same spirit we propose in this paper. The
derivations are however different: BALM’s Hessian is dense whereas our
derivation results in a block diagonal Hessian.

In BAREG[29], the authors use the eigendecomposition to demonstrate the
equivalence of the point-to-plane error to the plane-to-plane error, with
proper weights. This equivalence assumes that the plane parameters are optimal
which might be a too restrictive condition in practice. However, this
reduction greatly reduces the complexity and efficiency of the optimization
step. It also requires the plane estimation at each observation, which has
consequences when evaluating under strong point noise conditions.

In this work, we are interested also in assessing the quality of the map by
using a mapping metric. The ground truth on some datasets will not be accurate
enough for the mapping purposes, so these assessment tools are necessary.
Whereas there are many solutions for precise trajectory estimation from
sequences of point clouds, the topic of measuring map quality has not received
the same attention. The 3D graphics community proposes a set of full reference
metrics for solving tasks of point cloud compression and denoising, including
p2point metric, p2plane metric [32], angular similarity [33], projection based
methods [34] or SSIM [35, 36]. The majority of metrics based on point cloud
requires reference to ground truth, whereas in the evaluation scenario we
propose (KITTI), 3D reconstruction ground truth will not be available.
Therefore no-reference evaluation is used in this work. In particular, Mean
Map Entropy (MME) [37] and Mean Plane Variance (MPV) [38] evaluate noise in
the aggregated point cloud and Mutually Orthogonal Metric (MOM) [39]
correlates reconstruction noise with perturbation in the 3D poses.

##  III Eigen-Factors: Point to Plane Error

A plane $\pi$ in the 3D space is defined by a normal unit vector
$\eta\in\mathbb{S}^{2}$ and the plane distance to the origin $d\in\mathbb{R}$:

$\pi=\begin{bmatrix}\eta\\\ d\end{bmatrix}\in\mathbb{P}^{3},\qquad\text{where
}||\eta||=1.$ (1)

A point $p=[x,y,z]^{\top}$ lies on the plane $\pi$ if an only if:

$\pi^{\top}\begin{bmatrix}p\\\ 1\end{bmatrix}=\pi^{\top}\tilde{p}=0.$ (2)

With this constrain and given a set of noisy observed points from the same
planar surface, one can solve the plane estimation mainly by two approaches.
The centered method minimizes the overall point to plane error

$\min_{\pi}\sum_{n=1}^{N}||\eta^{\top}p_{n}+d||^{2}.$ (3)

The solution comes in two steps:

$\displaystyle\min_{\eta}\left\\{\sum_{n=1}^{N}\eta^{\top}\underbrace{(p_{n}-E\\{p\\})(p_{n}-E\\{p\\})}_{\Sigma_{p}}\eta\right\\}$
(4) $\displaystyle d=-E\\{\eta^{\top}p\\},$ (5)

first the eigendecomposition of the $3\times 3$ matrix $\Sigma_{p}$ solves the
value of $\eta$ and then $d$ is obtained.

The homogeneous method calculates the point to plane error of a plane $\pi$
without centering the data points:

$\min_{\pi}\sum_{n=1}^{N}||\pi^{\top}\tilde{p}_{n}||^{2}=\pi^{\top}\tilde{P}\tilde{P}^{\top}\pi$
(6)

where the matrix $\tilde{P}$ is the stacked vector of $N$ homogeneous points.
This arrangement of terms in the homogeneous plane estimation method allows us
to define the $S$ matrix:

$S=\tilde{P}\tilde{P}^{\top}=\sum_{n=1}^{N}\tilde{p}_{n}\tilde{p}_{n}^{\top},$
(7)

which is the Summation of the outer product of all points in the plane, hence
matrix $S$. The solution of (6) is calculated using the eigendecomposition of
the $4\times 4$ matrix $S$. In [40] there is a comparison of the two methods
yielding similar results.

Our ultimate objective is the optimization of the trajectory and the
geometrical features in the map. The disadvantage of the centered method is
that every time a new sample is added or modified, all calculations should be
carried out again. It is more natural to develop the homogeneous method to
account for poses in a highly efficient solution.

Then, the Eigen-Factors is a generalization of the homogeneous plane
estimation observed over a sequence of reference frames or poses
$T_{t}\in\operatorname*{\mathnormal{SE}(3)}3$ at different instants of time
$t\in\\{1,\ldots,H\\}$. It can also be understood as a trajectory
$\bm{T}=\\{T_{1},\ldots,T_{H}\\}$. Each of these poses transforms points to a
global reference frame $g$, such that ${}^{g}T_{t}$, where we will omit the
global frame reference for simplicity in following sections. The underlying
idea is that the same plane $\pi$ is observed from different poses $T_{t}$ but
the planar constraint must hold for all views. Accordingly, we can reformulate
the plane estimation problem for the plane $\pi$ as

$\min_{\pi}\sum_{t=1}^{H}\sum_{n=1}^{N_{t}}||\pi^{\top}T_{t}\tilde{p}_{t,n}||^{2}=\min_{\pi}\pi^{\top}\big{(}\sum_{t=1}^{H}T_{t}S_{t}T_{t}^{\top}\big{)}\pi,$
(8)

where each of the points $\tilde{p}_{t,n}$ denotes the reference frame it was
observed $t$ and the matrix $S_{t}$ includes all these $N_{t}$ points as in
(7). By re-arranging terms one can see as the plane $\pi$ is out of the
summation, so the plane solution boils down to the eigendecompsition of this
new matrix. Therefore, we define the variable $Q$ as the summation

$Q(\bm{T})=\sum_{t=1}^{H}T_{t}S_{t}T_{t}^{\top}=\sum_{t=1}^{H}Q_{t}(T_{t}),$
(9)

where each component $Q_{t}$ depends on the matrix $S_{t}$, calculated only
once, and the current 3D pose estimation $T_{t}$. The matrix $S_{t}$ expresses
points in the local coordinate system at time $t$ and is constant, while the
matrix $Q$ is expressed in the global reference frame and is recalculated for
each new updated trajectory.

The solution of the minimization problem is obtained directly by the
eigendecomposition:

$\displaystyle\pi^{*}=k\cdot v_{min}(Q(\bm{T}))$ (10) $\displaystyle\text{s.t.
}||\eta^{*}||=1.$

The minimum eigenvector $v_{min}$ is proportional to the plane solution. The
exact solution, in order to have geometrical meaning, must fulfill the
definition of a plane (1) with unit normal $\eta$. Therefore, it is multiplied
by the scalar value $k$. Without this correction, each of the errors would be
scaled differently and the solution would be biased. The summation of the
point to plane squared error (for plane $\pi$) is equal to:

$\min_{\pi}\pi^{\top}Q(\bm{T})\pi=k^{2}\lambda_{min}(Q(\bm{T}))=\lambda_{\pi}(\bm{T})$
(11)

where $\lambda_{min}$ is the minimum eigenvalue, which equals the minimum cost
or minimum error after estimation after scaling by $k^{2}$, easily derived
from (10). Note that the minimization of the plane is implicit after solving
the algebraic equation by the eigendecomposition, so the plane variables can
be removed from the state variables.

The function $\lambda_{\pi}(\cdot)$ obtained in (11) for the plane $\pi$ shows
some interesting properties. In the ideal case of no-noise, this function
should be zero and its values are always positive by construction
$\lambda_{\pi}(\cdot)\geq 0$. We could reformulate this as an implicit
observation function of the form $h(x,z)=0$, where in our particular case, the
state $x$ is the trajectory $\bm{T}$ and the observation $z$ are all points.
In contrast, an explicit observation is of the form $h(x)=z$ which when
optimizing tries to fit the state with the observations.

Unfortunately, this implicit observation function is not suitable for general
nonlinear least squares (NLLS). We will later derive a state estimation
problem from $\lambda_{\pi}(\cdot)$, where the total cost to minimize is

$C(\bm{T})=\sum_{m=1}^{M}\lambda_{\pi_{m}}(\bm{T})$ (12)

given that $M$ is the total number of observed planes. Some important remarks:

  * •

The cost $C(\bm{T})$ is equal to the total point-to-plane squared error
evaluated for each of the poses $T_{t}$.

  * •

The complexity is independent of the number of points, since the matrix
$S_{t}$ requires them to be summed only once and use thereafter without
modifications.

  * •

The plane parameters are estimated in closed form and the state variables are
the sequence of poses $\bm{T}$.

This is the essence of the EF, a function that seeks to minimizes the point to
plane error w.r.t. the trajectory. The next sections are dedicated to the
particular optimization of this function that requires specific development.

Our approach is an Alternating optimization, where the plane parameters $\pi$
are obtained in closed form by the eigendecomposition and then the reduced
problem is estimated by optimizing w.r.t. the trajectory $\bm{T}$ by assuming
$\pi$ to be a constant, such that

$\frac{\partial\pi}{\partial\bm{T}}\approx 0.$ (13)

Other alternatives to the calculation of separable variables include the
variable projection [41] or the general Wiberg approach [42] for which (13)
does not hold. We have chosen the alternating optimization due to its
simplicity and efficiency, since the main purpose of EF is map refinement to
improve the quality of the map.

##  IV Retraction and Optimization on the Manifold

This section introduces some tools from optimization on matrix manifolds [43]
required for the correct development of a second-order optimization method for
the Eigen Factors. The problem of interest is to further develop the form
$Q(T):\operatorname*{\mathnormal{SE}(3)}3\to\mathbb{R}^{4\times 4}$, so we
will propose systematic rules to differentiate any function involving
$\operatorname*{\mathnormal{SE}(3)}3$ in a simple way.

A Rigid Body Transformation (RBT), denoted as
${}^{f}T_{o}\in\operatorname*{\mathnormal{SE}(3)}3$, allows to transform
points from an origin reference system $o$ to a final reference system $f$. It
also transforms reference frames or simply expresses a 3D pose uniquely. More
formally, the Special Euclidean group is

$SE(3)=\left\\{T=\begin{bmatrix}\mathbf{R}&\bm{t}\\\
0&1\end{bmatrix}\,|\,\mathbf{R}\in
SO(3)\,\land\,\bm{t}\in\mathbb{R}^{3}\right\\}.$ (14)

In addition, RBT is a matrix group
$\operatorname*{\mathnormal{SE}(3)}3\subset\mathbb{R}^{4\times 4}$ and due to
its inner structure and properties is a Riemannian manifold. As such, there is
a smooth mapping (chart) from a subgroup of
$\operatorname*{\mathnormal{SE}(3)}3$ elements to the vector space
$\mathbb{R}^{6}$. It is well known that the dimensionality of a 3D pose is
six, however, there exist several possible mappings.

Our interest is on the optimization of a cost function whose input is an
element of $\operatorname*{\mathnormal{SE}(3)}3$ or more generally, an element
of a manifold $\mathcal{M}$. Let the function $f(x):\mathcal{M}\to\mathbb{R}$
be a mapping from a manifold to a real value. Then, if one applies directly
the definition of directional derivative

$\lim_{h\to 0}\frac{f(x+h\eta)-f(x)}{h},$ (15)

finds this operation not suitable for numerical calculations.

One can solve this problem by using a retraction mapping.

######  Definition 1.

A retraction $R_{T}(\xi)$ is a smooth mapping from the tangent space around
the $T$ element to the manifold:

$R_{T}(\xi):\mathbb{R}^{6}\to\operatorname*{\mathnormal{SE}(3)}3.$ (16)

Two conditions must be satisfied: i) $R_{T}(0_{T})=T$ and ii) local rigidity
(see Ch.4 [43]).

Accordingly, we can compose a new function when using a retraction:

$\hat{f}_{T}=f\circ R_{T},$ (17)

from the tangent space $\mathbb{R}^{6}$ to a real value. This new function
$\hat{f}_{T}$ admits a directional derivative (15). As a convention in this
paper, all retractions are denoted with the hat superscript.

As a result, an optimization problem $x^{\prime}=x+\Delta x$, which is not
well defined for elements of a manifold $x\in\mathcal{M}$, becomes well
defined when using a retraction

$T^{\prime}=R_{T}(\Delta\xi).$ (18)

In the case of first order optimization, the update is proportional to the
gradient direction $\Delta\xi=-\alpha\nabla_{\xi}\hat{f}_{T}$. Gradient decent
is one of the simplest methods for optimization: it requires a direction and a
magnitude. Second-order optimization methods, provide the direction of descent
$\Delta\xi=-\alpha(\nabla_{\xi}^{2}\hat{f}_{T})^{-1}\nabla_{\xi}\hat{f}_{T}$.
The advantage is that the Hessian already weights this descent direction,
resulting in improved convergence rates.

In the particular case of $\operatorname*{\mathnormal{SE}(3)}3$, its tangent
space has a structure of Euclidean space. Only under these conditions, the
gradient and the Hessian are well defined after using a retraction:

$\displaystyle\nabla_{T}f(T)=\nabla_{\xi}\hat{f}_{T}(0_{T})$
$\displaystyle\nabla_{T}^{2}f(T)=\nabla_{\xi}^{2}\hat{f}_{T}(0_{T}).$

These results are significant since they will allow to use the concept of
directional derivatives (15) by simply applying real analysis into the
calculation of the gradient and Hessian in Sec.V. By doing that, we will
implement a second-order optimization method that have almost no need for
hyper-parameters and quadratic convergence rate [44].

###  IV-A Differentiating over the Exponential Map in
$\operatorname*{\mathnormal{SE}(3)}3$

The exponential map is the most natural retraction choice: its definition is
highly related with the concept of perturbations over
$\operatorname*{\mathnormal{SE}(3)}3$ which will facilitate the calculation of
derivatives in the following sections. There are other options for choosing a
retraction mapping, as long as they satisfy the definition.

The exponential map defines the retraction $R_{T}$ in the following way:

$R_{T}(\xi)=\operatorname*{\mathrm{Exp}}(\xi)T.$ (19)

We follow a left-hand-side convention, but the right-hand-side is a valid
retraction and it would require a similar derivation of the terms. The
retraction around the identity can be expanded by using the matrix exponential
definition as a Taylor series

$R_{I}(\xi)=\operatorname*{\mathrm{Exp}}(\xi)I=I+\xi^{\land}+\frac{1}{2}(\xi^{\land})^{2}+o(|\xi|^{3}),$
(20)

where the vector
$\xi=[\theta\operatorname*{{}^{\top}},\rho\operatorname*{{}^{\top}}]\operatorname*{{}^{\top}}=[\xi_{1},\ldots,\xi_{6}]\operatorname*{{}^{\top}}\in\mathbb{R}^{6}$
and the matrix

$\displaystyle\xi^{\land}=\begin{bmatrix}\operatorname*{\theta^{\land}}&\rho\\\
0&0\end{bmatrix}=\begin{bmatrix}0&-\theta_{3}&\theta_{2}&\rho_{1}\\\
\theta_{3}&0&-\theta_{1}&\rho_{2}\\\ -\theta_{2}&\theta_{1}&0&\rho_{3}\\\
0&0&0&0\end{bmatrix}$ (21)

is the Lie Algebra matrix group $\mathfrak{se}(3)$ that represents the group
of infinitesimal RBT around the identity. See [45, 46] for a comprehensive
explanation on the topic of RBT and its Lie algebra. This matrix can be
rearranged such that it has a vector space structure

$\xi^{\land}=G_{1}\xi_{1}+G_{2}\xi_{2}+\ldots+G_{6}\xi_{6},$ (22)

where $G_{i}$ is a $4\times 4$ matrix of the i-$th$ generator and the basis
$\bm{G}=\\{G_{i}\\},\forall i=1,\ldots,6$.

After defining the retraction map (19), the next step is to obtain the
derivative of the retraction, i.e., the exponential map evaluated at the zero
element.

######  Lemma 1.

The derivative of the exponential map with respect to one of the Lie
coordinates $\xi_{i}$ is its corresponding matrix generator

$\frac{\partial\operatorname*{\mathrm{Exp}}(\xi)}{\partial\xi_{i}}=G_{i}.$
(23)

###### Proof.

One can expand analytically the definition of the matrix exponent in (20)

$\displaystyle\frac{\partial\operatorname*{\mathrm{Exp}}(\xi)}{\partial\xi_{i}}$
$\displaystyle=\frac{\partial}{\partial\xi_{i}}\Big{(}I+\xi^{\land}+\frac{1}{2}(\xi^{\land})^{2}+o(|\xi|^{3})\Big{)}\bigg{|}_{\xi_{i}=0}$
$\displaystyle=G_{i}+\frac{1}{2}G_{i}\xi^{\land}+\frac{1}{2}\xi^{\land}G_{i}+\ldots\bigg{|}_{\xi_{i}=0}=G_{i}$
(24)

where after differentiating each of the terms, the higher order terms vanish
when evaluated at $\xi_{i}=0$. ∎

######  Lemma 2.

The second derivative of the exponential with respect to the Lie coordinates
$\xi_{i}$ and $\xi_{j}$ is

$\frac{\partial^{2}\operatorname*{\mathrm{Exp}}(\xi)}{\partial\xi_{j}\partial\xi_{i}}=\frac{1}{2}G_{i}G_{j}+\frac{1}{2}G_{j}G_{i}.$
(25)

###### Proof.

one can repeat the same procedure to expand analytically the second-order
derivative:

$\displaystyle\frac{\partial^{2}\operatorname*{\mathrm{Exp}}(\xi)}{\partial\xi_{j}\partial\xi_{i}}$
$\displaystyle=\frac{\partial^{2}}{\partial\xi_{j}\partial\xi_{i}}\Big{(}I+\xi^{\land}+\frac{1}{2}(\xi^{\land})^{2}+o(|\xi|^{3})\Big{)}\bigg{|}_{\xi=0}$
$\displaystyle=\frac{\partial}{\partial\xi_{j}}\Big{(}G_{i}+\frac{1}{2}G_{i}\xi^{\land}+\frac{1}{2}\xi^{\land}G_{i}\Big{)}\bigg{|}_{\xi=0}$
$\displaystyle=\frac{1}{2}G_{i}G_{j}+\frac{1}{2}G_{j}G_{i}$ (26)

∎

This compact result is a $4\times 4$ matrix, corresponding to each of the
$(i,j)$ elements of the Hessian. The overall tensor
$\frac{\partial^{2}\operatorname*{\mathrm{Exp}}(\xi)}{\partial\xi^{2}}$ has
dimensions $4\times 4\times 6\times 6$ when considering each of its
coordinates $(i,j)$.

The idea of retraction functions plus the first and second-order derivatives
of the exponential map will be enormously useful when calculating derivatives
of more complex functions, such as the matrices in the EF, as we show in the
next section.

##  V Eigen-Factors Differentiation

We will start by redefining the key functions, now as a retraction function
(denoted by the $\hat{\cdot}$ indicator). The $Q$ matrix in (9) was defined as
a function of the trajectory $\bm{T}$. If only considering the transformation
$T_{t}$ at time $t$, then $Q=Q_{t}$ can be redefined as the retraction
function (17) of the form $\hat{Q}_{t}=Q_{t}\circ
R_{T_{t}}:\mathbb{R}^{6}\to\mathbb{R}^{4\times 4}$, such that

$\hat{Q}_{t}(\xi_{t})=\operatorname*{\mathrm{Exp}}(\xi_{t})\cdot
Q_{t}\cdot\operatorname*{\mathrm{Exp}}(\xi_{t})^{\top}.$ (27)

The point to plane error, defined as the $\lambda_{\pi}(\bm{T})$ in the EF
derivation in (11), can be as well redefined as the retraction function

$\hat{\lambda}_{\pi}(\xi_{t})=\lambda_{\pi}\circ\hat{Q}_{t}:\mathbb{R}^{6}\to\mathbb{R}$
(28)

which is now a well defined function to calculate directional derivatives
(15). One can evaluate the retraction around any other pose in a more general
form $\hat{\lambda}_{\pi}(\bm{\xi})$, however, for the development that
follows, we will use the current time index $t$ and generalize the obtained
expressions later to all the trajectory.

Now, it is necessary to propagate the gradient through the eigendecomposition.

######  Proposition 1.

The derivative of the eigenvalue $\lambda_{\pi}$ with respect to the pose
$T_{t}$ equals

$\nabla_{\xi_{t}}\hat{\lambda}_{\pi}=\frac{\partial\hat{\lambda}_{\pi}}{\partial\xi_{t}}=\pi^{\top}\frac{\partial\hat{Q}_{t}}{\partial\xi_{t}}\,\pi.$
(29)

###### Proof.

See the analytical derivation in the Appendix A-A. ∎

Accordingly, the calculation of the gradient boils down to calculate the
derivative of $\hat{Q}_{t}$ with respect to the tangent space around the
element $T_{t}$ of the manifold. We differentiate (27) such that

$\frac{\partial\hat{Q}_{t}}{\partial\xi_{t}}=\frac{\partial\text{Exp}(\xi_{t})}{\partial\xi_{t}}Q_{t}+Q_{t}\frac{\partial\text{Exp}(\xi_{t})^{\top}}{\partial\xi_{t}},$
(30)

which is obtained by applying the product rule.

By using Lemma 1, one simply substitutes to obtain the derivative wrt the
coordinate $i$:

$\frac{\partial\hat{Q}_{t}}{\partial\xi_{t,i}}=G_{i}\cdot Q_{t}+Q_{t}\cdot
G_{i}^{\top},$ (31)

where each of these elements is a $4\times 4$ matrix. The gradient is then
obtained after multiplying this result as in (29) which yields a scalar value;
in total a 6D vector for each of the $i=1,\ldots,6$ components.

This methodology is general for obtaining derivatives of functions wrt poses
in 3D and allows to have a compact result, even when the image of this
function is a matrix.

###  V-A Eigen Factors, Hessian

The objective of this section is to obtain an analytical solution for the
Hessian matrix $\nabla_{\xi}^{2}\hat{\lambda}$. Similarly, we follow the same
systematic procedure to derive the Hessian of the EF.

######  Proposition 2.

The Hessian of the Eigendecomposition with respect to the pose $T_{t}$ is

$\displaystyle\frac{\partial^{2}\hat{\lambda}_{\pi}}{\partial\xi_{t}^{2}}=\pi^{\top}\frac{\partial^{2}\hat{Q_{t}}}{\partial\xi_{t}^{2}}\pi$
(32)
$\displaystyle\frac{\partial^{2}\hat{\lambda}_{\pi}}{\partial\xi_{t}\partial\xi_{t^{\prime}}}=0,\qquad\forall
t\neq t^{\prime}.$ (33)

###### Proof.

See the analytical derivation in the Appendix A-B. ∎

An important property derived from Proposition 2 is that the Hessian matrix is
block diagonal, which implies a fast calculation of the optimization
iterations and a reduced complexity of the algorithm.

Now, we will derive the second derivative for the pair $(\xi_{t,i},\xi_{t,j})$
of components of the vector $\xi_{t}$. We calculate the derivative of
$\hat{Q}_{t}$ in (27) twice, by using the product rule:

$\displaystyle\frac{\partial^{2}\hat{Q}_{t}}{\partial\xi_{t,j}\partial\xi_{t,i}}=\frac{\partial}{\partial\xi_{t,j}}\left(\frac{\partial\text{Exp}(\xi_{t})}{\partial\xi_{t,i}}Q_{t}+Q_{t}\frac{\partial\text{Exp}(\xi_{t})^{\top}}{\partial\xi_{t,i}}\right)$
$\displaystyle=\underbrace{\frac{\partial^{2}\text{Exp}(\xi_{t})}{\partial\xi_{t,j}\partial\xi_{t,i}}Q_{t}+\frac{\partial\text{Exp}(\xi_{t})}{\partial\xi_{t,i}}\frac{\partial
Q_{t}}{\partial\xi_{t,j}}}_{B_{ij}}+$
$\displaystyle+\underbrace{\frac{\partial
Q_{t}}{\partial\xi_{t,i}}\cdot\frac{\partial\text{Exp}(\xi_{t})^{\top}}{\partial\xi_{t,j}}+Q_{t}\frac{\partial^{2}\text{Exp}(\xi_{t})^{\top}}{\partial\xi_{t,j}\partial\xi_{t,i}}}_{B_{ij}^{\top}}.$
(34)

Given the symmetry in $Q=Q^{\top}$, this expression is simplified and one has
to calculate only the expression

$B_{ij}=(G_{i}G_{j}+G_{j}G_{i})\frac{Q_{t}}{2}+G_{i}\frac{\partial
Q_{t}}{\partial\xi_{t,j}},$ (35)

and sum to its transpose in order to obtain the the components $(i,j)$ of the
Hessian. Note that the Hessian is symmetric as well, and only the upper-
triangular elements need to be calculated.

###  V-B Trajectory Optimization

In this section, we discuss different options for optimizing a trajectory
$\bm{T}=\\{T_{1},\ldots,T_{H}\\}$ or its corresponding set of vectors
$\bm{\xi}=\\{\xi_{1},\ldots,\xi_{H}\\}$. One can directly apply the EF
gradient and Hessian at each pose in the trajectory, resulting in a the joint
vector gradient

$\frac{\partial\hat{\lambda}_{\pi}}{\partial\bm{\xi}}=\left[\Big{(}\pi^{\top}\frac{\partial
Q}{\partial\xi_{1}}\pi\Big{)}^{\top},\ldots,\Big{(}\pi^{\top}\frac{\partial
Q}{\partial\xi_{H}}\pi\Big{)}^{\top}\right]^{\top}_{6H\times 1}.$ (36)

This aggregation corresponds to the Cartesian product of the manifolds of each
of the poses. In case there were not points observed from a pose, then the
gradient would zero since it is not contributing to the EF cost in any mean.
The same arrangement can be done for the Hessian:

$\frac{\partial^{2}\hat{\lambda}_{\pi}}{\partial\bm{\xi}^{2}}=\text{diag}\left[\Big{(}\pi^{\top}\frac{\partial^{2}Q}{\partial\xi_{1}^{2}}\pi\Big{)},\ldots,\Big{(}\pi^{\top}\frac{\partial^{2}Q}{\partial\xi_{H}^{2}}\pi\Big{)}\right]_{6H\times
6H}.$ (37)

The total EF cost (12) can be redefined in terms of retraction function w.r.t.
the joint vector of poses $\bm{\xi}$

$\hat{C}(\bm{\xi})=\sum_{m=1}^{M}\hat{\lambda}_{\pi_{m}}(\bm{\xi}).$ (38)

The optimal update of the joint vector of pose coordinates is solved by

$\displaystyle\Delta\bm{\xi}=-\alpha(\nabla_{\xi}^{2}\hat{C})^{-1}\nabla_{\xi}\hat{C}$
(39)
$\displaystyle\nabla_{\xi}\hat{C}=\sum_{m}^{M}\frac{\partial\hat{\lambda}_{\pi_{m}}}{\partial\bm{\xi}}$
(40)
$\displaystyle\nabla_{\xi}^{2}\hat{C}=\sum_{m}^{M}\frac{\partial^{2}\hat{\lambda}_{\pi_{m}}}{\partial\bm{\xi}^{2}}.$
(41)

###  V-C Error Invariance to the Reference Frame

Planes are elements of the projective space $\mathbb{P}^{3}$ and they can be
transformed from one reference frame to another, similarly to homogeneous
points by RBT:

${}^{t}\pi={}^{t}T_{g}^{-\top}\cdot{}^{g}\pi,$ (42)

where the plane ${}^{g}\pi$ in the global reference frame is transformed into
${}^{t}\pi$ in the local frame $t$.

If we explicitly account all reference frames in the point to plane error in
(8), where the point ${}^{t}p$ was expressed in a local coordinate and the
plane was expressed in the global coordinate $g$, then

${}^{g}\pi^{\top}\cdot{}^{g}T_{t}\cdot{}^{t}\tilde{p}={}^{g}\pi^{\top}\cdot{}^{g}\tilde{p}={}^{t}\pi^{\top}\cdot{}^{t}\tilde{p}.$
(43)

One can see that the error value is invariant to the choice of the reference:
we could either transform both elements to the local reference frame or to the
global frame. The EF development used the latter.

Not just that, but the planar constraint is fulfilled in any reference frame:

${}^{c}\pi^{\top}\cdot{}^{c}\tilde{p}=({}^{c}T_{g}^{-\top}{}^{g}\pi)^{\top}\cdot{}^{c}T_{g}{}^{g}\tilde{p}={}^{g}\pi^{\top}\cdot{}^{g}\tilde{p},$
(44)

where the transformations cancel out for any new reference frame $c$.

###  V-D Data Centering

From the last subsection, we can express the EF cost in a different reference
frame by multiplying a transformation ${}^{c}T_{g}$ (for simplicity we write
$T_{c}$) and its inverse, so the the error is unaltered:

$\displaystyle\lambda_{\pi}$
$\displaystyle=\pi^{\top}\hat{Q}\pi=\pi^{\top}\underbrace{T_{c}^{-1}T_{c}}_{I}\cdot\hat{Q}\cdot
T_{c}^{\top}\underbrace{T_{c}^{-\top}\pi}_{{}^{c}\pi}$
$\displaystyle={}^{c}\pi^{\top}T_{c}\cdot\hat{Q}\cdot T_{c}^{\top}{}^{c}\pi.$
(45)

The solution is the well known eigendecomposition of $T_{c}\cdot\hat{Q}\cdot
T_{c}^{\top}$. When applied to the gradient, yields

$\frac{\partial\hat{\lambda}_{\pi}}{\partial\xi}={}^{c}\pi^{\top}T_{c}\frac{\partial\hat{Q}}{\partial\xi}T_{c}^{\top}\,{}^{c}\pi.$
(46)

If the error $\lambda_{\pi}$ remains constant, then one can naturally ask if
there is a reference frame more convenient than others to calculate the error.
The matrix $Q$ from (9) is composed of the following terms:

$Q=\begin{bmatrix}Q_{p}&q\\\ q^{\top}&N\end{bmatrix}$ (47)

where there is a relation between $q$ and the mean of all points from the
plane, recall observed in the same global reference frame

$q=\sum_{t}^{H}\sum_{n}^{N_{t}}T_{t}\cdot
p_{t,n}=\sum_{t}^{H}q_{t}=N\cdot\mu.$ (48)

The vector $q$ is the summation of $q_{t}$ each the average point in global
coordinates. The $3\times 3$ block $Q_{p}$ is also related to the covariance

$\Sigma_{p}\cdot N=Q_{p}-\frac{1}{N}q\cdot q^{\top}$ (49)

Let the transformation matrix be a centering data transformation by a simple
translation of the mean value:

$T_{c}=\begin{bmatrix}I&-\mu\\\ 0&1\end{bmatrix}$ (50)

In that case, the result of this transformation is

$T_{c}\begin{bmatrix}Q_{p}&q\\\
q^{\top}&N\end{bmatrix}T_{c}^{\top}=\begin{bmatrix}Q_{p}-\frac{1}{N}q\cdot
q^{\top}&0\\\ 0&N\end{bmatrix},$ (51)

where top left block corresponds to the scaled covariance (49).

Alternatively, we could have asked: what is the transformation that applied to
the plane $\pi$ yields a $d$ component zero? Intuitively, the higher $d$, the
more sensible the eigendecomposion is. Then

$\begin{bmatrix}\eta\\\ d\end{bmatrix}=\begin{bmatrix}I&0\\\
-\mu^{\top}&1\end{bmatrix}\begin{bmatrix}\eta\\\ 0\end{bmatrix},$ (52)

which is valid since we already defined $d=-E\\{\eta^{\top}p\\}$ in (5).

We will make use of this form assuming $T_{c}$ is a constant. However, this is
not strictly the case, since $T_{c}(\bm{T})$. In Appendix B we discuss further
how the exact gradient is equal and how the Hessian is slightly different and
dense, however, this component depends on the number of poses. We have decided
to maintain the Hessian matrix diagonal as derived in previous sections.

In the evaluations, we will use the method just described in this section, the
EF-center, which will referred in short as EF.

##  VI Experiments

Experiments are targeted to evaluate performance of different plane-based SLAM
back-end approaches, formulated in the form of graph SLAM. Since the quality
of the SLAM full pipeline could also be affected significantly by detection
and association algorithms on the front-end, the question about quality of the
full pipeline remains out of the scope of this paper and instead we focus only
on the back-end. To achieve this, we consider pre-labeled data from different
sensor modalities by developing our synthetic generator and using labeled data
for ICL NUIM [47] and KITTI [48] sequences from EVOPS benchmark [49].

Figure 1: Evaluation results in synthetic environments with planes swiping
over such parameters as: pose perturbation, number of poses in trajectory,
number of points per plane, point-to-plane noise. _Top_ : RPE translation
error in different simulator configurations, _Down_ : RPE rotation error in
different simulator configurations

The following approaches are taken for evaluation as the main representatives
of methods on graph plane-based SLAM: BAREG [29], $\pi$-Factor [22], BALM [31,
30], EigenFactor (EF, ours). BAREG, $\pi$-Factor and EigenFactor are
implemented in our common framework mrob111https://github.com/prime-slam/mrob
for factor-graph optimization problems. In particular, we use a Levenberg-
Marquadt optimizer with tolerance of $10^{-2}$ for all the methods and the
exact same initial conditions. Since the BALM formulation is also graph SLAM,
we consider the author’s original implementation with default parameters [30].

###  VI-A Complexity Analysis

The theoretical complexity analysis is depicted in Table I. We have included
three categories that are required at each iteration step of the optimization,
assuming that planes are observed over the full sequence.

TABLE I: Complexity comparison of various plane SLAM methods. | # Planes Estim. | Evaluation | Optimization
---|---|---|---
BAREG | $MH$ | $O(HM)$ | $O(H)$
$\pi$-Factor | 0 | $O(HM)$ | $O(H^{3}+M)^{*}$
BALM2 | $M$ | $O(HM)$ | $O(H^{3})$
EF | $M$ | $O(HM)$ | $O(H)$

First column is the number of planes to be estimated. For instance,
$\pi$-Factor does not evaluate the planes, since they are part of their state
variables. Other methods evaluate each plane once and BAREG needs to evaluate
each plane at each pose in the trajectory. Second column — evaluation and it
refers to the complexity to evaluate the residuals for each method and
calculate the Jacobians and Hessian if needed. We see how all methods provide
the same complexity which depends on the number of planes $M$ and poses in the
trajectory $H$.

Finally, the optimization complexity, and this is related to the structure of
the problem. For instance BAREG and EF show linear complexity since they have
a block diagonal matrix to invert. $\pi$-Factor estimates a larger state
vector $H+M$ and requires to invert a sparse matrix (we show worst case) and
Schur trick could be used to marginalize the plane landmarks $M$. BALM2
requires to invert a dense matrix, so it shows cubic complexity.

###  VI-B Synthetic planes

In order to evaluate the quality of planar-based SLAM back-end, we have
developed a generator of synthetic environment with planes.

Setup. The generator provides set of point clouds with planes and reference
trajectory for them, and could be parameterized by number of poses, number of
planes, number of points per plane, point-to-plane noise, trajectory noise.
Example of generated synthetic environment is presented in Fig. 2 The default
parameters configuration is the following: number of poses — 10, number of
planes — 10, number of points per plane — 50, point-to-plane noise — 0.04
meters, trajectory perturbation per each pose — 0.05 meters and 5 degrees. In
evaluation, we chose one parameter to be swiped over different values whereas
other parameters remain constant. Relative Pose Error [50] is used to estimate
the accuracy of methods, rotation and translation part independently.

Figure 2: Synthetic Environment. Points are generated to emulate sampling over
planes.

Results and discussion Results of evaluation are depicted in Fig. 1. In
general, on all parameters variations methods could be ordered as following:
EF and $\pi$-Factor show almost the same quality, BAREG is on the second
place, but still at a relatively small distance, and BALM diverges on huge
perturbations. Since BALM pipeline is used with built-in front-end, the
divergence could be explained by this fact. It is worth to note that on
trajectory perturbation experiment, EF, $\pi$-Factor and BALM show consistent
results on all perturbation ranges, meaning that they show good convergence
and could be applied for global alignment, not only local trajectory
refinement. Swiping over such parameters as number of poses, number of points
per plane, point-to-plane noise the method’s quality is ranged in the same
order, but it could be noticed that in extreme conditions such as huge plane
noise or amount of poses, EF shows more stable results than $\pi$-Factor.

###  VI-C RGBD

Figure 3: _Top_ : Example of the map aggregated from 10 poses with amount of
noise 1 degree, 0.01 meter per pose. _Middle_ : Example of map aggregated from
10 poses with amount of noise 5 degree, 0.05 meter per pose. _Down_ : Result
of alignment using the EF algorithm.

The second stage of evaluation includes quality estimation on RGBD data. Since
evaluation considers only back-ends, we take RGBD sequences from the EVOPS
benchmark [49] with pre-labeled planes. Benchmark has labeled data for ICL
NUIM and TUM RGBD dataset. TUM RGBD has only small pieces of trajectories
covered by enough amount of planes, therefore only ICL NUIM data is used.

Setup For evaluation, we randomly sample 50 different subsequences, whose
length is 30 consecutive frames. For every subsequence, the poses in original
trajectory are perturbed in the range from 0.001-0.15 degrees and 0.001-0.15
meters. Two configuration of ICL NUIM dataset are used: original and one that
emulates Kinect noise. For evaluation, Relative Pose Error is used for
translation and rotation parts independently.

Results and discussion Qualitative example of pose perturbation and later
alignment by EF is depicted in Fig. 3. Since all methods give good map for
visual quality assessment, only one method is chosen for visualization.
Results of evaluation with metrics are presented in Fig. 4. On depth data
without noise, $\pi$-Factor and BAREG show the best results, EF is the second
but still with very small error and BALM shows good results only on small
perturbations, whereas later diverges, that could be explained by interference
with the front-end error. On data that emulates Kinect noise the situation
changes — EF shows the best performance, downgrading one order with respect to
data without noise, the second one result is demonstrated by $\pi$-Factor
which error is two orders more than on original synthetic data, BAREG is the
third one with 3 orders worse quality. The same as in previous synthetic
experiments, EF, $\pi$-Factor and BAREG demonstrate almost constant error
among all ranges of perturbation. To sum up, in this experiment setting EF
shows more stable results when the data is close to real in comparison to
other methods that degrade drastically when more noisy data is provided.

Figure 4: Left: Evaluation on ICL without noise in the point cloud with
respect to perturbations in the initial poses. Right: Evaluation adding noise
to the PC.

###  VI-D LiDAR

Finally, methods are evaluated on KITTI dataset with real LiDAR data.

Setup Original ground truth poses of KITTI dataset are not suitable for
evaluation using classical pose-based metrics, since they provide noisy data
[49]. Therefore, the following scenario is used: ground truth trajectory is
used as the initial conditions for all methods, then the methods optimize over
pre-labeled planes from the EVOPS benchmark. For evaluation we consider 100
submaps of length 30 poses, randomly sampled from map 00. To estimate the
quality of alignment, no-reference map metrics MME [37], MPV [38] and MOM [39]
are used.

Results and discussion No-reference map-metrics statistics over the methods is
presented in Tab. II, BALM is not reported in the results, since it doesn’t
provide any convergence on those sequences. Following results, EF shows
improvement of poses that gives better map quality on all metrics, qualitative
demonstration of this can be found in Fig. 5. Other methods sometimes show
improvement over ground truth poses, but not so drastically as EF. It is worth
to note, that KITTI LiDAR data is noisy and EF shows good results even in
those conditions.

Figure 5: Top: map aggregated from KITTI ground truth poses. Down: map aggregated from poses obtained by EF. TABLE II: No-reference map metrics on KITTI dataset | MME $\downarrow$ | MPV $\downarrow$ | MOM $\downarrow$
---|---|---|---
BAREG | 0.24 | 0.0276 | 0.0046
$\pi$-Factor | 0.26 | 0.0275 | 0.0028
GT | 0.23 | 0.0276 | 0.0029
EF (ours) | 0.17 | 0.0266 | 0.0023

##  VII Conclusions

We have presented the Eigen Factors method for back-end plane SLAM. The main
contribution is on calculating exactly the point error of all points in the
plane at $O(1)$ complexity when evaluating the optimization cost. We have
entirely reformulated the problem, from the perspective of manifold
optimization and we have derived an analytical solution for the gradient and
the Hessian matrix of the EF by using a simple procedure to calculate
derivatives. This derivation turns out to be very efficient for optimization,
since it requires an inversion of a block diagonal matrix (Hessian). The
proposed method is an alternating optimization where first the plane
parameters are calculated in closed form and then a Newton-Raphson
optimization is done on the resulting subproblem by considering the estimated
plane variables fixed.

For the evaluations, we have compared EF with other state-of-the-art methods
in back-end plane SLAM, showing solid results across synthetic and real data
for RGBD and LiDAR. EF presents overall good accuracy in all domains, and
excellent results for different number of poses, number of planes, point noise
and initial conditions. We have also evaluated no-reference metrics for the
quality of the 3D reconstruction, where EF shows superior results.

## Appendix A

###  A-A First Order Derivative over the Eigendecomposition

###### Proof.

By definition, the eigendecomposition is expressed as
$\hat{Q}v=\hat{\lambda}v$. The vector of parameters $v$ is a unit vector s.t.
$||v||^{2}=v^{\top}v=1$. One can exchange the unit vector $v$ for any other
vector $\pi\in\mathbb{P}^{3}$, since $\pi=kv$ and $||\pi||=k^{2}$, from
definition $\hat{\lambda}_{\pi}$ in (11).

$\hat{Q}\pi=\hat{\lambda}\pi$ is derived wrt each of $i$th coordinates of the
variable $\xi_{t,i}$, for simplicity it is written as $\xi_{i}$,

$\displaystyle\hat{Q}\frac{\partial\pi}{\partial\xi_{i}}+\frac{\partial\hat{Q}}{\partial\xi_{i}}\,\pi=\lambda\,\frac{\partial\pi}{\partial\xi_{i}}+\frac{\partial\hat{\lambda}}{\partial\xi_{i}}\,\pi$
(53)
$\displaystyle\text{s.t.}\quad\pi^{\top}\,\frac{\partial\pi}{\partial\xi_{i}}=0.$
(54)

Taking into account that the matrix $\hat{Q}$ is symmetric by construction,
then

$\hat{Q}=\hat{Q}^{\top}\quad\implies\quad\pi^{\top}\hat{Q}=\hat{\lambda}\pi^{\top}.$
(55)

One can pre-multiply the expression (53) by $\pi^{\top}$, substitute (55) and
do some manipulations:

$\displaystyle\pi^{\top}\hat{Q}\,\frac{\partial\pi}{\partial\xi_{i}}+\pi^{\top}\frac{\partial\hat{Q}}{\partial\xi_{i}}\,\pi$
$\displaystyle=\pi^{\top}\hat{\lambda}\,\frac{\partial\pi}{\partial\xi_{i}}+\pi^{\top}\frac{\partial\hat{\lambda}}{\partial\xi_{i}}\,\pi$
$\displaystyle\hat{\lambda}\underbrace{\pi^{\top}\,\frac{\partial\pi}{\partial\xi_{i}}}_{(\ref{eq_diff_svd_v})}+\pi^{\top}\frac{\partial\hat{Q}}{\partial\xi_{i}}\,\pi$
$\displaystyle=\hat{\lambda}\underbrace{\pi^{\top}\frac{\partial\pi}{\partial\xi_{i}}}_{(\ref{eq_diff_svd_v})}+\,\frac{\partial\hat{\lambda}}{\partial\xi_{i}}\underbrace{\pi^{\top}\pi}_{k^{2}}$
$\displaystyle\pi^{\top}\frac{\partial\hat{Q}}{\partial\xi_{i}}\,\pi$
$\displaystyle=\frac{\partial\hat{\lambda}k^{2}}{\partial\xi_{i}}=\frac{\partial\hat{\lambda}_{\pi}}{\partial\xi_{i}}.$
(56)

Recall from (9) that $\hat{Q}=\sum\hat{Q}_{\tau}$, so the only derivative that
is not zero is the one corresponding to the $\hat{Q}_{t}$ component. One can
re-arrange each of the $i$ coordinates in $\xi_{t}$ and finally write

$\frac{\partial\hat{\lambda}_{\pi}}{\partial\xi_{t}}=\pi^{\top}\frac{\partial\hat{Q}_{t}}{\partial\xi_{t}}\,\pi.$
(57)

∎

###  A-B Second Order Derivative over the Eigendecomposition

###### Proof.

The second derivative is obtained after deriving for each component of
$\xi_{t,i}$ and $\xi_{t,j}$, written as $\xi_{i}$ and $\xi_{j}$

$\displaystyle\frac{\partial}{\partial\xi_{j}}\left(\frac{\partial\hat{\lambda}_{\pi}}{\partial\xi_{i}}\right)=\frac{\partial}{\xi_{j}}\left(\pi^{\top}\frac{\partial\hat{Q}_{t}}{\partial\xi_{i}}\pi\right)$
(58)
$\displaystyle\frac{\partial^{2}\hat{\lambda}_{\pi}}{\partial\xi_{j}\partial\xi_{i}}=\frac{\partial\pi^{\top}}{\partial\xi_{j}}\frac{\partial\hat{Q}_{t}}{\partial\xi_{i}}\pi+\pi^{\top}\frac{\partial^{2}\hat{Q}_{t}}{\partial\xi_{j}\partial\xi_{i}}\pi+\pi^{\top}\frac{\partial\hat{Q}_{t}}{\partial\xi_{i}}\frac{\partial\pi}{\partial\xi_{j}}$
(59)

Expression (56) is simplified by the alternating optimization assumption in
(13):

$\frac{\partial^{2}\hat{\lambda}_{\pi}}{\partial\xi_{j}\partial\xi_{i}}=\pi^{\top}\frac{\partial^{2}\hat{Q}_{t}}{\partial\xi_{j}\partial\xi_{i}}\pi.$
(60)

The second condition can be easily verified by noting that this expression
only depends on the matrix $Q_{t}$ at pose $t$ and the contribution from other
poses at $t^{\prime}\neq t$ vanishes after the first differentiation. ∎

##  Appendix B Approximation of the EF Center

The retraction of the centering transformation $T_{c}$ for the particular pose
at $t$ can be written

$\displaystyle R_{T_{c}}(\xi_{t})$
$\displaystyle=\begin{bmatrix}I&-\mu(\xi_{t})\\\
0&1\end{bmatrix}=\begin{bmatrix}I&-\frac{1}{N}\left(\sum_{\tau\neq
t}^{H}q_{\tau}+\operatorname*{\mathrm{Exp}}(\xi_{t})q_{t}\right)\\\
0&1\end{bmatrix}=$
$\displaystyle=\begin{bmatrix}0&\frac{1}{N}(q_{t}-\operatorname*{\mathrm{Exp}}(\xi_{t}){q}_{t})\\\
0&1\end{bmatrix}T_{c}=\mathcal{E}(\xi_{t})T_{c}.$ (61)

where we have used the results from (48). Its derivative at any pose $t$ for
the $i$-th coordinate

$\frac{\partial\mathcal{E}(\xi)}{\partial\xi_{i}}=\begin{bmatrix}0&-\frac{1}{N}G_{i}q_{t}\\\
0&0\end{bmatrix}.$ (62)

Therefore, the total gradient starting from (46)

$\displaystyle\frac{\partial}{\partial\xi_{i}}\Big{(}\mathcal{E}(\xi_{t})T_{c}\operatorname*{\mathrm{Exp}}(\xi_{t})Q\operatorname*{\mathrm{Exp}}(\xi_{t})^{\top}T_{c}^{\top}\mathcal{E}(\xi_{t})^{\top}\Big{)}=$
$\displaystyle\frac{\partial\mathcal{E}(\xi)}{\partial\xi_{i}}T_{c}QT_{c}^{\top}+T_{c}G_{i}Q_{t}T_{c}^{\top}+(\ldots)^{\top}$
(63)

where the first term after the equal is zero (see below) and the second term
was already obtained before in (46). There are 2 more terms as a result of the
differentiation, but they happen to be the transposed of the first two. We
will expand the first term and evaluate its derivative after multiplying by
the plane ${}^{c}\pi$:

$\begin{bmatrix}\eta^{\top}&0\end{bmatrix}\begin{bmatrix}0&-\frac{1}{N}G_{i}q_{t}\\\
0&0\end{bmatrix}\begin{bmatrix}Q_{p}-\frac{1}{N}q\cdot q^{\top}&0\\\
0&N\end{bmatrix}\begin{bmatrix}\eta\\\ 0\end{bmatrix}=0.$ (64)

This expression is always zero $\forall i$, therefore it has null contribution
on the gradient, so the results is exactly the same as obtained in in (46).

The Hessian will be dense since there is a term from $T_{c}$ that depends on
all poses and it does not cancel out as in the gradient case. One could derive
it with the definition of retraction, but this term is weighted roughly by
$\frac{1}{H^{2}}$ so its contribution diminishes fast with the number of
poses.

We have tested empirically that for 2 poses, the relative difference between
the exact Hessian and the approximated Hessian (as a matrix norm) is around
30% and for 20 poses, the difference is less than 1%, so its approximation is
well justified.

## References

  * [1] F. Dellaert and M. Kaess, “Square root SAM: Simultaneous localization and mapping via square root information smoothing,” _The International Journal of Robotics Research_ , vol. 25, no. 12, pp. 1181–1203, 2006.
  * [2] B. Triggs, P. F. McLauchlan, R. I. Hartley, and A. W. Fitzgibbon, “Bundle adjustment—a modern synthesis,” in _International workshop on vision algorithms_ , 1999, pp. 298–372.
  * [3] J. Weingarten and R. Siegwart, “3D SLAM using planar segments,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2006, pp. 3062–3067.
  * [4] A. J. Trevor, J. G. Rogers, and H. I. Christensen, “Planar surface SLAM with 3D and 2D sensors,” in _2012 IEEE International Conference on Robotics and Automation_. IEEE, 2012, pp. 3041–3048.
  * [5] G. Ferrer, “Eigen-Factors: Plane Estimation for Multi-Frame and Time-Continuous Point Cloud Alignment,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2019, pp. 1278–1284.
  * [6] Y. Chen and G. Medioni, “Object modeling by registration of multiple range images,” in _IEEE International Conference on Robotics and Automation_. IEEE, 1991, pp. 2724–2729.
  * [7] P. J. Besl and N. D. McKay, “A method for registration of 3-D shapes,” _IEEE Transactions on pattern analysis and machine intelligence_ , vol. 14, no. 2, pp. 239–256, 1992.
  * [8] Z. Zhang, “Iterative point matching for registration of free-form curves and surfaces,” _International journal of computer vision_ , vol. 13, no. 2, pp. 119–152, 1994.
  * [9] S. Rusinkiewicz and M. Levoy, “Efficient variants of the ICP algorithm,” in _Proceedings Third International Conference on 3-D Digital Imaging and Modeling_. IEEE, 2001, pp. 145–152.
  * [10] B. K. Horn, “Closed-form solution of absolute orientation using unit quaternions,” _Journal of the Optical Society of America A_ , vol. 4, no. 4, pp. 629–642, 1987.
  * [11] K. S. Arun, T. S. Huang, and S. D. Blostein, “Least-squares fitting of two 3-D point sets,” _IEEE Transactions on pattern analysis and machine intelligence_ , no. 5, pp. 698–700, 1987.
  * [12] S. Umeyama, “Least-squares estimation of transformation parameters between two point patterns,” _IEEE Transactions on pattern analysis and machine intelligence_ , vol. 13, no. 4, pp. 376–380, 1991.
  * [13] A. Censi, “An icp variant using a point-to-line metric,” in _IEEE International Conference on Robotics and Automation_ , 2008, pp. 19–25.
  * [14] A. Segal, D. Haehnel, and S. Thrun, “Generalized-ICP,” in _Robotics: science and systems_ , vol. 2, no. 4, 2009.
  * [15] J. Serafin and G. Grisetti, “NICP: Dense normal based point cloud registration,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2015, pp. 742–749.
  * [16] R. Bergevin, M. Soucy, H. Gagnon, and D. Laurendeau, “Towards a general multi-view registration technique,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 18, no. 5, pp. 540–547, 1996.
  * [17] K. Pulli, “Multiview registration for large data sets,” in _Second international conference on 3-d digital imaging and modeling (cat. no. pr00062)_. IEEE, 1999, pp. 160–168.
  * [18] M. Kaess, “Simultaneous localization and mapping with infinite planes,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2015, pp. 4605–4611.
  * [19] P. Geneva, K. Eckenhoff, Y. Yang, and G. Huang, “Lips: Lidar-inertial 3d plane slam,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 123–130.
  * [20] X. Zhang, W. Wang, X. Qi, Z. Liao, and R. Wei, “Point-plane slam using supposed planes for indoor environments,” _Sensors_ , vol. 19, no. 17, 2019\.
  * [21] Y. Taguchi, Y.-D. Jian, S. Ramalingam, and C. Feng, “Point-plane SLAM for hand-held 3D sensors,” in _Robotics and Automation (ICRA), 2013 IEEE International Conference on_. IEEE, 2013, pp. 5182–5189.
  * [22] L. Zhou, D. Koppel, and M. Kaess, “Lidar slam with plane adjustment for indoor environment,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 7073–7080, 2021.
  * [23] J. Zhang and S. Singh, “Loam: Lidar odometry and mapping in real-time.” in _Robotics: Science and Systems_ , vol. 2, no. 9, 2014.
  * [24] ——, “Low-drift and real-time lidar odometry and mapping,” _Autonomous Robots_ , vol. 41, no. 2, pp. 401–416, 2017.
  * [25] H. Ye, Y. Chen, and M. Liu, “Tightly coupled 3d lidar inertial odometry and mapping,” in _International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 3144–3150.
  * [26] T. Whelan, R. F. Salas-Moreno, B. Glocker, A. J. Davison, and S. Leutenegger, “Elasticfusion: Real-time dense slam and light source estimation,” _The International Journal of Robotics Research_ , vol. 35, no. 14, pp. 1697–1716, 2016.
  * [27] J. Behley and C. Stachniss, “Efficient Surfel-Based SLAM using 3D Laser Range Data in Urban Environments,” in _Robotics: Science and Systems_ , 2018.
  * [28] L. Zhou, S. Wang, and M. Kaess, “$\pi$-lsam: Lidar smoothing and mapping with planes,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 5751–5757.
  * [29] H. Huang, Y. Sun, J. Wu, J. Jiao, X. Hu, L. Zheng, L. Wang, and M. Liur, “On bundle adjustment for multiview point cloud registration,” _IEEE Robotics and Automation Letters (RAL)_ , vol. 6, no. 4, pp. 8269–8276, 2021.
  * [30] Z. Liu, X. Liu, and F. Zhang, “Efficient and consistent bundle adjustment on lidar point clouds,” _arXiv preprint arXiv:2209.08854_ , 2022.
  * [31] Z. Liu and F. Zhang, “BALM: Bundle adjustment for lidar mapping,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 3184–3191, 2021\.
  * [32] D. Tian, H. Ochimizu, C. Feng, R. Cohen, and A. Vetro, “Geometric distortion metrics for point cloud compression,” in _IEEE International Conference on Image Processing (ICIP)_ , 2017, pp. 3460–3464.
  * [33] E. Alexiou and T. Ebrahimi, “Point cloud quality assessment metric based on angular similarity,” in _2018 IEEE International Conference on Multimedia and Expo (ICME)_ , 2018, pp. 1–6.
  * [34] E. M. Torlig, E. Alexiou, T. A. Fonseca, R. L. de Queiroz, and T. Ebrahimi, “A novel methodology for quality assessment of voxelized point clouds,” in _Applications of Digital Image Processing XLI_ , A. G. Tescher, Ed., vol. 10752, International Society for Optics and Photonics. SPIE, 2018, pp. 174–190. [Online]. Available: https://doi.org/10.1117/12.2322741
  * [35] G. Meynet, J. Digne, and G. Lavoué, “PC-MSDM: A quality metric for 3d point clouds,” in _Eleventh International Conference on Quality of Multimedia Experience (QoMEX)_ , 2019, pp. 1–3.
  * [36] E. Alexiou and T. Ebrahimi, “Towards a point cloud structural similarity metric,” in _2020 IEEE International Conference on Multimedia Expo Workshops (ICMEW)_ , 2020, pp. 1–6.
  * [37] D. Droeschel, J. Stückler, and S. Behnke, “Local multi-resolution representation for 6d motion estimation and mapping with a continuously rotating 3d laser scanner,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2014, pp. 5221–5226.
  * [38] J. Razlaw, D. Droeschel, D. Holz, and S. Behnke, “Evaluation of registration methods for sparse 3d laser scans,” in _European Conference on Mobile Robots (ECMR)_ , 2015, pp. 1–7.
  * [39] A. Kornilova and G. Ferrer, “Be your own benchmark: No-reference trajectory metric on registered point clouds,” in _European Conference on Mobile Robotics (ECMR)_ , 2021.
  * [40] K. Klasing, D. Althoff, D. Wollherr, and M. Buss, “Comparison of surface normal estimation methods for range sensing applications,” in _IEEE International Conference on Robotics and Automation_ , 2009, pp. 3206–3211.
  * [41] G. Golub and V. Pereyra, “Separable nonlinear least squares: the variable projection method and its applications,” _Inverse problems_ , vol. 19, no. 2, p. R1, 2003.
  * [42] D. Strelow, “General and nested Wiberg minimization: L2 and maximum likelihood,” in _European Conference on Computer Vision (ECCV)_. Springer, 2012, pp. 195–207.
  * [43] P.-A. Absil, R. Mahony, and R. Sepulchre, _Optimization algorithms on matrix manifolds_. Princeton University Press, 2009.
  * [44] J. Nocedal and S. Wright, _Numerical optimization_. Springer Science & Business Media, 2006.
  * [45] K. M. Lynch and F. C. Park, _Modern Robotics_. Cambridge University Press, 2017.
  * [46] J. Sola, J. Deray, and D. Atchuthan, “A micro lie theory for state estimation in robotics,” _arXiv preprint arXiv:1812.01537_ , 2018.
  * [47] A. Handa, T. Whelan, J. McDonald, and A. J. Davison, “A benchmark for rgb-d visual odometry, 3d reconstruction and slam,” in _2014 IEEE international conference on Robotics and automation (ICRA)_. IEEE, 2014, pp. 1524–1531.
  * [48] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The kitti dataset,” _The International Journal of Robotics Research_ , vol. 32, no. 11, pp. 1231–1237, 2013.
  * [49] A. Kornilova, D. Iarosh, D. Kukushkin, N. Goncharov, P. Mokeev, S. A., and G. Ferrer, “EVOPS Benchmark: Evaluation of Plane Segmentation from RGBD and LiDAR Data,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022.
  * [50] R. Kümmerle, B. Steder, C. Dornhege, M. Ruhnke, G. Grisetti, C. Stachniss, and A. Kleiner, “On measuring the accuracy of slam algorithms,” _Autonomous Robots_ , vol. 27, no. 4, pp. 387–407, 2009.