# Protected Fermionic Zero Modes in Periodic Gauge Fields

Võ Tiến Phong<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of
Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104,
U.S.A.  Department of Physics, Florida State University, Tallahassee, FL
32306, U.S.A.  The National High Magnetic Field Laboratory, Tallahassee, FL
32310, U.S.A.  Eugene J. Mele<EMAIL_ADDRESS>Department of Physics
and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, U.S.A.

(August 28, 2024)

###### Abstract

It is well-known that macroscopically-normalizable zero-energy wavefunctions
of spin-$\frac{1}{2}$ particles in a two-dimensional inhomogeneous magnetic
field are spin-polarized and exactly calculable with degeneracy equaling the
number of flux quanta linking the whole system. Extending this argument to
massless Dirac fermions subjected to magnetic fields that have zero net flux
but are doubly periodic in real space, we show that there exist only two
Bloch-normalizable zero-energy eigenstates, one for each spin flavor. This
result is immediately relevant to graphene multilayer systems subjected to
doubly-periodic strain fields, which at low energies, enter the Hamiltonian as
periodic pseudo-gauge vector potentials. Furthermore, we explore various
related settings including nonlinearly-dispersing band structure models and
systems with singly-periodic magnetic fields.

##  I Introduction

The motion of a charged particle in a uniform magnetic field is one of the
simplest well-studied elementary problems at both the classical and quantum
levels [1, 2]. For a spatially-varying field, the problem is more challenging
even in its classical treatment and arises in a variety of physical contexts,
ranging from strategies for trapping of ultracold atoms [3, 4, 5] to the
shapes of orbits for charged particles circulating around magnetic field lines
in plasmas [6, 7]. Surprisingly, the two-dimensional motion of a charge in an
inhomogeneous magnetic field remains an analytically-accessible problem in the
extreme quantum limit owing to a separation of two frequency (energy) scales.
In a strong magnetic field, the kinetic energy is effectively quenched by the
cyclotron motion and the guiding center of the cyclotron orbit can drift
slowly in the presence of additional potentials. Aharonov and Casher famously
found that an exact cancellation of the zero-point energy in the lowest
cyclotron orbit and the Zeeman splitting with gyromagnetic factor $g=2$
produces a spin-polarized zero-energy state that is macroscopically degenerate
[8, 9]. This is the exact analog of the lowest Landau level if the field were
made spatially uniform. Even with spatial variation, the degeneracy is
determined only by the total magnetic flux linking the macroscopic system and
not on the spatial distribution of the field [8, 9, 10].

A variant of this problem arises in two-dimensional materials that are
periodically patterned laterally [11, 12, 13]. For example, in the linearly-
dispersing bands for electrons in a single layer of graphene, a periodic
lattice strain is a momentum boost encoded as an effective pseudo-vector
potential [14]. This is completely analogous to the electromagnetic vector
potential except for its sign change in two time-reversed valleys. If the
strain pattern is made periodic, the total pseudo-flux that links the unit
cell is separately zero for each valley [15, 16, 17, 18, 19, 20, 21, 22]. A
naïve application of the Aharonov-Casher theorem would therefore exclude the
possibility of zero modes since the total flux is zero and there is no analog
to the Zeeman spin polarization energy with $g=2$.

In this article, we exploit the fact that neither of these conditions is
necessary when the pseudo-gauge field varies in space but is made periodic on
a superlattice. The wavefunctions in the periodic problem need only be Bloch
normalizable, with support on a finite real-space supercell instead of a
macroscopic two-dimensional domain. This weakens the normalizability condition
and enables the existence of exactly two zero-energy modes per valley. The
constraint of restricting these modes to a single pseudospin polarization is
thereby also eliminated: the analytic structure of these zero modes protects
one member from both pseudospins (sublattice polarizations) in each valley.
These zero-energy states are parts of dispersive low-energy bands whose
bandwidths can be estimated from the velocity at these zero-energy crossings.
The velocity depends on the strength of the periodic pseudo-field and can be
significantly smaller than the backfolding energy scale produced by the
periodic superlattice. For a general period and a general field strength, one
finds a manifold of spectrally-isolated low-energy bands that possess a
nontrivial quantum geometry [18, 19, 20, 23, 24]. These zero-energy
eigenfunctions have recently been studied in Dirac systems subjected to real
periodic magnetic fields [25, 26] and strain-induced pseudo-magnetic fields
[20]. In the following, we develop these ideas further for the linearly
dispersing Dirac model. We then show that these results generalize to other
nonlinearly-dispersing long-wavelength band structure models.

##  II Aharonov-Casher Argument

For completeness, we begin with a brief summary of Aharonov and Casher’s
construction [8, 9, 27, 28]. We consider a two-dimensional Dirac Hamiltonian
in the presence of a spatially-dependent (not yet assumed to be periodic)
magnetic field
$\mathbf{B}(\mathbf{r})=B(\mathbf{r})\hat{e}_{z}=\left[\partial_{x}A_{y}(\mathbf{r})-\partial_{y}A_{x}(\mathbf{r})\right]\hat{e}_{z}$
of the form

$\mathcal{H}_{1}=\hbar
v_{F}\left(-i\nabla_{\mathbf{r}}+\frac{e}{\hbar}\mathbf{A}(\mathbf{r})\right)\cdot\bm{\sigma},$
(1)

where $v_{F}$ is the Dirac velocity and the Pauli matrices $\bm{\sigma}$ act
on generalized spin space 111One can equivalently consider a Pauli
Hamiltonian, which is $\mathcal{H}_{1}^{2}$ in our notation. The argument
proceeds in the exact same manner for the Pauli Hamiltonian.. Though not
strictly necessary, for simple analytic control, we assume that
$B(\mathbf{r})$ has compact support so that the total flux
$\Phi=\int_{\mathbb{R}^{2}}B(\mathbf{r})d^{2}\mathbf{r}$ is finite:
$\lfloor{|\Phi|/\Phi_{0}\rfloor}=N,$ where $\Phi_{0}=h/e$ is the flux quantum
and $N\in\mathbb{Z}_{>0}$ is a positive integer. Assuming Lorenz gauge
$\nabla\cdot\mathbf{A}=0,$ we can choose a scalar potential $\phi(\mathbf{r})$
such that $\partial_{x}\phi(\mathbf{r})=A_{y}(\mathbf{r})$ and
$\partial_{y}\phi(\mathbf{r})=-A_{x}(\mathbf{r}).$ This scalar potential
satisfies Poisson’s equation sourced by the magnetic field,
$\Delta\phi(\mathbf{r})=B(\mathbf{r}),$ with formal solution

$\phi(\mathbf{r})=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}B(\mathbf{r}^{\prime})\ln|\mathbf{r}-\mathbf{r}^{\prime}|d^{2}\mathbf{r}^{\prime}.$
(2)

By writing the zero-energy eigenstates for Hamiltonian (1) as

$\psi(\mathbf{r})=\begin{pmatrix}e^{+e\phi(\mathbf{r})/\hbar}f_{+}(\mathbf{r})\\\
e^{-e\phi(\mathbf{r})/\hbar}f_{-}(\mathbf{r})\end{pmatrix},$ (3)

we find that $\left(\partial_{x}\pm
i\partial_{y}\right)f_{\pm}\left(\mathbf{r}\right)=0.$ This implies that
$f_{\pm}$ are both entire functions 222More precisely,
$\partial_{\bar{z}}f_{+}=0$ and $\partial_{z}f_{-}=0.$ So $f_{+}$ is
holomorphic with respect to $z=x+iy$ and $f_{-}$ is holomorphic with respect
to $\bar{z}=x-iy$, which is usually called antiholomorphic in complex
analysis. Therefore, both functions have power-series representations
$\begin{split}f_{+}(z)=\sum_{n=0}^{\infty}a_{n}z^{n}\quad\text{and}\quad
f_{-}(\bar{z})=\sum_{n=0}^{\infty}b_{n}\bar{z}^{n},\end{split}$ (4) where
$a_{n}\in\mathbb{C}$ and $b_{n}\in\mathbb{C}$ are constants. The forms of
these functions are constrained by the normalizability of the wavefunction. To
ascertain these constraints, we observe that as
$|\mathbf{r}|\rightarrow\infty,$ the scalar potential tends to
$\phi(\mathbf{r})\rightarrow\Phi\ln\left|\mathbf{r}\right|/2\pi=\ln\left|\mathbf{r}\right|^{\Phi/2\pi}.$
The exponentials in Eq. (3) have asymptotic behaviors $e^{\pm
e\phi(\mathbf{r})/\hbar}\rightarrow|\mathbf{r}|^{\pm\Phi/\Phi_{0}}.$ Since
entire functions do not decay globally, we only admit $e^{-\eta
e\phi(\mathbf{r})/\hbar}$ to ensure normalizability, where
$\eta=\text{sign}\left(\Phi\right)$. Finally, we require
$\lim_{|\mathbf{r}|\rightarrow\infty}|\mathbf{r}||f_{-\eta}(\mathbf{r})|e^{-\eta
e\phi(\mathbf{r})/\hbar}=0,$ which implies that $f_{+}(z=x+iy)$ or
$f_{-}(\bar{z}=x-iy)$ is a polynomial of degree at most $N-1.$ Thus, the $N$
independent solutions are

$\psi_{n}(\mathbf{r})=\begin{pmatrix}\Theta\left[-\eta\right]\\\
\Theta\left[+\eta\right]\end{pmatrix}e^{-\eta e\phi(\mathbf{r})/\hbar}(x-i\eta
y)^{n},$ (5)

for $n=0,1,2,...,N-1.$ Here, $\Theta\left[\eta\right]$ is the Heaviside theta
function. In brief, Aharonov and Casher showed that electrons in a magnetic
field with total flux $\lfloor|\Phi|\rfloor=N\Phi_{0}$ have $N$ zero-energy
eigenstates that are spin polarized. In particular, the wavefunctions can be
written as products of analytic functions of $z$ or $\bar{z}$ and exponentials
of the scalar function $\phi(\mathbf{r})$. This analysis is quite general
since it does not assume any particular form of the magnetic field except that
it is compactly supported. So the degeneracy of the zero modes is not mandated
by a spatial symmetry. However, chiral symmetry
$\sigma_{z}\mathcal{H}\sigma_{z}=-\mathcal{H}$ is crucial to the existence of
these modes as mass terms $m\sigma_{z}$ would lift them away from zero energy.

##  III Zero Modes in a Doubly-Periodic Magnetic Field with No Net Magnetic
Flux

The above analysis suggests that a magnetic field with zero flux cannot induce
any zero mode. This is only true in the space of normalizable wavefunctions in
the entire plane for which the preceding analysis applies [9]. However, in
certain situations, the relevant domain for the wavefunctions is not the
entire plane. The spectrum of Dirac and Pauli operators in the presence of a
magnetic field have been studied in various different domains and with fields
of different regularities [31]. Of particular interest to us are systems with
a magnetic field that is periodic in two independent directions with primitive
lattice vectors $\mathbf{L}_{1}$ and $\mathbf{L}_{2}:$
$B\left(\mathbf{r}+n_{1}\mathbf{L}_{1}+n_{2}\mathbf{L}_{2}\right)=B\left(\mathbf{r}\right),$
where $n_{1},n_{2}$ are integers. We can write the magnetic field as a Fourier
series

$B(\mathbf{r})=\sum_{\mathbf{G}}\tilde{B}_{\mathbf{G}}e^{i\mathbf{G}\cdot\mathbf{r}}.$
(6)

For physically-relevant fields, it is often enough to approximate the magnetic
field with a finite number of Fourier harmonics. Therefore, we can assume that
the magnetic field is defined by its finite Fourier series 333Though this
assumption may seem strict, it is often employed in numerical studies due to
the impracticality of including an infinite number of wavevectors. This
restriction can be relaxed considerably, but such a generalization is of
secondary importance for us here. We focus on the case of zero magnetic flux,
$\tilde{B}_{\mathbf{0}}=0.$ Because of that, the vector potential can be
written explicitly as

$\begin{split}A_{x}(\mathbf{r})&=\sum_{\mathbf{G}\neq\mathbf{0}}\frac{iG_{y}}{|\mathbf{G}|^{2}}\tilde{B}_{\mathbf{G}}e^{i\mathbf{G}\cdot\mathbf{r}},\\\
A_{y}(\mathbf{r})&=\sum_{\mathbf{G}\neq\mathbf{0}}-\frac{iG_{x}}{|\mathbf{G}|^{2}}\tilde{B}_{\mathbf{G}}e^{i\mathbf{G}\cdot\mathbf{r}}.\end{split}$
(7)

Thus, we conclude that $\mathbf{A}(\mathbf{r})$ has the same periodicity as
$B(\mathbf{r}).$ It is worth emphasizing that this fact follows from the
vanishing magnetic flux. If $\tilde{B}_{\mathbf{0}}$ were not zero, there
would have been a non-periodic component to the vector potential, such as
$A_{y}(\mathbf{r})=\tilde{B}_{\mathbf{0}}x+...$ This is exactly like a
translationally-invariant constant magnetic field having magnetic vector
potentials that are not translationally-invariant. With
$\mathbf{A}(\mathbf{r})$ proven to be periodic when $\mathbf{B}(\mathbf{r})$
carries no flux, it is clear that the Hamiltonian (1) is spatially periodic.
It thus follows from Bloch’s theorem that the eigenfunctions must have
periodic norm. The appropriate solution space now becomes that of normalizable
wavefunctions on a compact torus, not the entire plane.

We now adapt the argument of Ref. [8] to the restricted setting of a torus.
Still writing the zero-energy eigenstates as in Eq. (3), where
$\phi(\mathbf{r})=-\sum_{\mathbf{G}\neq\mathbf{0}}\tilde{B}_{\mathbf{G}}e^{i\mathbf{G}\cdot\mathbf{r}}/|\mathbf{G}|^{2}$
is also a continuous periodic function, we still find that $f_{+}(z)$ and
$f_{-}(\bar{z})$ are analytic functions of $z$ and $\bar{z}$ respectively.
Now, by requiring that the wavefunction be continuous and have periodic norm,
it follows that the wavefunction components $e^{\pm
e\phi(\mathbf{r})/\hbar}f_{\pm}(\mathbf{r})$ must be globally bounded as well.
Since the exponential factors are globally bounded as $\phi(\mathbf{r})$ is
periodic, it must be the case that $|f_{\pm}(z)|$ is globally bounded. By
Liouville’s theorem that globally-bounded, entire functions are constants, we
find that $f_{\pm}$ must be spatially uniform. Thus, there are only two
independent zero-mode solutions with periodic norm, henceforth called Bloch
zero modes:

$\psi_{+}(\mathbf{r})=\frac{1}{A_{+}}\begin{pmatrix}e^{e\phi(\mathbf{r})/\hbar}\\\
0\end{pmatrix}\text{ and
}\psi_{-}(\mathbf{r})=\frac{1}{A_{-}}\begin{pmatrix}0\\\
e^{-e\phi(\mathbf{r})/\hbar}\end{pmatrix},$ (8)

where $A_{\pm}$ are normalization constants are given by

$A_{\pm}=\left[\int_{\Omega}e^{\pm
2e\phi(\mathbf{r})/\hbar}d^{2}\mathbf{r}\right]^{\frac{1}{2}},$ (9)

and $\Omega$ is the unit cell. These same solutions were also studied in Ref.
[20]. Eq. (8) shows that even when the net magnetic flux is zero, there are
still two zero modes, but these modes are only normalizable within a unit
cell. Furthermore, on a torus, we have zero modes of both spin flavors,
contrary to the original formulation where the zero modes are spin polarized.
However, these two Bloch zero modes feature spatial spin isolation because
$\psi_{+}$ is enhanced precisely where $\psi_{-}$ is suppressed due to the
different signs in the exponential. Examples are shown in Fig. 1.

These Bloch zero modes can be interpreted in the context of a band structure.
In the absence a magnetic field, the energies of Hamiltoninan (1) form two
linear branches: $\mathcal{E}_{\pm}=\pm\hbar v_{F}|\mathbf{k}|,$ where
$\mathbf{k}$ is the wavevector. These branches cross exactly at
$\mathbf{k}=\mathbf{0}.$ In the presence of a periodic magnetic field, the
spectrum consists of bands defined within a Brillouin zone. Generically, one
would expect the degeneracy point at $\mathcal{E}=0$ to be gapped out by a
general periodic field without any symmetry. Our analysis proves the contrary,
that the degeneracy point remains intact no matter the form of the magnetic
field. Therefore, the bands near $\mathcal{E}=0$ must at minimum form a
doublet set. The Dirac velocity is, however, renormalized by the magnetic
field. Using first-order perturbation theory, the renormalized velocity
$v_{\text{renorm}}$ is given by

$\frac{v_{\text{renorm}}}{v_{F}}=\frac{|\Omega|}{A_{+}A_{-}}.$ (10)

By the Cauchy-Schwarz inequality, $|\Omega|\leq A_{+}A_{-}.$ So the velocity
is always renormalized downward as expected. In order for the velocity to
vanish, $A_{\pm}\rightarrow\infty.$ However, as long as $e^{\pm
e\phi(\mathbf{r})/\hbar}$ is integrable, which we assume, this condition is
never exactly satisfied. So, for physical magnetic fields, the bands can be
made very narrow, but never exactly flat, at least to first order in
perturbation theory.

Figure 1: Zero modes in a periodic magnetic field. The magnetic field for the
top panels forms a triangular lattice
$B(\mathbf{r})=B_{0}\sum_{i=1}^{3}\cos\left(\mathbf{G}\cdot\mathbf{r}\right),$
where $\mathbf{G}_{1}=G\left(0,1\right),$
$\mathbf{G}_{2}=G\left(-\sqrt{3}/2,-1/2\right),$ and
$\mathbf{G}_{1}=G\left(\sqrt{3}/2,-1/2\right).$ The magnetic field for the
bottom panels is $B(\mathbf{r})=B_{0}\left(\cos 2\pi x+\cos 2\pi y\right).$ We
emphasize that the two eigenstates are localized on dual lattices.

It is straightforward to show the existence of Bloch zero modes in a variety
of other settings. To start, let us consider a $4\times 4$ Hamiltonian
inspired by Bernal bilayer graphene of the form

$\mathcal{H}_{2}=\hbar v_{F}\begin{pmatrix}0&\Pi_{-}&0&0\\\
\Pi_{+}&0&\gamma_{1}/\hbar v_{F}&0\\\ 0&\gamma_{1}/\hbar v_{F}&0&\Pi_{-}\\\
0&0&\Pi_{+}&0\end{pmatrix},$ (11)

where $\Pi_{\pm}=\left(-i\partial_{x}+\frac{e}{\hbar}A_{x}\right)\pm
i\left(-i\partial_{y}+\frac{e}{\hbar}A_{y}\right)$ and $\gamma_{1}$ is a
constant. Again, we assume that $B(\mathbf{r})$ is periodic and
$\phi(\mathbf{r})$ is a scalar potential defined as before. Then, we write the
zero-energy eigenstates for Hamiltonian (11) as

$\psi(\mathbf{r})=\begin{pmatrix}e^{+e\phi(\mathbf{r})/\hbar}f_{1,+}(\mathbf{r})\\\
e^{-e\phi(\mathbf{r})/\hbar}f_{1,-}(\mathbf{r})\\\
e^{+e\phi(\mathbf{r})/\hbar}f_{2,+}(\mathbf{r})\\\
e^{-e\phi(\mathbf{r})/\hbar}f_{2,-}(\mathbf{r})\\\ \end{pmatrix}.$ (12)

We find the following conditions: $\partial_{z}f_{1,-}=0$ and
$\partial_{\bar{z}}f_{2,+}=0,$ which imply that $f_{1,-}$ and $f_{2,+}$ must
be constants. If $f_{1,-}\neq 0$ and $f_{2,+}\neq 0$, then the remaining two
functions satisfy, for complex constants $c_{1}=\frac{\gamma_{1}}{2i\hbar
v_{F}}f_{2,+}$ and $c_{2}=\frac{\gamma_{1}}{2i\hbar v_{F}}f_{1,-},$
$\partial_{\bar{z}}f_{1,+}=c_{1}$ and $\partial_{z}f_{2,-}=c_{2}.$ This
implies that there are functions $\mathcal{F}_{1,+}=f_{1,+}-c_{1}\bar{z}$ and
$\mathcal{F}_{2,-}=f_{2,-}-c_{2}z$ satisfying
$\partial_{\bar{z}}\mathcal{F}_{1,+}=0$ and $\partial_{z}\mathcal{F}_{2,-}=0.$
So, $\mathcal{F}_{1,+}$ and $\mathcal{F}_{2,-}$ are holomorphic with respect
to $z$ and $\bar{z}$ respectively. Consequently, the original functions can be
written as

$\begin{split}f_{1,+}=c_{1}\bar{z}+\mathcal{F}_{1,+},\quad
f_{2,-}=c_{2}z+\mathcal{F}_{2,-}.\\\ \end{split}$ (13)

Now imposing a global bound $B_{1},$ we observe via the reverse triangle
inequality that for $f_{1,+}$

$|\mathcal{F}_{1,+}|-|c_{1}\bar{z}|\leq|c_{1}\bar{z}+\mathcal{F}_{1,+}|<B_{1},$
(14)

which implies that $\mathcal{F}_{1,+}$ has at most linear growth,
$|\mathcal{F}_{1,+}|<B_{1}+|c_{1}||z|.$ Now, because $\mathcal{F}_{1,+}$ is
entire, the generalized Liouville’s theorem states that we can write
$\mathcal{F}_{1,+}=a_{0}+a_{1}z.$ A similar reasoning applies to
$\mathcal{F}_{2,-}.$ In the end, we can in general write

$\begin{split}f_{1,+}=a_{0}+a_{1}z+c_{1}\bar{z},\quad
f_{2,-}=b_{0}+b_{1}\bar{z}+c_{2}z.\\\ \end{split}$ (15)

Finally, imposing periodicity on $|f_{1,+}|$ and $|f_{2,-}|$ eliminates
$a_{1},c_{1},b_{1},c_{2}.$ Thus, we arrive at the conclusion that there are
only two Bloch zero modes, which can be written explicitly as

$\psi_{+}(\mathbf{r})=\frac{1}{A_{+}}\begin{pmatrix}e^{+e\phi(\mathbf{r})/\hbar}\\\
0\\\ 0\\\ 0\\\ \end{pmatrix}\text{ and
}\psi_{-}(\mathbf{r})=\frac{1}{A_{-}}\begin{pmatrix}0\\\ 0\\\ 0\\\
e^{-e\phi(\mathbf{r})/\hbar}\\\ \end{pmatrix}.$ (16)

It is worth pointing out the formal similarity between Eq. (16) and Eq. (8).
The only difference between the two is the number of internal degrees of
freedom. The band structure of Bernal bilayer graphene under a periodic
pseudo-magnetic field was studied in Ref. [23], wherein these same zero-modes
were found.

The argument above can be extended to show that for any number of layers for
chirally-stacked multilayer, there are also two zero modes given explicitly by
formulas similar to Eq. (16). The argument is simple but tedious; essentially,
it is just a recursion of the steps done in the bilayer graphene. So we leave
it for Appendix A. Here, we present an alternative, much quicker, method to
obtain the same result. In the absence of a magnetic field, it is well-known
that the low-energy spectrum of chirally-stacked multilayer graphene is a
polynomial two-band crossing of the form
$\mathcal{E}_{\pm}\propto|\mathbf{k}|^{\ell},$ where $\ell$ is the number of
layers [33]. The appropriate Hamiltonian describing only these two bands is
also of chiral form

$\mathcal{H}_{\text{eff},\ell}\propto\begin{pmatrix}0&\Pi_{-}^{\ell}\\\
\Pi_{+}^{\ell}&0\end{pmatrix}.$ (17)

There are two Bloch zero modes to this Hamiltonian (17) given by Eq. (8). This
follows immediately by noting the properties
$\Pi_{-}e^{-e\phi(\mathbf{r})/\hbar}f_{-}=-2ie^{-e\phi(\mathbf{r})/\hbar}\partial_{z}f_{-}$
and
$\Pi_{+}e^{+e\phi(\mathbf{r})/\hbar}f_{+}=-2ie^{+e\phi(\mathbf{r})/\hbar}\partial_{\bar{z}}f_{+}.$
So the exponentials can be pulled past the derivatives, which then annihilate
the remaining constant, nulling the whole function as desired. Therefore, the
two Bloch zero modes in Eq. (8) satisfy band degeneracy of any order, not just
linear band crossings. As presented, this method does not exclude the
possibility that there may be more than two zero modes, but with some more
work, one can probably eliminate that possibility as well.

It is worth mentioning briefly that the aforementioned considerations
immediately imply that odd-layer $ABA$ multilayer graphene in a periodic gauge
field with zero flux also possesses zero modes. This is because odd-layer
$ABA$ multilayer graphene, due to a layer-exchange symmetry, can be decomposed
into direct sum of chiral sectors which host zero modes [34]. As an example,
we show this for $ABA$ trilayer graphene. Then, the generalization to any
number of layers should be straightforward. The Hamiltonian for an $ABA$
trilayer is

$\mathcal{H}_{3}=\hbar v_{F}\begin{pmatrix}0&\Pi_{-}&0&0&0&0\\\
\Pi_{+}&0&\gamma_{1}/\hbar v_{F}&0&0&0\\\ 0&\gamma_{1}/\hbar
v_{F}&0&\Pi_{-}&0&\gamma_{1}/\hbar v_{F}\\\ 0&0&\Pi_{+}&0&0&0\\\
0&0&0&0&0&\Pi_{-}\\\ 0&0&\gamma_{1}/\hbar v_{F}&0&\Pi_{+}&0\\\ \end{pmatrix}.$
(18)

Upon a unitary transformation, this Hamiltonian can be brought into a direct
sum of a monolayer and a bilayer

$\tilde{\mathcal{H}_{3}}=\hbar v_{F}\begin{pmatrix}0&\Pi_{-}&0&0&0&0\\\
\Pi_{+}&0&0&0&0&0\\\ 0&0&0&\Pi_{-}&0&0\\\
0&0&\Pi_{+}&0&\sqrt{2}\gamma_{1}/\hbar v_{F}&0\\\
0&0&0&\sqrt{2}\gamma_{1}/\hbar v_{F}&0&\Pi_{-}\\\ 0&0&0&0&\Pi_{+}&0\\\
\end{pmatrix}.$ (19)

So the zero modes analyzed before remain valid in this situation. In this
particular example, the count of zero modes is four: two for the monolayer
sector and two for the bilayer sector.

##  IV Zero Modes in a Singly-Periodic Magnetic Field

For a final generalization, we consider a magnetic field that is periodic
along one direction only. A different, but closely related, problem was
studied in Ref. [9] using doubly-periodic Weierstrass sigma functions. Here,
we do not assume that the magnetic field is also periodic in the second
direction. Without loss of generality, let $B(x,y)=B(x+n,y),$ where $n$ is an
integer, repeat in the $x$-direction but be a general function in the
$y$-direction. Other cases can be similarly obtained via a rotation and
scaling of coordinates. The dimensions of $x,y$ are suppressed. We do not
assume that this magnetic field has vanishing flux. We write the scalar
potential as [35]

$\begin{split}\phi(x,y)&=\int_{\mathbb{R}}dy^{\prime}\int_{-\frac{1}{2}}^{\frac{1}{2}}dx^{\prime}B(x^{\prime},y^{\prime})G(x-x^{\prime},y-y^{\prime}),\\\
G(x,y)&=\frac{1}{4\pi}\ln\left[\cosh\left(2\pi y\right)-\cos\left(2\pi
x\right)\right].\end{split}$ (20)

A brief derivation of the Green’s function is presented in Appendix B. It is
straightforward to check that $\phi(x,y)=\phi(x+n,y)$ as desired. It
immediately follows that the corresponding $\mathbf{A}(\mathbf{r})$ is also
periodic in the $x$-direction. As a consequence, the Hamiltonian is defined on
a cylinder. In the limit $|y|\rightarrow\infty,$
$\phi(\mathbf{r})\rightarrow\Phi\ln\left[\cosh\left(2\pi|y|\right)\right]/4\pi,$
where $\Phi=\int_{\mathbb{R}}dy\int_{-\frac{1}{2}}^{\frac{1}{2}}dxB(x,y).$ So
the exponential factors tend to the following limits: $e^{\pm
e\phi/\hbar}\rightarrow\left[\cosh\left(2\pi|y|\right)\right]^{\pm\Phi/2\Phi_{0}}.$
For large $|y|,$ $\cosh(2\pi|y|)\rightarrow e^{2\pi|y|}/2.$ So, we have
$e^{\pm e\phi/\hbar}\rightarrow e^{\pm\pi|y|\Phi/\Phi_{0}}.$ Because we insist
the $f_{\pm}$ functions have periodic norm in the $x$-direction, they must
grow (unless they are constants) in the $y$-direction. Therefore, we must
again choose $e^{-e\eta\phi(\mathbf{r})/\hbar}.$ We initially consider
$f_{\pm}$ functions with period $1$ in the $x$-direction: $f_{+}(z)=e^{2\pi
miz}$ and $f_{-}(\bar{z})=e^{-2\pi mi\bar{z}},$ where $m$ are integers. We
need to determine restrictions on $m$ to ensure normalizability. If $\Phi>0,$
then we have $e^{-\phi(\mathbf{r})/\hbar}f_{-}(\bar{z})\rightarrow e^{-2\pi
mix-\left(2\pi my+\pi|y|\Phi/\Phi_{0}\right)}\rightarrow 0$ for both positive
and negative large $y$ if $|m|<\Phi/2\Phi_{0}.$ If $\Phi<0,$ then we have
$e^{+\phi(\mathbf{r})/\hbar}f_{+}(z)\rightarrow e^{+2\pi mix-\left(2\pi
my-\pi|y|\Phi/\Phi_{0}\right)}\rightarrow 0$ for both positive and negative
large $y$ if $|m|<-\Phi/2\Phi_{0}.$ If we lift the requirement that $f_{\pm}$
has period $1$ in $x$ but instead has integer period $M>1,$ then everything
stays the same in the above analysis except for the replacement $m\rightarrow
m/M.$ This extension in the period is allowed by Bloch’s theorem because Bloch
eigenstates do not need to be periodic, only their norms need to be. We can
write this new $m$ as $m=qM+p,$ where $q\in\mathbb{Z}$ and
$p\in\left[0,M-1\right].$ Using this, we can write eigenstates explicitly in
Bloch form $\psi_{k_{x},q}(\mathbf{r})=e^{-i\eta
k_{x}x}u_{k_{x},q}(\mathbf{r}),$ where

$u_{k_{x},q}(\mathbf{r})=\begin{pmatrix}\Theta\left[-\eta\right]\\\
\Theta\left[+\eta\right]\end{pmatrix}e^{-\eta e\phi(\mathbf{r})/\hbar}e^{-2\pi
q\eta ix-\left(2\pi q+k_{x}\right)y},$ (21)

where $k_{x}=2\pi p/M=\left[0,2\pi\right)$ and
$u_{k_{x},q}(x,y)=u_{k_{x},q}(x+n,y).$ In the limit $M\rightarrow\infty,$
$k_{x}$ becomes a continuous variable. The indices are still subject to the
constraint $|q+k_{x}/2\pi|<|\Phi|/2\Phi_{0}.$

As an example, we take $B(x,y)=B_{0}+B_{1}\cos\left(2\pi x/L\right),$ where
$L$ is the period [36]. Strictly speaking, the preceding analysis does not
apply to this magnetic profile because it is not compactly supported and the
magnetic flux can be infinite. However, as we will show, Eq. (21) still
produces the correct zero-energy solutions. The corresponding scalar potential
is

$\phi(x,y)=\frac{B_{0}y^{2}}{2}-\frac{B_{1}L^{2}}{4\pi^{2}}\cos\left(\frac{2\pi}{L}x\right).$
(22)

Assuming $B_{0}>0,$ the eigenstates are

$\psi_{k_{x}}(\mathbf{r})=\frac{e^{-\frac{eB_{0}y^{2}}{2\hbar}+\frac{eB_{1}L^{2}}{4\pi^{2}\hbar}\cos\left(\frac{2\pi}{L}x\right)-ik_{x}x-k_{x}y}}{A_{k_{x}}}\begin{pmatrix}0\\\
1\end{pmatrix},$ (23)

where the normalization is given explicitly as

$A_{k_{x}}=\left(\frac{\pi
L^{2}\hbar}{B_{0}e}\right)^{\frac{1}{4}}\exp\left(\frac{k_{x}^{2}\hbar}{2eB_{0}}\right)\left(I_{0}\left[\frac{eB_{1}L^{2}}{2\pi^{2}\hbar}\right]\right)^{\frac{1}{2}},$
(24)

where $I_{0}(x)$ is the modified Bessel function. From the normalization
factor, we see that the eigenstates exist for all real $k_{x};$ this is
because the Gaussian factor $e^{-y^{2}}$ decays much faster than the remaining
factors. If $B_{1}=0,$ we recover exactly the lowest Landau level. This
example illustrates a very general procedure. If the magnetic field has a
constant component, then one can choose an axis along which the magnetic
vector potential is non-periodic. Then, the remaining orthogonal direction can
be periodic or not. If it is, then the final generalization just discussed can
be applied, as was done in Ref. [25]. If it is not, then we are back to the
original Aharonov-Casher setup.

##  V Discussion and Conclusion

In closing, it is worth emphasizing again that our arguments apply only to
fermions described by the Dirac equation (or in the multilayer graphene case,
described by the $2\ell\times 2\ell$ Hamiltonians, where $\ell$ is the number
of layers). The complementary problem concerning the spectrum of Schrodinger
fermions in a periodic magnetic field governed by the Hamiltonian
$\mathcal{H}\propto(\mathbf{p}+e\mathbf{A})^{2}$ has a very different
structure and is considerably more complicated [37, 38, 39, 40, 41, 42, 43,
44]. Whereas Schrodinger fermions are historically relevant to physics in two-
dimensional electron gases, Dirac fermions are prominent in modern two-
dimensional materials [45, 46, 47]. Graphene is probably the most well-known
member of this family. In addition to a real magnetic field, the same physics
can be obtained in graphene by subjecting it to a strain field since such a
field behaves effectively as a pseudo-magnetic vector potential necessarily
with zero flux. Therefore, our analysis is especially relevant to graphene and
its multilayer cousins. Beyond graphene, Dirac fermions can also be found at
boundaries of topological insulators [48, 25]. Though these boundary spectra
generically disperse linearly, they can have nonlinear dispersions as well
like in topological crystalline insulators [49]. In these cases, subjecting
the surfaces with topologically-nontrivial boundary states to a patterned
periodic magnetic field should induce the indicated manifolds of zero modes
that are localized on the boundaries.

We acknowledge funding from the U.S. Department of Energy under grant DE-
FG02-84ER45118.

##  Appendix A Zero Modes in $\ell$-Layer Chirally-Stacked Multilayer Graphene

In this section, we prove that there are two zero-modes for $\ell$-layer
chirally-stacked multilayer graphene in doubly-periodic magnetic fields. This
is analogous to the method presented for $\ell=2$ in the main text. The
Hamiltonian is

$\mathcal{H}_{\ell}=\hbar v_{F}\begin{pmatrix}0&\Pi_{-}&0&0&\ldots&0&0&0&0\\\
\Pi_{+}&0&\gamma_{1}/\hbar v_{F}&0&\ldots&0&0&0&0\\\ 0&\gamma_{1}/\hbar
v_{F}&0&\Pi_{-}&\ldots&0&0&0&0\\\ 0&0&\Pi_{+}&0&\ldots&0&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&0&\ldots&0&\Pi_{-}&0&0\\\ 0&0&0&0&\ldots&\Pi_{+}&0&\gamma_{1}/\hbar
v_{F}&0\\\ 0&0&0&0&\ldots&0&\gamma_{1}/\hbar v_{F}&0&\Pi_{-}\\\
0&0&0&0&\ldots&0&0&\Pi_{+}&0\\\ \end{pmatrix}_{2\ell\times 2\ell},$ (25)

where $\Pi_{\pm}=\left(-i\partial_{x}+\frac{e}{\hbar}A_{x}\right)\pm
i\left(-i\partial_{y}+\frac{e}{\hbar}A_{y}\right)=-i\left(\partial_{x}\mp\frac{e}{\hbar}A_{y}\right)+\left(\pm\partial_{y}+\frac{e}{\hbar}A_{x}\right)=-i\left(\partial_{x}\mp\frac{e}{\hbar}\partial_{x}\phi\right)+\left(\pm\partial_{y}-\frac{e}{\hbar}\partial_{y}\phi\right).$
It is clear that for any function of the form
$e^{+e\phi(\mathbf{r})/\hbar}f_{+}(\mathbf{r})$ or
$e^{-e\phi(\mathbf{r})/\hbar}f_{-}(\mathbf{r}),$ we have

$\begin{split}\Pi_{+}\left[e^{+e\phi(\mathbf{r})/\hbar}f_{+}(\mathbf{r})\right]&=e^{+e\phi(\mathbf{r})/\hbar}\left(-i\partial_{x}+\partial_{y}\right)f_{+}(\mathbf{r})=-ie^{+e\phi(\mathbf{r})/\hbar}\left(\partial_{x}+i\partial_{y}\right)f_{+}(\mathbf{r})=-2ie^{+e\phi(\mathbf{r})/\hbar}\partial_{\bar{z}}f_{+}(\mathbf{r}),\\\
\Pi_{-}\left[e^{-e\phi(\mathbf{r})/\hbar}f_{-}(\mathbf{r})\right]&=e^{-e\phi(\mathbf{r})/\hbar}\left(-i\partial_{x}-\partial_{y}\right)f_{-}(\mathbf{r})=-ie^{-e\phi(\mathbf{r})/\hbar}\left(\partial_{x}-i\partial_{y}\right)f_{-}(\mathbf{r})=-2ie^{-e\phi(\mathbf{r})/\hbar}\partial_{z}f_{-}(\mathbf{r}).\end{split}$
(26)

Then, writing the zero-mode eigenfunctions as

$\psi(\mathbf{r})=\begin{pmatrix}e^{+e\phi(\mathbf{r})/\hbar}f_{1,+}&e^{-e\phi(\mathbf{r})/\hbar}f_{1,-}&e^{+e\phi(\mathbf{r})/\hbar}f_{2,+}&e^{-e\phi(\mathbf{r})/\hbar}f_{2,-}&\ldots&e^{+e\phi(\mathbf{r})/\hbar}f_{\ell,+}&e^{-e\phi(\mathbf{r})/\hbar}f_{\ell,-}\end{pmatrix}^{T},$
(27)

we are led to the following conditions:

$\begin{split}\partial_{z}f_{1,-}&=0,\\\ \gamma_{1}f_{1,-}-2i\hbar
v_{F}\partial_{z}f_{2,-}&=0,\\\ \gamma_{1}f_{2,-}-2i\hbar
v_{F}\partial_{z}f_{3,-}&=0,\\\ &\vdots\\\ \gamma_{1}f_{\ell-1,-}-2i\hbar
v_{F}\partial_{z}f_{\ell,-}&=0,\\\ \end{split}$ (28)

and

$\begin{split}-2i\hbar v_{F}\partial_{\bar{z}}f_{1,+}+\gamma_{1}f_{2,+}&=0,\\\
-2i\hbar v_{F}\partial_{\bar{z}}f_{2,+}+\gamma_{1}f_{3,+}&=0,\\\ &\vdots\\\
-2i\hbar v_{F}\partial_{\bar{z}}f_{\ell-1,+}+\gamma_{1}f_{\ell,+}&=0,\\\
\partial_{\bar{z}}f_{\ell,+}&=0,\\\ \end{split}$ (29)

A symmetry is clear: the $-$ series does not couple to the $+$ series. This is
essential to the argument. Now, let us focus on the $+$ series. The last
condition requires that $f_{\ell,+}$ be an entire function. Since it is
bounded, it must be a constant. Then, we have
$\partial_{\bar{z}}f_{\ell-1,+}=\frac{\gamma_{1}}{2i\hbar
v_{F}}f_{\ell,+}=c_{\ell,+}.$ This means that we can write
$f_{\ell-1,+}=\mathcal{F}_{\ell-1,+}+c_{\ell,+}\bar{z},$ where
$\mathcal{F}_{\ell-1,+}$ is holomorphic. Now, since $f_{\ell-1,+}$ has
periodic norm, we can write its bound as $B_{\ell-1,+}.$ Then, by using the
reverse triangle inequality, we have

$|\mathcal{F}_{\ell-1,+}|-|c_{\ell,+}\bar{z}|\leq|\mathcal{F}_{\ell-1,+}+c_{\ell,+}\bar{z}|<B_{\ell-1,+}\rightarrow|\mathcal{F}_{\ell-1,+}|<B_{\ell-1,+}+|c_{\ell,+}z|.$
(30)

So $\mathcal{F}_{\ell-1,+}$ must be at most a linear function of $z.$ But
since $f_{\ell-1,+}$ has periodic norm, this implies that $c_{\ell,+}=0$ and
$\mathcal{F}_{\ell-1,+}$ is actually a constant. Thus, we conclude that
$f_{\ell,+}=0$ and $f_{\ell-1,+}$ is a constant. Now, using that, we obtain
that $\partial_{\bar{z}}f_{\ell-2,+}=\frac{\gamma_{1}}{2i\hbar
v_{F}}f_{\ell-1,+}=c_{\ell-1,+}.$ Then repeating the line of reasoning above,
we obtain that $c_{\ell-1,+}=0\rightarrow f_{\ell-1,+}=0$ and $f_{\ell-2,+}$
is a constant. This recursive procedure continues until we get to $f_{1,+},$
where we can only conclude that it is a constant but not zero. Next, we study
the $-$ series. This is essentially the same process in reverse with
$z\mapsto\bar{z}.$ From the first condition, $\partial_{z}f_{1,-}=0,$ we get
that $f_{1,-}$ is a constant. Then, $\partial_{z}f_{2,-}=c_{1,-},$ which by
the same argument above, implies that $f_{1,-}=0$ and $f_{2,-}$ is a constant.
The same recursive argument then applies to all other terms showing that
$f_{1,-}=f_{2,-}=f_{3,-}=...=f_{\ell-1,-}=0$ and $f_{\ell,-}$ is a non-zero
constant. Therefore, we conclude that there are only two zero modes with
periodic norm, which can be written explicitly as

$\psi_{+}(\mathbf{r})=\frac{1}{A_{+}}\begin{pmatrix}e^{+e\phi(\mathbf{r})/\hbar}\\\
0\\\ \vdots\\\ 0\\\ 0\end{pmatrix}\text{ and
}\psi_{-}(\mathbf{r})=\frac{1}{A_{-}}\begin{pmatrix}0\\\ 0\\\ \vdots\\\ 0\\\
e^{-e\phi(\mathbf{r})/\hbar}\end{pmatrix}.$ (31)

##  Appendix B Green’s Function of Two-Dimensional Periodic Laplacian

In this section, we provide a brief derivation of the Green’s function of the
two-dimensional periodic Laplacian. This is a textbook result [35]; we only
provide it here to make the manuscript self-contained, and do not claim any
originality in this derivation. We seek the solution for
$G(x,x^{\prime},y,y^{\prime})$

$\frac{\partial^{2}G(x,x^{\prime},y,y^{\prime})}{\partial
x^{2}}+\frac{\partial^{2}G(x,x^{\prime},y,y^{\prime})}{\partial
y^{2}}=\sum_{p=-\infty}^{\infty}\delta^{2}\left(\mathbf{r}-p\hat{x}-\mathbf{r}^{\prime}\right)=\delta(y-y^{\prime})\sum_{p=-\infty}^{\infty}\delta(x-p-x^{\prime}).$
(32)

We perform Fourier transformation using the following convention:

$\begin{split}G(x,x^{\prime},y,y^{\prime})&=\sum_{n=-\infty}^{\infty}\tilde{G}_{n}(y,y^{\prime})e^{2\pi
ni(x-x^{\prime})}\text{ and
}\tilde{G}_{n}(y,y^{\prime})=\int_{-\frac{1}{2}}^{\frac{1}{2}}dxG(x,x^{\prime},y,y^{\prime})e^{-2\pi
ni(x-x^{\prime})},\end{split}$ (33)

which leads to the following ordinary differential equation in reciprocal
space:

$\frac{\partial^{2}\tilde{G}_{n}(y,y^{\prime})}{\partial
y^{2}}-4\pi^{2}n^{2}\tilde{G}_{n}(y,y^{\prime})=\delta(y-y^{\prime}).$ (34)

This can be solved for $y-y^{\prime}<0$ and $y-y^{\prime}>0$ separately, and
then matched at $y-y^{\prime}=0$ for continuity in the function and
discontinuity in the derivative of the function

$\begin{split}\tilde{G}_{0}(y,y^{\prime})&=\frac{1}{2}|y-y^{\prime}|+c,\\\
\tilde{G}_{n\neq
0}(y,y^{\prime})&=-\frac{1}{4\pi|n|}e^{-2\pi|n||y-y^{\prime}|},\end{split}$
(35)

where we have exploited the symmetry
$(y-y^{\prime})\rightarrow-(y-y^{\prime}).$ Now, inverting the Fourier
transform, we obtain

$\begin{split}G(x,x^{\prime},y,y^{\prime})&=\frac{1}{2}|y-y^{\prime}|+c-\sum_{n\neq
0}\frac{1}{4\pi|n|}e^{-2\pi|n||y-y^{\prime}|}e^{2\pi
ni(x-x^{\prime})}.\end{split}$ (36)

Using the following summation identity $-\sum_{n\neq
0}\frac{1}{4\pi|n|}e^{-2\pi|n||y|}e^{2\pi nix}=\frac{\ln
2}{4\pi}-\frac{|y|}{2}+\frac{1}{4\pi}\ln\left[\cosh\left(2\pi|y|\right)-\cos\left(2\pi
x\right)\right],$ we obtain
$G(x,x^{\prime},y,y^{\prime})=\frac{1}{4\pi}\ln\left[\cosh\left(2\pi|y-y^{\prime}|\right)-\cos\left(2\pi|x-x^{\prime}|\right)\right]+c^{\prime}.$
For simplicity, we set $c^{\prime}=0.$ Because the hyperbolic cosine is an
even function, we can drop the absolute value on $y$ to write

$G(x,x^{\prime},y,y^{\prime})=\frac{1}{4\pi}\ln\left[\cosh\left(2\pi\left(y-y^{\prime}\right)\right)-\cos\left(2\pi\left(x-x^{\prime}\right)\right)\right].$
(37)

Obviously, this Green’s function is singular at
$(x+n,y)=(x^{\prime},y^{\prime}).$ Direct calculation confirms that
$\Delta_{\mathbf{r}}G(x,x^{\prime},y,y^{\prime})=0$ everywhere else. For a
function $\phi(x,y)$ satisfying the Poisson’s equation

$\frac{\partial^{2}\phi(x,y)}{\partial
x^{2}}+\frac{\partial^{2}\phi(x,y)}{\partial y^{2}}=B(x,y),$ (38)

we can write its formal solution has a convolution with the Green’s function

$\phi(x,y)=\int_{-\infty}^{\infty}dy^{\prime}\int_{-\frac{1}{2}}^{\frac{1}{2}}dx^{\prime}B(x^{\prime},y^{\prime})G(x,x^{\prime},y,y^{\prime}).$
(39)

We do not worry much about the boundary condition of $\phi(x,y).$ We only
require that it be periodic in $x,$ which is clear since
$G(x+n,x^{\prime},y,y^{\prime})=G(x,x^{\prime},y,y^{\prime}):$

$\phi(x+n,y)=\int_{-\infty}^{\infty}dy^{\prime}\int_{-\frac{1}{2}}^{\frac{1}{2}}dx^{\prime}B(x^{\prime},y^{\prime})G(x+n,x^{\prime},y,y^{\prime})=\int_{-\infty}^{\infty}dy^{\prime}\int_{-\frac{1}{2}}^{\frac{1}{2}}dx^{\prime}B(x^{\prime},y^{\prime})G(x,x^{\prime},y,y^{\prime})=\phi(x,y).$
(40)

## References

  * Ashcroft  _et al._ [2016] N. Ashcroft, N. Mermin, and D. Wei,  _Solid State Physics_ (Cengage Learning, 2016).
  * Jackson [1998] J. Jackson,  _Classical Electrodynamics_ (Wiley, 1998).
  * Paul [1990] W. Paul, Electromagnetic traps for charged and neutral particles, Rev. Mod. Phys. 62, 531 (1990).
  * Wieman  _et al._ [1999] C. E. Wieman, D. E. Pritchard, and D. J. Wineland, Atom cooling, trapping, and quantum manipulation, Rev. Mod. Phys. 71, S253 (1999).
  * Fortágh and Zimmermann [2007] J. Fortágh and C. Zimmermann, Magnetic microtraps for ultracold atoms, Rev. Mod. Phys. 79, 235 (2007).
  * Chen [2013] F. Chen,  _Introduction to Plasma Physics and Controlled Fusion: Volume 1: Plasma Physics_ (Springer US, 2013).
  * Bellan [2008] P. Bellan,  _Fundamentals of Plasma Physics_ (Cambridge University Press, 2008).
  * Aharonov and Casher [1979] Y. Aharonov and A. Casher, Ground state of a spin-½ charged particle in a two-dimensional magnetic field, Phys. Rev. A 19, 2461 (1979).
  * Dubrovin and Novikov [1980] B. A. Dubrovin and S. P. Novikov, Ground states of a two-dimensional electron in periodic magnetic field, Zh. Eksp. Teor. Fiz. 79, 1006 (1980).
  * Rozenblum and Shirokov [2006] G. Rozenblum and N. Shirokov, Infiniteness of zero modes for the pauli operator with singular magnetic field, J. Funct. Anal. 233, 135 (2006).
  * Mi  _et al._ [2015] H. Mi, S. Mikael, C.-C. Liu, J.-H. Seo, G. Gui, A. L. Ma, P. F. Nealey, and Z. Ma, Creating periodic local strain in monolayer graphene with nanopillars patterned by self-assembled block copolymer, Appl. Phys. Lett. 107, https://doi.org/10.1063/1.4932657 (2015).
  * Banerjee  _et al._ [2020] R. Banerjee, V.-H. Nguyen, T. Granzier-Nakajima, L. Pabbi, A. Lherbier, A. R. Binion, J.-C. Charlier, M. Terrones, and E. W. Hudson, Strain modulated superlattices in graphene, Nano Lett. 20, 3113 (2020).
  * Mao  _et al._ [2020] J. Mao, S. P. Milovanović, M. Anđelković, X. Lai, Y. Cao, K. Watanabe, T. Taniguchi, L. Covaci, F. M. Peeters, A. K. Geim, _et al._ , Evidence of flat bands and correlated states in buckled graphene superlattices, Nature 584, 215 (2020).
  * Vozmediano  _et al._ [2010] M. A. Vozmediano, M. Katsnelson, and F. Guinea, Gauge fields in graphene, Phys. Rep. 496, 109 (2010).
  * Milovanović  _et al._ [2020] S. P. Milovanović, M. Anđelković, L. Covaci, and F. M. Peeters, Band flattening in buckled monolayer graphene, Phys. Rev. B 102, 245427 (2020).
  * Manesco  _et al._ [2020] A. L. Manesco, J. L. Lado, E. V. Ribeiro, G. Weber, and D. Rodrigues Jr, Correlations in the elastic landau level of spontaneously buckled graphene, 2D Mater. 8, 015011 (2020).
  * Manesco and Lado [2021] A. L. Manesco and J. L. Lado, Correlation-induced valley topology in buckled graphene superlattices, 2D Mater. 8, 035057 (2021).
  * Phong and Mele [2022] V. T. Phong and E. J. Mele, Boundary modes from periodic magnetic and pseudomagnetic fields in graphene, Phys. Rev. Lett. 128, 176406 (2022).
  * De Beule  _et al._ [2023a] C. De Beule, V. T. Phong, and E. J. Mele, Network model for periodically strained graphene, Phys. Rev. B 107, 045405 (2023a).
  * Gao  _et al._ [2023] Q. Gao, J. Dong, P. Ledwith, D. Parker, and E. Khalaf, Untwisting moiré physics: Almost ideal bands and fractional chern insulators in periodically strained monolayer graphene, Phys. Rev. Lett. 131, 096401 (2023).
  * Mahmud  _et al._ [2023] M. T. Mahmud, D. Zhai, and N. Sandler, Topological flat bands in strained graphene: Substrate engineering and optical control, Nano Lett. 23, 7725 (2023).
  * Andrade  _et al._ [2023] E. Andrade, F. López-Urías, and G. G. Naumis, Topological origin of flat bands as pseudo-landau levels in uniaxial strained graphene nanoribbons and induced magnetic ordering due to electron-electron interactions, Phys. Rev. B 107, 235143 (2023).
  * Wan  _et al._ [2023] X. Wan, S. Sarkar, K. Sun, and S.-Z. Lin, Nearly flat chern band in periodically strained monolayer and bilayer graphene, Phys. Rev. B 108, 125129 (2023).
  * De Beule  _et al._ [2023b] C. De Beule, V. T. Phong, and E. J. Mele, Rose patterns in the nonperturbative current response of two-dimensional superlattices, arXiv preprint arXiv:2305.03013 (2023b).
  * Dong  _et al._ [2022] J. Dong, J. Wang, and L. Fu, Dirac electron under periodic magnetic field: Platform for fractional chern insulator and generalized wigner crystal, arXiv preprint arXiv:2208.10516 (2022).
  * Mao and Chowdhury [2023] D. Mao and D. Chowdhury, Upper bounds on superconducting and excitonic phase-stiffness for interacting isolated narrow bands, arXiv preprint arXiv:2304.07318 (2023).
  * Erdős and Vougalter [2002] L. Erdős and V. Vougalter, Pauli operator and aharonov–casher theorem for measure valued magnetic fields, Comm. Math. Phys. 225, 399 (2002).
  * Kailasvuori [2009] J. Kailasvuori, Pedestrian index theorem à la aharonov-casher for bulk threshold modes in corrugated multilayer graphene, EPL 87, 47008 (2009).
  * Note [1] One can equivalently consider a Pauli Hamiltonian, which is $\mathcal{H}_{1}^{2}$ in our notation. The argument proceeds in the exact same manner for the Pauli Hamiltonian.
  * Note [2] More precisely, $\partial_{\bar{z}}f_{+}=0$ and $\partial_{z}f_{-}=0.$ So $f_{+}$ is holomorphic with respect to $z=x+iy$ and $f_{-}$ is holomorphic with respect to $\bar{z}=x-iy$, which is usually called antiholomorphic in complex analysis. Therefore, both functions have power-series representations  $\begin{split}f_{+}(z)=\sum@\slimits@_{n=0}^{\infty}a_{n}z^{n}\hskip 10.00002pt\text{and}\hskip 10.00002ptf_{-}(\bar{z})=\sum@\slimits@_{n=0}^{\infty}b_{n}\bar{z}^{n},\end{split}$ (41) where $a_{n}\in\mathbb{C}$ and $b_{n}\in\mathbb{C}$ are constants.
  * [31] Please see Refs [50, 51, 52, 53, 54, 55] and the references contained therein.
  * Note [3] Though this assumption may seem strict, it is often employed in numerical studies due to the impracticality of including an infinite number of wavevectors.
  * Min and MacDonald [2008] H. Min and A. H. MacDonald, Electronic structure of multilayer graphene, Prog. Theor. Phys. Supp. 176, 227 (2008).
  * Khalaf  _et al._ [2019] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Magic angle hierarchy in twisted graphene multilayers, Phys. Rev. B 100, 085109 (2019).
  * Ammari  _et al._ [2018] H. Ammari, B. Fitzpatrick, H. Kang, M. Ruiz, S. Yu, and H. Zhang,  _Mathematical and computational methods in photonics and phononics_ , Vol. 235 (American Mathematical Society, 2018).
  * Ye  _et al._ [1995] P. D. Ye, D. Weiss, R. R. Gerhardts, M. Seeger, K. von Klitzing, K. Eberl, and H. Nickel, Electrons in a periodic magnetic field induced by a regular array of micromagnets, Phys. Rev. Lett. 74, 3013 (1995).
  * Hunziker [1980] W. Hunziker, Schrödinger operators with electric or magnetic fields, in  _Mathematical Problems in Theoretical Physics_ , edited by K. Osterwalder (Springer Berlin Heidelberg, Berlin, Heidelberg, 1980) pp. 25–44.
  * Xue and Xiao [1992] D. P. Xue and G. Xiao, Magnetotransport properties of two-dimensional electron gases under a periodic magnetic field, Phys. Rev. B 45, 5986 (1992).
  * Peeters and Vasilopoulos [1993] F. M. Peeters and P. Vasilopoulos, Quantum transport of a two-dimensional electron gas in a spatially modulated magnetic field, Phys. Rev. B 47, 1466 (1993).
  * Chang and Niu [1994] M. C. Chang and Q. Niu, Electron band structure in a two-dimensional periodic magnetic field, Phys. Rev. B 50, 10843 (1994).
  * Carmona  _et al._ [1995] H. A. Carmona, A. K. Geim, A. Nogaret, P. C. Main, T. J. Foster, M. Henini, S. P. Beaumont, and M. G. Blamire, Two dimensional electrons in a lateral magnetic superlattice, Phys. Rev. Lett. 74, 3009 (1995).
  * Ibrahim and Peeters [1995] I. S. Ibrahim and F. M. Peeters, Two-dimensional electrons in lateral magnetic superlattices, Phys. Rev. B 52, 17321 (1995).
  * Krakovsky [1996] A. Krakovsky, Electronic band structure in a periodic magnetic field, Phys. Rev. B 53, 8469 (1996).
  * Li  _et al._ [1996] T.-Z. Li, S.-W. Gu, X.-H. Wang, and J.-P. Peng, Electrical properties of a two-dimensional electron gas under a general one-dimensional periodic magnetic field, J. Phys. Condens. Matter 8, 313 (1996).
  * Miró  _et al._ [2014] P. Miró, M. Audiffred, and T. Heine, An atlas of two-dimensional materials, Chem. Soc. Rev. 43, 6537 (2014).
  * Wehling  _et al._ [2014] T. O. Wehling, A. M. Black-Schaffer, and A. V. Balatsky, Dirac materials, Adv. Phys. 63, 1 (2014).
  * Wang  _et al._ [2015] J. Wang, S. Deng, Z. Liu, and Z. Liu, The rare two-dimensional materials with dirac cones, Natl. Sci. Rev. 2, 22 (2015).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Fu [2011] L. Fu, Topological crystalline insulators, Phys. Rev. Lett. 106, 106802 (2011).
  * Shigekawa [1991] I. Shigekawa, Spectral properties of schrödinger operators with magnetic fields for a spin 12 particle, J. Funct. Anal. 101, 255 (1991).
  * Arai [1993] A. Arai, Properties of the dirac–weyl operator with a strongly singular gauge potential, J. Math. Phys. 34, 915 (1993).
  * Hirokawa and Ogurisu [2001] M. Hirokawa and O. Ogurisu, Ground state of a spin-1/2 charged particle in a two-dimensional magnetic field, J. Math. Phys 42, 3334 (2001).
  * Persson [2008] M. Persson, Zero modes for the magnetic pauli operator in even-dimensional euclidean space, Lett. Math. Phys. 85, 111 (2008).
  * Bony  _et al._ [2019] J.-F. Bony, N. Espinoza, and G. Raikov, Spectral properties of 2d pauli operators with almost-periodic electromagnetic fields, Publ. Res. Inst. Math. Sci. 55, 453 (2019).
  * Fialová [2023] M. Fialová, Aharonov-casher theorems for manifolds with boundary and aps boundary condition, arXiv preprint arXiv:2304.13373 https://arxiv.org/abs/2304.13373 (2023).


# Horospherical Decision Boundaries for Large Margin Classification in
Hyperbolic Space

Xiran Fan
Department of Statistics
University of Florida
<EMAIL_ADDRESS>
&Chun-Hao Yang
Institute of Statistics and Data Science
National Taiwan University
<EMAIL_ADDRESS>
Baba C. Vemuri
Department of CISE
University of Florida
<EMAIL_ADDRESS>

###### Abstract

Hyperbolic spaces have been quite popular in the recent past for representing
hierarchically organized data. Further, several classification algorithms for
data in these spaces have been proposed in the literature. These algorithms
mainly use either hyperplanes or geodesics for decision boundaries in a large
margin classifiers setting leading to a non-convex optimization problem. In
this paper, we propose a novel large margin classifier based on horospherical
decision boundaries that leads to a geodesically convex optimization problem
that can be optimized using any Riemannian gradient descent technique
guaranteeing a globally optimal solution. We present several experiments
depicting the competitive performance of our classifier in comparison to SOTA.

##  1 Introduction

Hyperbolic space, a non-Euclidean space with constant negative curvature, has
been shown [25, 23, 29, 24] to be effective for representing hierarchically
organized data. For example, authors in [25] showed that a tree can be
embedded in a hyperbolic space with arbitrarily small distortion. The main
reason for this is that a hyperbolic space can be regarded as a continuous
version of trees – the volume of the space grows _exponentially_ as one moves
away from the center in hyperbolic space. This matches the growth pattern of
the number of nodes in a tree which grows _exponentially_ as the depth of the
tree increases. Hyperbolic space embedding has been shown to be a promising
approach for representing data with a (latent) hierarchical structure [23, 24,
29, 15].

Recently, representation of data in hyperbolic space for the fundamental tasks
of unsupervised and supervised learning has been popularized in various
contexts, e.g., dimensionality reduction [10, 13], clustering [22], large-
margin classifier [12, 33, 11], regression [20], etc. Existing ’linear’
classifiers in hyperbolic spaces are predominantly based on _geodesics_ i.e.,
using geodesics as decision boundaries. In [12], the decision boundary is
chosen to be the intersection of the hyperboloid model and a hyperplane in the
ambient space, which in this case is the Minkowski space. Then, the support
vector machine (SVM) in hyperbolic space is formulated as a nonconvex
optimization problem. In [33], authors followed the same parameterization of
the hyperbolic geodesic decision plane and provided a series of algorithms to
provably learn large margin classifiers in hyperbolic space. However, as
pointed out by [11], the algorithm in [33] fails to converge in practice.
Authors in [11] used the Poincaré ball model and parameterized the geodesic
decision plane as a hyperplane mapped using the exponential map from the
_tangent space_ at some reference point. They first constructed convex hulls
for each data cluster in the hyperbolic space and the reference point is then
chosen to be the midpoint between different convex hulls. Then they apply a
_Euclidean_ perceptron/SVM algorithm to data lifted in to the tangent space at
the aforementioned reference point. Although the optimization problem in the
tangent space is convex, the procedure of the tangent space approximation
introduces inaccuracies and distortions. Moreover, convex hull learning is
highly unstable and their implementation is only applicable to the
2-dimensional hyperbolic space.

Finally, it is worth mentioning that linear classification within hyperbolic
space, which can be considered as the last layer, referred to as the
hyperbolic logistic regression (LR), is a fundamental component of hyperbolic
neural networks (HNNs) [16, 27]. The calculation of logits in this layer is
based on the distances between samples and the geodesic decision boundary.
Notably, this hyperbolic LR employs a geodesic decision boundary but is not a
large-margin classifier.

Figure 1: A 2-d Poincaré disk model $\mathbb{B}^{2}$ and its boundary
$\partial\mathbb{B}^{2}=\mathbb{S}^{1}$. Given an ideal point
$\boldsymbol{\omega}\in\partial\mathbb{B}^{2}$, the black lines and curves are
hyperbolic geodesics starting (ending) at $\boldsymbol{\omega}$ and the red
circles are horocycles centered at $\boldsymbol{\omega}$.

###  1.1 Horospherical Decision Boundaries for Classification in Hyperbolic
Sapce

_Horospheres_ , which are the level sets of the _Busemann function_ in
hyperbolic spaces, are the analogs of Euclidean hyperplanes [4]. Horospheres
(horocycles) are contained in the Poincaré ball (disk) and are tangential to
the ball (disk) at an ideal point as shown in Figure 1. A collection of
horospheres centered at the same ideal point are parallel to each other and
the lengths of geodesic segments between two horospheres are all equal, just
as the lengths of line segments between parallel hyperplanes in Euclidean
space are all equal. This property of horospheres was explored by [9] to
develop a dimensionality reduction method for data in hyperbolic space. By
using horospherical projection, they are able to preserve the distance
information in the original data. However, there is no literature on
constructing a ‘linear’ classifier in hyperbolic space using horospheres as
the decision boundaries although the horospheres are the hyperbolic equivalent
of Euclidean hyperplanes. Therefore, it is natural to consider the use of
horospheres as decision boundaries for classification in hyperbolic spaces. In
this work, we propose a novel _hyperbolic large-margin classifier using
horospheres as decision boundaries in the Poincaré model._ We term this
classifier as a _HoroSVM_. The horospheres are well-defined in the Poincaré
ball model. A toy example as shown in Figure 2 demonstrates the advantage of
horospheres decision boundaries over geodesic decision boundaries. As the
tree-structured data grows in depth, leaf nodes are embedded closer to each
other within a subtree and among different subtrees. One of the classification
problems in hyperbolic space is to determine whether a node belongs to a
chosen subtree given the embedding. For comparison purposes, the decision
boundaries of HoroSVM (Figure 2(a)) and hyperboloid SVM [12] (Figure 2(b)) are
shown in the figure. As evident, the horospheres decision boundary perfectly
separates (the root node is excluded in training) the data while the geodesic
decision boundary makes several mistakes on both positive and negative
samples. We present a novel formulation of the classification problem in the
hyperbolic space as a geodesically convex optimization problem on a Riemannian
manifold. This optimization problem can be easily solved using any Riemannian
gradient descent technique guaranteeing global optimality. Gradient-based
optimizations for geodesically convex problems guaranteeing global optimal
solutions are the topic of investigation in optimization literature and we
refer the reader to [36] for detailed convergence analysis of several such
optimization methods. Further, we empirically validate our method on several
real and synthetic data sets.

(a) Horocycle decision boundary

(b) Geodesic decision boundary

Figure 2: A balanced tree with depth 6 and spread 4 embedded in a 2-d Poincaré
ball model using [15] is depicted in the figure. The orange plus node is the
root of a chosen subtree, the blue dots are positive samples (nodes of the
subtree) and the red dots are negative samples. (a) depicts a HoroSVM
performance on the classification of positive and negative samples/nodes along
with a zoomed-in version on its right. (b) depicts the geodesic boundary from
the competing method, Hyperboloid SVM [12], along with the zoomed-in version
to its right.

It should be noted that a horosphere decision boundary has been used in some
recent works [32, 28] in constructing HNNs. For example, authors in [32]
proposed hyperbolic neuron models using the Busemann function as a
generalization of the Euclidean inner product to extract horosphere features
from data. Authors in [28] proposed a shallow fully-connected continuous
network spanned by (hyperbolic) neurons, on noncompact symmetric space
(including hyperbolic space) using the Helgason-Fourier transform. Authors in
[34] produced Euclidean features from hyperbolic embeddings via the
eigenfunctions of the Laplace operator in the hyperbolic space where the
eigenfunctions involve horosphere features. Note that none of the above works
developed a large margin classifier using the horosphere as a decision
boundary. To the best of our knowledge, our work is the first in the
literature to present a convex optimization formulation of a large-margin
classifier using a horosphere decision boundary in a hyperbolic space.

The rest of this paper is organized as follows. In Section 2, we present some
background on hyperbolic geometry pertinent to the work presented here. In
Section 3, we present our horospherical boundary-based classification methods.
Experimental results are presented in Section 4 to demonstrate the advantage
of our HoroSVM over competing hyperbolic classifiers. Finally, we conclude in
Section 5.

##  2 Background

In this section, we review some basic concepts of hyperbolic geometry
including the generalization of the Euclidean hyperplane to the hyperbolic
space namely, the horosphere.

###  2.1 Hyperbolic Space and the Poincaré Ball Model

There are five isometric models of the hyperbolic space: the Poincaré ball
model, the Lorentz model, the Klein model, the upper-half space model, and the
Hemisphere model [7]. We choose the Poincaré Ball model in this paper as it is
easy to visualize and the Busemann function has a nice closed-form expression
in this model. Note that the decision boundary of choice in our work is the
level-set of the Busemann function namely, the horosphere.

An $n$-dimensional Poincaré Ball model, denoted by
$(\mathbb{B}^{n},g_{\mathbb{B}})$, consists of all points in an open ball of
radius 1, i.e.,
$\mathbb{B}^{n}=\\{\boldsymbol{x}\in\mathbb{R}^{n}:\lVert\boldsymbol{x}\rVert<1\\}$,
and equipped with the Riemannian metric
$g_{\mathbb{B}}(\boldsymbol{x})=4(1-\lVert\boldsymbol{x}\rVert^{2})^{-2}g_{\mathbb{R}}$,
where $\lVert\cdot\rVert$ is the Euclidean $L_{2}$ norm, and $g_{\mathbb{R}}$
is the Euclidean metric. The geodesic distance between points
$\boldsymbol{x},\boldsymbol{y}\in\mathbb{B}^{n}$ is
$d_{\mathbb{B}}(\boldsymbol{x},\boldsymbol{y})=\cosh^{-1}\left(1+2\frac{\lVert\boldsymbol{x}-\boldsymbol{y}\rVert^{2}}{(1-\lVert\boldsymbol{x}\rVert^{2})(1-\lVert\boldsymbol{y}\rVert^{2})}\right)$.

###  2.2 Horospheres

#### Geodesics, geodesic rays, and ideal points

The shortest path that connects two points in the Poincaré Ball model is
called a _geodesic segment_. A _geodesic ray_ is a geodesic segment that can
be infinitely extended in one direction. We call the endpoint at infinity of a
geodesic ray an _ideal point_. For $\mathbb{B}^{n}$, ideal points form the
boundary of the ball:
$\partial\mathbb{B}^{n}=\mathbb{S}^{n-1}=\\{x\in\mathbb{R}^{n}:\lVert\boldsymbol{x}\rVert=1\\}$,
where $\mathbb{S}^{n-1}$ is the $(n-1)$-dimensional hypersphere. The
hypersphere $(\mathbb{S}^{n-1},g_{\mathbb{S}})$ is a Riemannian manifold
equipped with the Riemannian metric
$g_{\mathbb{S}}=4(1+\lVert\boldsymbol{x}\rVert^{2})^{-2}g_{\mathbb{R}}$.

#### Busemann function [6]

Let $\boldsymbol{\omega}\in\partial\mathbb{B}^{n}$ be an ideal point and
$\gamma_{\boldsymbol{\omega}}:[0,\infty)\rightarrow\mathbb{B}^{n}$ a geodesic
ray pointing $\boldsymbol{\omega}$. The _Busemann function_ is defined as

$b_{\boldsymbol{\omega}}(\boldsymbol{x})=\lim_{t\rightarrow\infty}(d(\gamma_{\boldsymbol{\omega}}(t),\boldsymbol{x})-t),\quad\boldsymbol{x}\in\mathbb{B}^{n}.$
(1)

In the Poincaré Ball model, Eq. (1) has a closed form :
$b_{\boldsymbol{\omega}}(\boldsymbol{x})=-\log\frac{1-\lVert\boldsymbol{x}\rVert^{2}}{\lVert\boldsymbol{\omega}-\boldsymbol{x}\rVert^{2}}.$

#### Horospheres [6]

In $\mathbb{\mathbb{B}}^{n}$, a horosphere is a $(n-1)$-dimensional sphere
that is internally tangent to $\partial\mathbb{B}^{n}$ at an ideal point. For
a given $\boldsymbol{\omega}\in\partial\mathbb{B}^{n}$, the level sets of
Busemann function $b_{\boldsymbol{\omega}}(\boldsymbol{x})$ in the Poincaré
ball model is a series of horospheres tangent at $\boldsymbol{\omega}$.
Horospheres are hyperbolic hyperplanes in the sense that the corresponding
construction in Euclidean space gives a hyperplane.

#### Poincaré inner product

The Busemann function can be regarded as the inner product between a direction
vector and a point in the Poincaré ball model. We call it Poincaré inner
product and write it as
$\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}=\log\frac{1-\lVert\boldsymbol{x}\rVert^{2}}{\lVert\boldsymbol{\omega}-\boldsymbol{x}\rVert^{2}}$
111Note that the sign is opposite as in the Busemann function.  .

Given $\boldsymbol{\omega}\in\partial\mathbb{B}^{n}$, the Poincaré inner
product $\langle\boldsymbol{\omega},\cdot\rangle_{\mathbb{B}}$ is constant
over horosphere tangent at $\boldsymbol{\omega}$. Note that the Euclidean
inner product $\langle\boldsymbol{w},\cdot\rangle$ is constant over a
hyperplane that is perpendicular to a given direction $\boldsymbol{w}$. Let
$\Pi$ denote the set of horospheres of $\mathbb{B}^{n}$. A horosphere
$\pi\in\Pi$ can be parameterized by $0<\mu\in\mathbb{R}^{+}$,
$\boldsymbol{\omega}\in\mathbb{S}^{n-1}$, and $b\in\mathbb{R}$ as

$\pi_{\mu,\boldsymbol{\omega},b}:=\\{\boldsymbol{z}\in\mathbb{B}^{n}|\mu\langle\boldsymbol{\omega},\boldsymbol{z}\rangle_{\mathbb{B}}-b=0\\}.$
(2)

We will use $\pi$ or $\pi_{\cdot}$ to represent a horosphere in different
parameterizations hereafter.

##  3 Horosphere Boundary-based Classification

In this section, we present the key theoretical contributions of our work
namely, a horosphere-based SVM classifier that involves formulating and
solving a geodesically convex optimization problem. First, we present some
preliminary facts and results about the horospheres. Then, we present the
optimization problems for the horospherical perceptron and SVM along with
analysis.

###  3.1 Point to Horocycle Distance

While the distance from the origin $\boldsymbol{o}$ to a horosphere has been
known for decades [18] (Introduction 4.1, p.31), we present a natural
generalization of previous results by providing a closed-form expression for
measuring the hyperbolic distance from any arbitrary point
$\boldsymbol{x}\in\mathbb{B}^{n}$ to a given horosphere
$\pi_{\mu,\boldsymbol{\omega},b}$. The following remark provides this result.

######  Proposition 3.1.

Let $\pi_{\mu,\boldsymbol{\omega},b}$ be a horosphere. The hyperbolic distance
of a point $\boldsymbol{x}\in\mathbb{B}^{n}$ to a horosphere
$\pi_{\mu,\boldsymbol{\omega},b}$ is given by

$d_{\mathbb{B}}(\boldsymbol{x},\pi_{\mu,\boldsymbol{\omega},b})=\frac{\lvert\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b\rvert}{\mu}.$
(3)

Notice that it shares a similarity to the Euclidean distance of a point to a
hyperplane. Before presenting the proof for Proposition 3.1, we recall the
following Fact 3.2 and Lemma 3.3 from [32].

######  Fact 3.2.

Given an ideal point $\boldsymbol{\omega}$ and a point
$\boldsymbol{x}\in\mathbb{B}^{n}$, there is a unique horosphere passing
through $\boldsymbol{x}$ and tangent at $\boldsymbol{\omega}$.

######  Lemma 3.3.

[32] Let $\Pi_{\boldsymbol{\omega}}$ be the set of horocycles of
$\mathbb{B}^{n}$ tangent at $\boldsymbol{\omega}$. Given
$\lambda\in\mathbb{R}$, let $\pi_{\lambda,\boldsymbol{\omega}}$ be the unique
horosphere that passes through $\tanh(\lambda/2)\cdot\boldsymbol{\omega}$ and
tangent at $\boldsymbol{\omega}$. Note that
$\Pi_{\boldsymbol{\omega}}=\cup_{\lambda\in\mathbb{R}}\\{\pi_{\lambda,\boldsymbol{\omega}}\\}$.
We have the following two results: (i) the hyperbolic lengths of geodesic
(that pass through $\boldsymbol{\omega}$) segments between
$\pi_{\lambda_{1},\boldsymbol{\omega}}$ and
$\pi_{\lambda_{2},\boldsymbol{\omega}}$ are equal to
$\lvert\lambda_{1}-\lambda_{2}\rvert$; (ii)
$\langle\boldsymbol{\omega},x\rangle_{\mathbb{B}}=\lambda$ for any
$\boldsymbol{x}\in\pi_{\lambda,\boldsymbol{\omega}}$.

Figure 3 shows a 2D Poincaré disk model $\mathbb{B}^{2}$ and its boundary
$\mathbb{S}^{1}$. The point $o$ is the origin of the disk,
$\boldsymbol{x}\in\mathbb{B}^{2}$ is a point, and
$\boldsymbol{\omega}\in\mathbb{S}^{1}$ is a point at infinity (an ideal
point). Two geodesics ending at the same $\boldsymbol{\omega}$ from
$\boldsymbol{x}$ and $o$ respectively are shown in the figure (black solid
line/curve). The circle $\pi_{\mu,\boldsymbol{\omega},b}$ is a given horocycle
tangent (red solid circle) at $\boldsymbol{\omega}$. The hyperbolic distance
$d_{\mathbb{B}}(\boldsymbol{x},\pi_{\mu,\boldsymbol{\omega},b})$ between
$\boldsymbol{x}$ and $\pi_{\mu,\boldsymbol{\omega},b}$ is identified as the
distance between $\boldsymbol{x}$ and $\boldsymbol{y}_{\boldsymbol{x}}$ where
$\boldsymbol{y}_{\boldsymbol{x}}$ is the projection of $\boldsymbol{x}$ to
$\pi_{\mu,\boldsymbol{\omega},b}$ along the geodesic ending at
$\boldsymbol{\omega}$. Let $\pi^{\boldsymbol{x}}$ (red dashed circle) be the
unique horocycle that passes through $\boldsymbol{x}$ and is tangent at
$\boldsymbol{\omega}$. Note that the lengths of all geodesic segments between
two horocycles are the same. That is,
$d_{\mathbb{B}}(\boldsymbol{x},\pi_{\mu,\boldsymbol{\omega},b})=d_{\mathbb{B}}(\boldsymbol{x},\boldsymbol{y}_{\boldsymbol{x}})=d_{\mathbb{B}}(\boldsymbol{x}_{0},\boldsymbol{y}_{0})$,
where $\boldsymbol{x}_{0},\boldsymbol{y}_{0}$ are _horocyclic projections_ [9]
of $\boldsymbol{x},\boldsymbol{y}_{\boldsymbol{x}}$ along
$\pi^{\boldsymbol{x}}$ and $\pi_{\mu,\boldsymbol{\omega},b}$ respectively.

Figure 3:  The relationship between horocycles, geodesic, and the horocyclic
projections in $\mathbb{B}^{2}$.

###### Proof of Proposition 3.1.

Given a point $\boldsymbol{x}\in\mathbb{B}^{n}$ and an ideal point
$\boldsymbol{\omega}\in\mathbb{S}^{n-1}$, let
$\pi\in\Pi_{\boldsymbol{\omega}}$ be a horosphere tangent at
$\boldsymbol{\omega}$ and $\pi^{\boldsymbol{x}}\in\Pi_{\boldsymbol{\omega}}$
be the unique horosphere that passes through $\boldsymbol{x}$. Write
$\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}=\lambda_{\boldsymbol{x}}$
and then
$\boldsymbol{x}_{0}=\tanh(\lambda_{\boldsymbol{x}}/2)\cdot\boldsymbol{\omega}$
is the horospherical projection of $\boldsymbol{x}$ along
$\pi_{\boldsymbol{x}}$. For consistency in notation, we write
$\pi^{\boldsymbol{x}}$ as
$\pi_{\lambda_{\boldsymbol{x}},\boldsymbol{\omega}}$. Since the horosphere
$\pi_{\mu,\boldsymbol{\omega},b}$ can be reparameterized as
$\pi_{\lambda,\boldsymbol{\omega}}$, where $\lambda=\frac{b}{\mu}$, the
hyperbolic distance from $\boldsymbol{x}$ to $\pi_{\mu,\boldsymbol{\omega},b}$
is the hyperbolic distance between
$\pi_{\lambda_{\boldsymbol{x}},\boldsymbol{\omega}}$ and
$\pi_{\lambda,\boldsymbol{\omega}}$ which is
$\lvert\lambda_{\boldsymbol{x}}-\lambda\rvert=\frac{\lvert\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b\rvert}{\mu}$
(by Lemma 3.3). This completes the proof. ∎

###  3.2 Horospherical Decision Boundaries

We consider classification problems in hyperbolic space of the following form:
$\mathcal{X}\subset\mathbb{B}^{n}$ denotes the feature space and
$\mathcal{Y}=\\{\pm 1\\}$ denotes the binary label space. In the following, we
denote the training set by $S\subset\mathcal{X}\times\mathcal{Y}$. The
decision rule using a horosphere as its decision boundary can be written as
the following function $f:\mathcal{X}\mapsto\mathcal{Y}$ where

$f(\boldsymbol{x};\mu,\boldsymbol{\omega},b)=\text{sign}\left(\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b\right).$
(4)

The positive samples are expected to lie inside a horosphere while the
negative samples are expected to lie outside a horosphere. This is analogous
to the linear decision boundary in Euclidean space and we will build a
horospherical perceptron and a horospherical SVM based on this decision
boundary.

It should be noted that in $\mathbb{R}^{n}$ the hyperplane
$\xi_{a,\boldsymbol{w},b}=\\{\boldsymbol{z}\in\mathbb{R}^{n}|a\langle\boldsymbol{w},\boldsymbol{z}\rangle-b=0\\}$
where $a\in\mathbb{R}^{+}$,$b\in\mathbb{R}$,
$\boldsymbol{w}\in\mathbb{S}^{n-1}$ is the hyperplane
$\xi_{a,-\boldsymbol{w},-b}$. However, the horospheres
$\pi_{\mu,\boldsymbol{\omega},b}$ and $\pi_{\mu,-\boldsymbol{\omega},-b}$
respectively represent two distinct horospheres, centered at
$\boldsymbol{\omega}$ and $-\boldsymbol{\omega}$ respectively. Let
$\Pi^{+}=\\{\pi_{\mu,\boldsymbol{\omega},b}\in\Pi|b>0\\}$ and
$\Pi^{-}=\\{\pi_{\mu,\boldsymbol{\omega},b}\in\Pi|b<0\\}$. Thus,
$\Pi^{+},\Pi^{-}\subset\Pi$ and the radius of $\pi\in\Pi^{+}$ is less than
$1/2$ and the radius of $\pi\in\Pi^{-}$ is greater than $1/2$. In most cases
of classification in hyperbolic space, the positive samples are clustered near
the boundaries. Hence, we restrict ourselves to finding a horosphere
$\pi\in\Pi^{+}$ that separates data, instead of searching over $\Pi$.
Intuitively, we are looking for a ‘small’ horosphere that captures the
positive samples. We are now ready to present the Horospherical Perceptron
followed by the Horospherical SVM.

###  3.3 Horospherical Perceptron

The loss function for the proposed horospherical perceptron is given by

$l(\mu,\boldsymbol{\omega},b;\boldsymbol{x},y)=\max(0,-y\cdot\left(\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b)\right),\quad(\mu,\boldsymbol{\omega},b)\in\mathbb{R}^{+}\times\mathbb{S}^{n-1}\times\mathbb{R}^{+},$
(5)

which is zero when the instance is classified correctly and is proportional to
the signed distance of the instance from the horosphere when it is
misclassified. The empirical loss for a given data set $S$ is

$L(\mu,\boldsymbol{\omega},b)=\frac{1}{|S|}\sum_{\\{\boldsymbol{x},y\\}\in
S}l(\mu,\boldsymbol{\omega},b;\boldsymbol{x},y)$ (6)

Hence, the optimal horosphere is learned by solving the above optimization
problem on the manifold
$\mathbb{R}^{+}\times\mathbb{S}^{n-1}\times\mathbb{R}^{+}$, i.e.,

$\mu^{*},\boldsymbol{\omega}^{*},b^{*}=\arg\min_{(\mu,\boldsymbol{\omega},b)}L(\mu,\boldsymbol{\omega},b)$
(7)

To further analyze this optimization problem, we first recall some facts about
geodesic convexity.

######  Definition 3.4.

(Geodesically convex sets [30]). Let $(\mathcal{M},g)$ be a Riemannian
manifold. A set $\boldsymbol{A}\subseteq\mathcal{M}$ is said to be a
geodesically convex set if, for any two points
$\boldsymbol{p},\boldsymbol{q}\in\boldsymbol{A}$, the geodesic
$\gamma_{\boldsymbol{p}\boldsymbol{q}}$ that connects them is contained in
$\boldsymbol{A}$.

######  Definition 3.5.

(Geodesically convex/concave functions [30]) Let
$\boldsymbol{A}\subseteq\mathcal{M}$ be a geodesically convex set. A function
$f:\boldsymbol{A}\rightarrow\mathbb{R}$ is said to be a geodesically convex
function if, for any $\boldsymbol{p},\boldsymbol{q}\in\boldsymbol{A}$ the
composition
$f\circ\gamma_{\boldsymbol{p}\boldsymbol{q}}:[0,1]\rightarrow\mathbb{R}$ is a
convex function, where
$\gamma_{\boldsymbol{p}\boldsymbol{q}}:[0,1]\rightarrow\mathcal{M}$ is a
geodesic that connects $\boldsymbol{p},\boldsymbol{q}$. $f$ is said to be a
geodesically concave function if $-f$ is a geodesically convex function.

######  Theorem 3.6.

[30] Let $\boldsymbol{A}\subseteq\mathcal{M}$ be a geodesically convex set. A
function $f:\boldsymbol{A}\rightarrow\mathbb{R}$ is geodesically convex if and
only if its epigraph
$\text{epi}(f)=\\{(\boldsymbol{p},c)|f(\boldsymbol{p})\leq
c\\}\subset\boldsymbol{A}\times\mathbb{R}$ is a convex set.

Now we present the main theoretical result of this paper.

######  Theorem 3.7.

For a given training data sample $\\{\boldsymbol{x},y\\}\in
S,0<\lVert\boldsymbol{x}\rVert<R<1$ (the hyperbolic feature $\boldsymbol{x}$
neither lie on the center nor lie on the boundary),
$l(\mu,\boldsymbol{\omega},b;\boldsymbol{x},y)$ is a geodesically convex
function on $\mathbb{R}^{+}\times\boldsymbol{A}\times\mathbb{R}^{+}$ and is a
geodesically concave function on
$\mathbb{R}^{+}\times\boldsymbol{B}\times\mathbb{R}^{+}$, where

$\boldsymbol{A}=\left\\{\boldsymbol{\nu}\in\mathbb{S}^{n-1}\bigg{|}y\cdot\frac{\boldsymbol{x}^{T}\boldsymbol{\nu}}{\lVert\boldsymbol{x}\rVert}>0\right\\}\subset\mathbb{S}^{n-1},\quad\boldsymbol{B}=\left\\{\boldsymbol{\nu}\in\mathbb{S}^{n-1}\bigg{|}y\cdot\frac{\boldsymbol{x}^{T}\boldsymbol{\nu}}{\lVert\boldsymbol{x}\rVert}<0\right\\}\subset\mathbb{S}^{n-1}.$
(8)

Note that both $\boldsymbol{A}$ and $\boldsymbol{B}$ are geodesically convex
sets.

###### Proof.

Let
$l(\mu,\boldsymbol{\omega},b;\boldsymbol{x},y)=\max(0,-g(\mu,\boldsymbol{\omega},b;\boldsymbol{x},y))$
where
$g(\mu,\boldsymbol{\omega},b;\boldsymbol{x},y)=y\cdot\left(\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b\right)$.
Since $\max(0,-a)$ is a convex function in $a\in\mathbb{R}$ and $g(\cdot)$ is
linear in $\mu$ and $b$, we only need to show that $g(\cdot)$, as a function
of $\boldsymbol{\omega}\in\mathbb{S}^{n-1}$, is geodesically convex (concave).
Without loss of generality, let $\mu=1$ and $b=0$. Also note that $y$ is the
label of data that takes values from $\\{-1,+1\\}$ which may flip the
inequality. It suffices to validate results for the positive sample, i.e,
$y=1$.

With a slight abuse of notation, let
$g(\boldsymbol{\omega};\boldsymbol{x})=g(1,\boldsymbol{\omega},0;\boldsymbol{x},1)=\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}=\ln\frac{1-\lVert\boldsymbol{x}\rVert^{2}}{\lVert\boldsymbol{\omega}-\boldsymbol{x}\rVert^{2}}$
defined on
$\boldsymbol{A}=\left\\{\boldsymbol{\nu}\in\mathbb{S}^{n-1}\big{|}\frac{\boldsymbol{x}^{T}\boldsymbol{\nu}}{\lVert\boldsymbol{x}\rVert}>0\right\\}\subset\mathbb{S}^{n-1}$.
Since $-\ln(\cdot)$ is decreasing and convex, we only need to check
$h(\boldsymbol{\omega};\boldsymbol{x})=\lVert\boldsymbol{\omega}-\boldsymbol{x}\rVert^{2}$
is geodesically convex on $\boldsymbol{A}$, i.e, check that the epigraph of
$h$ is a convex set. Note that
$\frac{\boldsymbol{x}^{T}\boldsymbol{\omega}}{\lVert\boldsymbol{x}\rVert}=\cos(\theta_{\boldsymbol{\omega}})$
for $\boldsymbol{\omega}\in\mathbb{S}^{n-1}$ where
$\theta_{\boldsymbol{\omega}}=\angle(\boldsymbol{\omega},\boldsymbol{x})$.

$\displaystyle\text{epi}(h)$
$\displaystyle=\\{(\boldsymbol{\omega},c)\in\boldsymbol{A}\times\mathbb{R}|\lVert\boldsymbol{\omega}-\boldsymbol{x}\rVert^{2}\leq
c\\}$ (9)
$\displaystyle=\left\\{(\boldsymbol{\omega},c)\bigg{|}\frac{\boldsymbol{x}^{T}\boldsymbol{\omega}}{\lVert\boldsymbol{x}\rVert}\geq\frac{1}{2\lVert\boldsymbol{x}\rVert}(1+\lVert\boldsymbol{x}\rVert^{2}-c)\right\\}$
$\displaystyle=\left\\{(\boldsymbol{\omega},c)|\cos(\theta_{\boldsymbol{\omega}})\geq
d(c)\right\\}=\begin{cases}\boldsymbol{A}\times[d(c),\infty)&\text{if
}d(c)\leq 0\\\ \boldsymbol{A}_{d}\times[d(c),\infty)&\text{if
}d(c)>0\end{cases}$

where $d(c)$ is a real number depending on $c$ and
$\lVert\boldsymbol{x}\rVert$ and
$\boldsymbol{A}_{d}=\left\\{(\boldsymbol{\omega},c)|\cos(\theta_{\boldsymbol{\omega}})\geq
d(c)\right\\}$ is the collection of unit vectors where the angle between the
vectors given the data $\boldsymbol{x}$ is small i.e., restricted to a small
region on the sphere. The last equality follows from the definition of
$\boldsymbol{A}$ and $\boldsymbol{A}_{d}$: if $d(c)\leq 0$, then
$\\{\omega:\cos(\theta_{\boldsymbol{\omega}})\geq
d(c)\\}\cap\boldsymbol{A}=\boldsymbol{A}$. Similarly, if $d(c)>0$,
$\\{\omega:\cos(\theta_{\boldsymbol{\omega}})\geq
d(c)\\}\cap\boldsymbol{A}=\\{\boldsymbol{\omega}\in\mathbb{S}^{n-1}|\cos(\theta_{\boldsymbol{\omega}})\geq
d(c)\\}\coloneqq\boldsymbol{A}_{d(c)}$. Both $\boldsymbol{A}$ and
$\boldsymbol{A}_{d}$ are geodesically convex sets and this completes the
proof. ∎

The convex sets $\boldsymbol{A},\boldsymbol{B}$ are the hemispheres of
$\mathbb{S}^{n-1}$ separated by the hyperplane
$\\{\boldsymbol{z}\in\mathbb{R}^{n}|\boldsymbol{x}^{T}\boldsymbol{z}=0\\}$
(the hyperplane has $\boldsymbol{x}$ as its normal vector) in its ambient
space $\mathbb{R}^{n}$. Theorem 3.7 tells us that given one data sample
$(\boldsymbol{x},y)$, the optimal value described in Eq. (7) exists in
$\mathbb{R}^{+}\times\boldsymbol{A}\times\mathbb{R}^{+}$, and it is globally
optimal. For a collection of training samples
$S=\\{(\boldsymbol{x}_{i},y_{i})\\}_{i=1}^{N}$, let
$\boldsymbol{A}_{i}=\left\\{\boldsymbol{\nu}\in\mathbb{S}^{n-1}\big{|}y_{i}\cdot\frac{\boldsymbol{x}_{i}^{T}\boldsymbol{\nu}}{\lVert\boldsymbol{x}_{i}\rVert}>0\right\\}$,
if the data are separable by a horosphere, it indicates that
$\cap_{i=1}^{N}\boldsymbol{A}_{i}\neq\emptyset$, then the global optimal can
be obtained using any gradient-based optimization on the Riemannian manifold
$\mathbb{R}^{+}\times\cap_{i=1}^{N}\boldsymbol{A}_{i}\times\mathbb{R}^{+}$.
Numerically, we apply a Riemannian gradient descent method on the entire space
$\mathbb{R}^{+}\times\mathbb{S}^{n-1}\times\mathbb{R}^{+}$ since
$g(\mu,\boldsymbol{\omega},b;\boldsymbol{x},y)$ is continuous.

###  3.4 Horospherical SVM

Given a horospherical decision boundary $\pi_{\mu,\boldsymbol{\omega},b}$
parameterized by
$\boldsymbol{\omega}\in\mathbb{S}^{n-1},\mu\in\mathbb{R}^{+}$, and
$b\in\mathbb{R}$, the margin $\gamma$ is the minimal distance from training
samples $S$ to the decision boundary:

$\gamma(\mu,\boldsymbol{\omega},b)=\inf_{\\{\boldsymbol{x},y\\}\in S}y\cdot
f(\boldsymbol{x};\mu,\boldsymbol{\omega},b)\cdot
d_{\mathbb{B}}(\boldsymbol{x},\pi_{\mu,\boldsymbol{\omega},b})=\inf_{\\{\boldsymbol{x},y\\}\in
S}y\cdot\frac{(\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b)}{\mu}.$
(10)

The maximum margin classifier can be obtained by solving the following
optimization problem:

$\displaystyle\max_{\mu,\boldsymbol{\omega},b}\quad$
$\displaystyle\gamma(\mu,\boldsymbol{\omega},b)\qquad s.t.\quad
y\cdot\frac{(\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b)}{\mu}\geq\gamma\quad\text{for
all }(\boldsymbol{x},y)\in S.$ (11)

######  Theorem 3.8.

The maximum margin classification problem in hyperbolic space with horosphere
as its decision boundary described in Eq. (11) is equivalent to the following
optimization problem:

$\min_{\mu,\boldsymbol{\omega},b}\quad\frac{1}{2}\mu^{2}\qquad s.t.\quad
y\cdot(\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b)\geq
1\quad\text{for all }(\boldsymbol{x},y)\in S.$ (12)

The proof is analogous to that in the Euclidean case [3]. Note that the margin
is unchanged if we apply the following scale transformation:
$\mu\rightarrow\mu/\gamma$ and $b\rightarrow b/\gamma$. We can also build a
soft-margin horospherical SVM, dubbed HoroSVM, by minimizing the following
loss function:

$l(\mu,\boldsymbol{\omega},b;\boldsymbol{x},y)=\frac{1}{2}\mu^{2}+C\sum\limits_{i=1}^{\lvert
S\rvert}\max(0,1-y_{i}\cdot(\mu\langle\boldsymbol{\omega},\boldsymbol{x}_{i}\rangle_{\mathbb{B}}-b)).$
(13)

where $C$ is a hyperparameter that controls the tradeoff between minimizing
misclassification and maximizing margin.

It is easy to see that the same result as in Theorem 3.7 holds for the loss
function in Eq. (13). That is, the loss function is a geodesically convex
function on one geodesically convex subset of the parameter space and a
geodesically concave function on the other geodesically convex subset of the
parameter space. Recall that the idea behind proving that the loss function in
the horospherical perceptron, $\max(0,-g(\cdot))$, where
$g(\cdot)=y\cdot(\mu\langle\boldsymbol{\omega},\boldsymbol{x}\rangle_{\mathbb{B}}-b)$,
is geodesically convex is based on the important fact that $\max(0,-a)$ is a
convex function in $a\in\mathbb{R}$. Similarly, the same idea applies to
HoroSVM, where the hinge loss is used, and the loss function becomes
$\frac{1}{2}\mu^{2}+\max(0,1-g(\cdot))$. Note that $\max(0,1-a)$ is also a
convex function in $a\in\mathbb{R}$ and $\frac{1}{2}\mu^{2}$ is a convex
function in $\mu$, these facts complete the proof of the desired property for
the HoroSVM, which is analogous to Theorem 3.7 for horospherical perceptron.

We can then apply any Riemannian gradient descent optimization methods for
updating the parameters in HoroSVM since the problem is an optimization
problem over a product space of Riemannian manifolds,
$\mathbb{R}^{+}\times\mathbb{S}^{n-1}\times\mathbb{R}^{+}$. We refer the
readers to [5] and [1] for more details about optimization techniques on
Riemannian manifolds.

##  4 Experiments

In this section, we present several experimental results obtained from an
application of our HoroSVM to synthetic data as well as real data sets used in
published literature. Our implementation is based on Pymanopt [19] using the
Riemannian conjugate gradient method [26] on Intel(R) Xeon(R) CPU E5-2683 v3 @
2.00GHz.

###  4.1 Network Data Set

Here, we follow the experimental setup in [12], and evaluate our HoroSVM over
four real-world network data sets used by [8]: karate [35] (2 classes, 34
nodes ), polblogs [2] (2 classes, 1224 nodes ), polbooks 222http://www-
personal.umich.edu/~mejn/netdata/ (3 classes, 105 nodes ), and football [17]
(12 classes, 115 nodes ).

The network data is embedded in a 2D hyperbolic space. Given the hyperbolic
embeddings, we compare our HoroSVM with three other competing large margin
classifiers: Euclidean SVM (even though it violates the hyperbolic geometry),
hyperboloid SVM [12], and Poincaré SVM [11]. A one-verses-rest strategy is
applied for multiclass classification. We conducted a five-fold cross-
validation on each data set, where we chose the hyperparameter $C$ from
$\\{1,5,10\\}$ during the cross-validation procedure. The mean of the F1 score
followed by the standard deviation over five trials are summarized in Table 1.
As evident from the table, our method yields the best results on all the data
sets.

HoroSVM outperformed other methods on all four data sets. The data in karate
are well-separated and thus both Euclidean SVM and hyperboloid SVM performed
equally well. Our method outperforms the others since the horospheres have
several nice properties, the most important of which is that the Busemann
function whose level sets are the horospheres is a convex function that
guarantees global optimality in the optimization. Notice that the performance
of Poincaré SVM on karate is inferior to others by a significant amount. The
reason is that the performance of Poincaré SVM is sensitive to the choice of
the reference point. We demonstrate our performance gain on the remaining data
sets, and our method is more consistent, compared to Euclidean SVM and
hyperboloid SVM, in terms of lower standard deviation, on football data set
where data exhibit a larger variance/spread.

Table 1: F1 scores for node classification on network datasets. Boldface indicates best performance. Methods | Karate | Polblogs | Polbooks | Football
---|---|---|---|---
Euclidean SVM | 0.95 $\pm$ 0.06 | 0.92 $\pm$ 0.02 | 0.83 $\pm$ 0.03 | 0.29 $\pm$ 0.12
Hyperboloid SVM | 0.95 $\pm$ 0.06 | 0.92 $\pm$ 0.01 | 0.83 $\pm$ 0.03 | 0.30 $\pm$ 0.14
Poincaré SVM | 0.78 $\pm$ 0.16 | 0.92 $\pm$ 0.02 | 0.84 $\pm$ 0.03 | 0.32 $\pm$ 0.04
HoroSVM (Ours) | 0.98 $\pm$ 0.04 | 0.93 $\pm$ 0.01 | 0.85 $\pm$ 0.04 | 0.34 $\pm$ 0.06
Table 2: F1 scores for subtree classification on four subtrees of WordNet. Boldface indicates best performance on 2D embeddings of each dataset.  Methods | animal.n.01 | group.n.01 | worker.n.01 | mammal.n.01
---|---|---|---|---
3218/798 | 6649/1727 | 861/254 | 953/228
Hyperboloid SVM (D = 2) | 0.53 $\pm$ 0.07 | 0.52 $\pm$ 0.01 | 0.54 $\pm$ 0.04 | 0.39 $\pm$ 0.03
Hyperbolic LR (D = 2) | 0.46 $\pm$ 0.08 | 0.52 $\pm$ 0.04 | 0.54 $\pm$ 0.07 | 0.32 $\pm$ 0.10
Hyperbolic LR (D = 5) | 0.95 $\pm$ 0.03 | 0.76 $\pm$ 0.07 | 0.80 $\pm$ 0.08 | 0.78 $\pm$ 0.04
Hyperbolic LR (D = 10) | 0.96 $\pm$ 0.01 | 0.86 $\pm$ 0.05 | 0.84 $\pm$ 0.04 | 0.94 $\pm$ 0.04
Euclidean SVM (D = 2) | 0.39 $\pm$ 0.01 | 0.39 $\pm$ 0.00 | 0.32 $\pm$ 0.02 | 0.20 $\pm$ 0.01
Euclidean SVM (D = 5) | 0.95 $\pm$ 0.00 | 0.79 $\pm$ 0.01 | 0.38 $\pm$ 0.02 | 0.44 $\pm$ 0.01
Euclidean SVM (D = 10) | 0.97 $\pm$ 0.00 | 0.91 $\pm$ 0.00 | 0.46 $\pm$ 0.04 | 0.72 $\pm$ 0.05
HoroSVM (D = 2) | 0.57 $\pm$ 0.07 | 0.65 $\pm$ 0.01 | 0.62 $\pm$ 0.01 | 0.42 $\pm$ 0.01
HoroSVM (D = 5) | 0.93 $\pm$ 0.01 | 0.88 $\pm$ 0.00 | 0.82 $\pm$ 0.04 | 0.88 $\pm$ 0.01
HoroSVM (D = 10) | 0.95 $\pm$ 0.02 | 0.91 $\pm$ 0.01 | 0.86 $\pm$ 0.01 | 0.93 $\pm$ 0.02

###  4.2 Subtree Classification in WordNet

A task of considerable interest in hyperbolic space classification problems is
to determine whether a node belongs to a given subtree in the hyperbolic
embedding. We obtained hyperbolic embeddings in various dimensions using the
approach in [15] for WordNet[31] noun hierarchy (82,115 nodes). We consider
four subtrees whose roots are the following synsets: ANIMAL.N.01, GROUP.N.01,
WORKER.N.01, and MAMMAL.N.01.

We split all nodes in a subtree into positive training (80%) and test (20%)
nodes and applied the same process to the remaining WordNet nodes to create
negative training and test sets. The average F1 scores and the standard
deviations over 3 trials are shown in Table 2. The number of positive
training/test samples of each data set are listed as well. We exclude Poincaré
SVM from the comparisons in this task. The reason being, data are highly
imbalanced in this task and the positive samples are clustered near the
boundary. The reference point learned in Poincaré SVM will be close to the
boundary where the tangent approximation of data at this reference point is
highly distorted, as opposed to the original hyperbolic embeddings. It is
therefore hard to locate a hyperplane in the tangent space that separates the
lifted (mapped) data. In addition, the learning of the reference point is only
applicable to 2D hyperbolic space.

As well known in the Euclidean SVM literature, a vanilla (unweighted)
implementation of SVM performs poorly on extremely imbalanced data, we
observed the same behavior in training HoroSVM on this task. We preprocessed
the data by downsampling the majority class of samples (the negative samples)
to train a robust model. Since there is no protocol for dealing with
unbalanced data in training a hyperboloid SVM, we presented the Euclidean SVM
results using the same preprocessed data for reference. In addition, we
presented results of hyperbolic logistic regression (LR) [15], which is not a
large-margin classifier on this task, where the imbalanced data is handled by
sampling the equal number of negative and positive nodes in each mini-batch of
size 16 during training.

Now, we highlight several results in Table 2. The superior performance of
hyperboloid SVM over hyperbolic LR is expected, as both methods use geodesic
decision boundaries but hyperboloid SVM aims to maximize the margin. However,
the training of hyperboloid SVM is highly unstable as we mentioned earlier due
to the non-convex optimization process. Euclidean SVM under-performs as it
does not take into account the hyperbolic geometry. Our HoroSVM exhibits a
significant improvement in predicting words in a subtree, as evidenced by
higher F1 scores across all the subtrees. The small number of nodes within a
subtree, compared to the whole WordNet, causes the nodes to cluster near the
boundary in their hyperbolic embeddings. Thus, horosphere is an ideally suited
decision boundary (in comparison to the geodesic boundary in [12]) to isolate
the subtree.

###  4.3 Synthetic Data with Noisy Labels

Figure 4: Training (left) and test (right) F1 scores of several methods on
synthetic data with noisy labels at different noise levels.

To demonstrate the robustness of our HoroSVM, we apply it to synthetic data
with noisy labels at varying levels/amounts of noise. Specifically, we
generated 100 synthetic datasets by sampling from a Gaussian mixture model
defined on the Poincaré disk model as in [12]. The isotropic Gaussian
distribution in hyperbolic space is referred to as the _Riemannian normal_
distribution, and we used the sampling method presented in [21]. For each
dataset, we sampled two centroids from a zero-mean Riemannian normal
distribution with a variance of 1.5. We then sampled 200 data points from a
unit-variate Riemannian normal distribution centered at each centroid,
resulting in a dataset of 400 points classified into positive and negative
classes. We split the dataset into training and test sets with 100
positive/negative samples in the training set and 100 positive/negative
samples in the test set. We then generated datasets with noisy labels at noise
levels: $\eta\in\\{0,0.05,0.1,\dots,0.5\\}$ by flipping the labels of a
proportion $\eta$ of the training (not test) samples, with an equal number of
positive and negative samples flipped. The train/test average F1 scores of
each method across all datasets at varying noise levels are shown in Fig 4. We
compared our HoroSVM with hyperboloid SVM, hyperbolic LR, and a two-layer HNN
[15] (with a hidden dimension of 5). While all methods depict decreasing
training F1 scores as the noise level increases, HoroSVM outperforms the
others consistently throughout the training process. Hyperbolic LR and HNN
exhibit the least resistance to label noise, with test F1 scores dropping
(faster) with increasing noise level. Both hyperboloid SVM and HoroSVM
demonstrate consistent performance across different noise levels, owing to the
inherent robustness of large-margin classifiers. However, the training of
hyperboloid SVM is highly unstable resulting its inferior performance. HoroSVM
demonstrate its superiority in accuracy and robustness as evidenced in the
results.

##  5 Discussion and Conclusions

In this paper, we presented a novel large margin classifier, dubbed HoroSVM,
whose decision boundaries are horospheres that are the level sets of a
Busemann function. We presented a novel formulation leading to the
optimization of a geodesically convex loss performed using a Riemannian
gradient-based method and guaranteeing a globally optimal solution. We
demonstrated superior to competitive performance of the HoroSVM over SOTA
large margin classifiers.

In Euclidean space, a kernel SVM is usually favored over the linear SVM due to
its ability to cope with non-linearly separable data. In Hyperbolic space, the
challenge lies in developing valid positive definite kernels (see [14] for
details on validity of kernels on Riemannian manifolds). The only reported
work on KSVM in hyperbolic space that we are aware of is [12], which uses a
kernel that violates the positive definiteness property of RKHS kernels. Thus
the problem of interest is primarily defining a valid family of kernels in
hyperbolic space. We will address this in our future work.

## Acknowledgements

This research was in part funded by the NSF grant IIS 1724174 and the NIH
NINDS and NIA grant RF1NS121099 to Vemuri.

## References

  * [1] P-A Absil, Robert Mahony, and Rodolphe Sepulchre.  Optimization algorithms on matrix manifolds.  In Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2009.
  * [2] Lada A Adamic and Natalie Glance.  The political blogosphere and the 2004 us election: divided they blog.  In Proceedings of the 3rd international workshop on Link discovery, pages 36–43, 2005.
  * [3] Christopher M Bishop and Nasser M Nasrabadi.  Pattern recognition and machine learning, volume 4.  Springer, 2006.
  * [4] Roberto Bonola.  Non-Euclidean geometry: A critical and historical study of its development.  Courier Corporation, 1955.
  * [5] Nicolas Boumal.  An introduction to optimization on smooth manifolds.  Available online, May, 3, 2020.
  * [6] Martin R Bridson and André Haefliger.  Metric spaces of non-positive curvature, volume 319.  Springer Science & Business Media, 2013.
  * [7] James W. Cannon, William J. Floyd, Richard Kenyon, and Walter R. Parry.  Hyperbolic Geometry, volume 31.  MSRI Publications, 1997.
  * [8] Benjamin Paul Chamberlain, James Clough, and Marc Peter Deisenroth.  Neural embeddings of graphs in hyperbolic space.  arXiv preprint arXiv:1705.10359, 2017.
  * [9] Ines Chami, Albert Gu, Dat P Nguyen, and Christopher Ré.  Horopca: Hyperbolic dimensionality reduction via horospherical projections.  In International Conference on Machine Learning, pages 1419–1429. PMLR, 2021.
  * [10] Ines Chami, Zhitao Ying, Christopher Ré, and Jure Leskovec.  Hyperbolic graph convolutional neural networks.  Advances in Neural Information Processing Systems, 32:4868–4879, 2019.
  * [11] Eli Chien, Chao Pan, Puoya Tabaghi, and Olgica Milenkovic.  Highly scalable and provably accurate classification in poincaré balls.  In 2021 IEEE International Conference on Data Mining (ICDM), pages 61–70. IEEE, 2021.
  * [12] Hyunghoon Cho, Benjamin DeMeo, Jian Peng, and Bonnie Berger.  Large-margin classification in hyperbolic space.  In The 22nd international conference on artificial intelligence and statistics, pages 1832–1840. PMLR, 2019.
  * [13] Xiran Fan, Chun-Hao Yang, and Baba C Vemuri.  Nested hyperbolic spaces for dimensionality reduction and hyperbolic nn design.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 356–365, 2022.
  * [14] Aasa Feragen, Francois Lauze, and Soren Hauberg.  Geodesic exponential kernels: When curvature and linearity conflict.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3032–3042, 2015.
  * [15] Octavian Ganea, Gary Bécigneul, and Thomas Hofmann.  Hyperbolic entailment cones for learning hierarchical embeddings.  In International Conference on Machine Learning, pages 1646–1655. PMLR, 2018.
  * [16] Octavian-Eugen Ganea, Gary Bécigneul, and Thomas Hofmann.  Hyperbolic neural networks.  Advances in Neural Information Processing Systems 31, pages 5345–5355, 2019.
  * [17] Michelle Girvan and Mark EJ Newman.  Community structure in social and biological networks.  Proceedings of the national academy of sciences, 99(12):7821–7826, 2002.
  * [18] Sigurdur Helgason.  Groups and geometric analysis: integral geometry, invariant differential operators, and spherical functions, volume 83.  American Mathematical Society, 2022.
  * [19] Niklas Koep and Sebastian Weichwald.  Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation.  Journal of Machine Learning Research, 17:1–5, 2016.
  * [20] Gian Marconi, Carlo Ciliberto, and Lorenzo Rosasco.  Hyperbolic manifold regression.  In International Conference on Artificial Intelligence and Statistics, pages 2570–2580. PMLR, 2020.
  * [21] Emile Mathieu.  Charline le lan, chris j maddison, ryota tomioka, and yee whye teh. continuous hierarchical representations with poincaré variational auto-encoders.  Advances in Neural Information Processing Systems, pages 12544–12555, 2019.
  * [22] Nicholas Monath, Manzil Zaheer, Daniel Silva, Andrew McCallum, and Amr Ahmed.  Gradient-based hierarchical clustering using continuous representations of trees in hyperbolic space.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 714–722, 2019.
  * [23] Maximillian Nickel and Douwe Kiela.  Poincaré embeddings for learning hierarchical representations.  Advances in neural information processing systems, 30:6338–6347, 2017.
  * [24] Frederic Sala, Chris De Sa, Albert Gu, and Christopher Ré.  Representation tradeoffs for hyperbolic embeddings.  In International conference on machine learning, pages 4460–4469. PMLR, 2018.
  * [25] Rik Sarkar.  Low distortion delaunay embedding of trees in hyperbolic plane.  In International Symposium on Graph Drawing, pages 355–366. Springer, 2011.
  * [26] Hiroyuki Sato.  Riemannian conjugate gradient methods: General framework and specific algorithms with convergence analyses.  SIAM Journal on Optimization, 32(4):2690–2717, 2022.
  * [27] Ryohei Shimizu, YUSUKE Mukuta, and Tatsuya Harada.  Hyperbolic neural networks++.  In International Conference on Learning Representations, 2020.
  * [28] Sho Sonoda, Isao Ishikawa, and Masahiro Ikeda.  Fully-connected network on noncompact symmetric space and ridgelet transform based on helgason-fourier analysis.  In International Conference on Machine Learning, pages 20405–20422. PMLR, 2022.
  * [29] Alexandru Tifrea, Gary Becigneul, and Octavian-Eugen Ganea.  Poincare glove: Hyperbolic word embeddings.  In International Conference on Learning Representations, 2018.
  * [30] Constantin Udriste.  Convex functions and optimization methods on Riemannian manifolds, volume 297.  Springer Science & Business Media, 2013.
  * [31] Pinceton University.  About wordnet, 2010.
  * [32] Ming-Xi Wang.  Laplacian eigenspaces, horocycles and neuron models on hyperbolic spaces.  2021\.
  * [33] Melanie Weber, Manzil Zaheer, Ankit Singh Rawat, Aditya K Menon, and Sanjiv Kumar.  Robust large-margin learning in hyperbolic space.  Advances in Neural Information Processing Systems, 33:17863–17873, 2020.
  * [34] Tao Yu and Christopher De Sa.  Hyla: Hyperbolic laplacian features for graph learning.  arXiv preprint arXiv:2202.06854, 2022.
  * [35] Wayne W Zachary.  An information flow model for conflict and fission in small groups.  Journal of anthropological research, 33(4):452–473, 1977.
  * [36] Hongyi Zhang and Suvrit Sra.  First-order methods for geodesically convex optimization.  In Conference on Learning Theory, pages 1617–1638. PMLR, 2016.