# Angle estimation using mmWave RSS measurements with enhanced multipath
information

Neharika Valecha1, Jesus Omar Lacruz2, Michael Lentmaier1, Joerg Widmer2, and
Fredrik Tufvesson1 This project has received funding from the European Union’s
EU Framework Programme for Research and Innovation Horizon 2020 under Grant
Agreement No. 861222. 1 Department of Electrical and Information Technology,
Lund University, Lund, Sweden  2 IMDEA Networks Institute, Madrid, Spain

###### Abstract

mmWave communication has come up as the unexplored spectrum for 5G services.
With new standards for 5G NR positioning, more off-the-shelf platforms and
algorithms are needed to perform indoor positioning. An object can be
accurately positioned in a room either by using an angle and a delay estimate
or two angle estimates or three delay estimates. We propose an algorithm to
jointly estimate the angle of arrival (AoA) and angle of departure (AoD),
based only on the received signal strength (RSS). We use mm-FLEX, an
experimentation platform developed by IMDEA Networks Institute that can
perform real-time signal processing for experimental validation of our
proposed algorithm. Codebook-based beampatterns are used with a uniquely
placed multi-antenna array setup to enhance the reception of multipath
components and we obtain an AoA estimate per receiver thereby overcoming the
line-of-sight (LoS) limitation of RSS-based localization systems. We further
validate the results from measurements by emulating the setup with a simple
ray-tracing approach.

##  I Introduction

The ever-increasing need for larger bandwidths and higher data rates has led
the telecom industry to move to hitherto unexplored higher frequency
millimeter-wave (mmWave) band (20-300 GHz). This range of frequencies has
quasi-optical properties that make it ideal for positioning applications. From
the Internet of Things (IoT) to Vehicle-to-everything (V2X), most emerging
applications need location awareness [1]. In the past, Global Navigation
Satellite Systems have been used for UE positioning, but that is preferable
for outdoor systems. When trying to position UEs indoors, range-based or
range-free algorithms are used.

mmWave spectrum frequencies (over 24 GHz), due to their shorter wavelengths
have longer, narrower beams [2] and are ideal for positioning. Consequently,
the higher path loss due to the shorter wavelengths makes the use of massive
antenna arrays mandatory. Furthermore, using beamforming in combination with a
large number of antennas leads to highly directional beams. Fully digital
beamforming, which supports massive MIMO arrays (upwards of $128\times 256$)
has emerged as a ’paradise for positioning’. With calibrated systems, it
offers the use of super-resolution algorithms that can give position estimates
with centimeter-level accuracy. However, the hardware and power constraints
hinder its commercialization. Thus, hybrid beamforming is a solution, as it
provides a trade-off between power consumption and position resolution [3].

Localization using radio frequencies has its own challenges. 5G systems
augment the use of massive antenna arrays which allow for increased
resolvability in the angular domain. Most systems rely strongly on line-of-
sight (LoS) information for positioning but these signal components are often
blocked in indoor environments. Range-based localization hinges on time
synchronization between the transmitter (TX) and receiver (RX) and requires
calibrated antennas and coherent communication between TX and RX. Multipath
components (MPCs) can be exploited for location information if effectively
resolved [4], [5], [6]. Received Signal Strength (RSS) based localization is
thus effective and can exploit hybrid beamforming setups for accurate
positioning [7].

In this paper, we propose an RSS-based least squares (RSS-LS2D) algorithm to
jointly estimate the Angle of Arrival (AoA) and the Angle of Departure (AoD)
by minimizing over both TX and RX beampatterns. Our algorithm is an extension
of the algorithm proposed in [8], referred to as RSS-LS1D henceforth. It is 1D
because either the AoA is estimated at the RX or the AoD at the TX, at one
time. The algorithms are applied to channel measurements obtained using the
mm-FLEX testbed described in [8]. As the RSS measurements rely on the strength
of the signal, we use a special receiver with 4 RX antenna arrays placed such
that they receive MPCs over a $360^{\circ}$ range. We show that our RSS-LS2D
algorithm estimates AoA/AoD accurately even in the absence of LoS. Further, we
validate our results by emulating the measurement environment using ray
tracing.

Figure 1: Left: The 4 RX antenna arrays are connected to the SIVERS FPGA to
perform real-time processing. Right: Room at IMDEA with TX and RX placed at
positions $\mathbf{p}_{1}=[0\,\mathrm{m},2.4\,\mathrm{m}]$ and
$\mathbf{p}_{\mathrm{RX}}=[0,0]$, respectively.

##  II Setup

###  II-A mm-FLEX Testbed

We use Millimeter Wave Mobile FulL-Bandwidth EXperimentation Platform (mm-
FLEX), an open-source research platform with a 60 GHz RF front end that
supports bandwidths up to 2 GHz. The TX and RX each have the AMC599 board,
integrating a Xilinx Kintex Ultrascale FPGA, highspeed AD/DA converters, and
DDR4 memory banks [9] and an AMC726 board with an Intel Core i7 processor
[10]. It uses the EVK06002 development kit from Sivers IMA [11] with 60 GHz
up/down converters as the RF front-end for mm-FLEX. The kit includes 2D analog
beam-forming capabilities through phase shifters with a 6-bit resolution for
each antenna element. In addition, the kit incorporates fast beam-switching
capabilities through simple pulses via a GPIO interface, that allows switching
beam patterns every $10\,\mathrm{ns}$. However, the kit has a maximum RF
settling time of $35\,\mathrm{ns}$ when changing beam patterns. The RX setup
is as shown in Fig. 1. The TX uses a single RF linear antenna array with 16
antennas and the receiver has 4 separate arrays of the same type. The 4 RX
arrays are placed orthogonally such that they receive multipath information
and provide a quasi-omnidirectional view of the room. Each TX and RX array has
64 codebook-based beampatterns with main lobes at angles ranging from -45∘ to
45∘, illustrated in Fig. 2. The TX array and each of the RX arrays are loaded
with the same codebook and have the same beampatterns, each covering an
angular span of approximately $160^{\circ}$. In particular, $\mathrm{RX}_{1}$
has an angular span from $-79.5^{\circ}$ to $80^{\circ}$, $\mathrm{RX}_{2}$
has an angular span from $10.5^{\circ}$ to $170^{\circ}$, $\mathrm{RX}_{3}$
has an angular span from $100.5^{\circ}$ to $260^{\circ}$ and
$\mathrm{RX}_{4}$ has an angular span from $190.5^{\circ}$ to $350^{\circ}$.
Thus, all 4 RX arrays together cover a $360^{\circ}$ range with some overlap.
The TX array has the same angular span as $\mathrm{RX}_{3}$.

Figure 2: Illustration of a selection of measured beam patterns used in mm-
FLEX.

###  II-B Simulation and Measurement Environment

We validate the proposed RSS-LS2D algorithm with real measured data as well as
synthetically generated data. The measurements were carried out at IMDEA, the
room with the TX and RX is shown in Fig. 1. The room is $5\,\mathrm{m}\times
8\,\mathrm{m}$. As described in [8], an IEEE 802.11 compatible frame is
transmitted with Golay (Ga128) sequences in the preamble and Zadoff-Chu
sequences as payload, having a bandwidth of 1.76 GHz. The hardware includes a
square-root-raised-cosine (SRRC) filter which was switched off during our
experiments. Instead, the frame was up/downsampled by using Fast Fourier
Transform (FFT) as in 4G LTE systems. The frame structure can be seen in Fig.
3.

Figure 3: Structure of the transmitted frame with Golay sequences of length
128 bits and Zadoff-Chu sequences of length 384 bits.

The preamble consists of the Short Training Field (STF), Channel Estimation
Field (CEF), and the C-PHY header containing information about the MCS, data
length, and some other data. The STF contains 25 repeated Golay (Gb128)
sequences of length 128 samples each and is needed for packet detection and
synchronization. The CEF are four pairs of complementary Golay sequences.
Zadoff-Chu sequences of length $K=384$ samples111 The size of the payload is
limited by the storage available in the hardware. As the channel processing
happens in real-time, the total data over all beampatterns should not exceed
the DDR4 memory bank capacity. are used as payload data, which can be
generated as

$z_{k}=e^{-j\pi\frac{uk(k+c+2q)}{K}}$

where $k=0,1,\ldots,K-1$, $u=0,1,\ldots,K-1$ and $\gcd(u,K)=1$,
$c=\mathrm{mod}(K,2)$ and $q\in\mathbb{Z}$. The frame is calibrated to receive
the same Zadoff-Chu sequence over the $N_{\mathrm{BP},\mathrm{RX}}$
beampatterns. We obtain this frame at each of the four RX. The transmission
and reception takes place over all $N_{\mathrm{BP},\mathrm{TX}}\times
N_{\mathrm{BP},\mathrm{RX}}$ beampatterns, unlike in [8] where a beam
refinement protocol (BRP) was used to find aligning beampatterns and the frame
was transmitted over only the strongest TX beam pattern.

The position of the RX is stationary and known, while the TX is moved to
different positions as shown in Fig. 4. We place the receiver’s ‘center’ at
the origin with each of the 4 receivers placed at
$\mathbf{p}_{\mathrm{RX},1}=[0.25\,\mathrm{m},0\,\mathrm{m}]$,
$\mathbf{p}_{\mathrm{RX},2}=[0\,\mathrm{m},0.25\,\mathrm{m}]$,
$\mathbf{p}_{\mathrm{RX},3}=[-0.25\,\mathrm{m},0\,\mathrm{m}]$,
$\mathbf{p}_{\mathrm{RX},4}=[0\,\mathrm{m},-0.25\,\mathrm{m}]$. The TX is
placed at 10 different positions,
$\mathbf{p}_{1}=[2.4\,\mathrm{m},0\,\mathrm{m}]$,
$\mathbf{p}_{2}=[2.4\,\mathrm{m},0.6\,\mathrm{m}]$,
$\mathbf{p}_{3}=[2.4\,\mathrm{m},1.2\,\mathrm{m}]$,
$\mathbf{p}_{4}=[2.4\,\mathrm{m},-0.3\,\mathrm{m}]$,
$\mathbf{p}_{5}=[2.4\,\mathrm{m},-0.6\,\mathrm{m}]$,
$\mathbf{p}_{6}=[3\,\mathrm{m},0\,\mathrm{m}]$,
$\mathbf{p}_{7}=[3\,\mathrm{m},0.6\,\mathrm{m}]$,
$\mathbf{p}_{8}=[3\,\mathrm{m},1.2\,\mathrm{m}]$,
$\mathbf{p}_{9}=[3\,\mathrm{m},-0.3\,\mathrm{m}]$,
$\mathbf{p}_{10}=[3\,\mathrm{m},-0.6\,\mathrm{m}]$. The positions were chosen
to study the effect on the received signal laterally as well as
longitudinally. The receiver has 4 antenna arrays placed orthogonally, while
the transmitter faces one of the receiver arrays as shown in Fig. 4. Thus, we
ensure LoS between the TX and $\mathrm{RX}_{1}$ but also nonline-of-sight
(NLoS) on the other three RX. We use the measured RSS to estimate the AoA/AoD
in azimuth only. The angles are measured w.r.t the plane parallel to the
ground and counter-clockwise for each receiver. The synthetic environment is
based on the measurement environment.

Figure 4: Floor plan of the room with RX placed at the origin,
$\mathbf{p}_{1}$ to $\mathbf{p}_{10}$ represent the different TX static
positions.

##  III System Model

Consider an indoor office scenario with a TX equipped with a Uniform Linear
Array(ULA) and a quasi-omnidirectional RX connected to 4 identical ULAs placed
perpendicular to each other. In terms of RF chains, this forms a SIMO setup,
which is chosen to enhance the number of MPCs. The transmitted signal
$x[k]=[x_{1},x_{2},\ldots,x_{K}]$ is a time-varying, repeating Zadoff-Chu
sequence as mentioned above. The received signal measured at
$\mathrm{RX}_{i}$, $i=1,\dots,4$, is

$y_{i,t,r}[k]=\sum_{l=0}^{L-1}a_{l}b_{r}(\theta_{l})b_{t}(\phi_{l})x[k-\tau_{l}]+\eta[k]\
,$ (1)

where $L$ is the number of MPCs, $a_{l}$ is the path loss with amplitude
$\alpha$ and phase $\varphi\in[0,2\pi]$ and $b_{r}(\theta_{l})$ and
$b_{t}(\phi_{l})$ are beampattern weights for receiver and transmitter,
respectively, with AoA $\theta$ and AoD $\phi\in[0,360^{\circ}]$. $\tau_{l}$
is the time delay and $\eta$ is the noise with distribution
$\mathcal{N}(x;\mu,\mathrm{N_{0}})$, $t$ and $r$ correspond to the beampattern
numbers at the TX and RX, respectively. As can be observed from Fig. 2, each
beampattern contributes to every sample of $\theta$ and/or $\phi$.

##  IV The RSS-based Angle Estimation

We first compute the RSS at each $\mathrm{RX}_{i}$ as,

$\mathrm{RSS}^{(i)}_{t,r}=\sum_{k=1}^{K}y[k]y[k]^{*}\ .$

A first rough estimate of the AoA/AoD can be obtained by observing the RSS in
the beampattern domain. Since the beampatterns are packed quite closely to
each other, merely observing the power received across beampatterns provides
us with an estimate for the general direction of the strongest beam of the
signal. AoA can be assigned to the dominating angle of the beampattern with
the maximum gain.

Figure 5: RSS plot for $\mathrm{RX}_{1}$ at position $\mathbf{p}_{5}$ showing
received power across TX and RX beampatterns. Here we see the strongest power
between beampatterns numbered $t=15-25$ and $r=10-20$.

(a) Receiver 1 (Measured) (b) Receiver 1 (Simulated) (c) Receiver 2 (Measured)
(d) Receiver 2 (Simulated)

(e) Receiver 3 (Measured) (f) Receiver 3 (Simulated) (g) Receiver 4 (Measured)
(h) Receiver 4 (Simulated)

Figure 6: Measured RSS data and simulated RSS data for all RX when TX is at
position $\mathbf{p}_{1}$

(a) Receiver 1 (Measured) (b) Receiver 1 (Simulated) (c) Receiver 2 (Measured)
(d) Receiver 2 (Simulated)

(e) Receiver 3 (Measured) (f) Receiver 3 (Simulated) (g) Receiver 4 (Measured)
(h) Receiver 4 (Simulated)

Figure 7: Measured RSS data and simulated RSS data for all RX when TX is at
position $\mathbf{p}_{5}$

An example of RSS is shown in Fig. 5. Observing position $\mathbf{p}_{5}$ in
Fig. 4, we see a LoS path between $\mathrm{TX}$ and $\mathrm{RX}_{1}$. The
beampatterns move counterclockwise from 0 to 64 for both TX and RX, the
aligning beampatterns should then be approx. 10-15 at TX and 15-20 at RX. This
is seen in Fig. 5, where the beampatterns with the strongest received power
are between 10-20. We pick beampattern 14 with maximum RSS for TX and
beampattern 19 at RX. Now, we find the angle where the beampattern has its
main lobe and assign that as AoD/AoA which gives AoD as $-15^{\circ}$ and AoA
as $-23.5^{\circ}$, respectively.

###  IV-A RSS-LS2D Algorithm

Let us now build a concrete algorithm that extends the single-angle estimation
algorithm proposed in [8] to joint AoA/AoD estimation. Since our measurements
bypass BRP, we obtain a full $N_{\mathrm{\mathrm{BP},{TX}}}\times
N_{\mathrm{\mathrm{BP},{RX}}}$ beampattern grid across the TX and each of the
four RX arrays. So, instead of just identifying one dominant beampattern, we
formulate a least-squares problem utilizing the beampattern information at
both the TX and RX. The joint AoD and AoA estimation can then be written as:

$\displaystyle[\hat{\theta},\hat{\phi}]=$
$\displaystyle\operatorname*{arg\,min}_{\theta,\phi}\min_{\alpha}\sum_{t,r}{}||{\mathrm{RSS}^{(i)}_{t,r}-\alpha
b_{r}(\theta)b_{t}(\phi)}||^{2}$ (2) $\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{\theta,\phi}\min_{\alpha}\Big{\\{}|\alpha|^{2}\sum_{t,r}||b_{r}(\theta)b_{t}(\phi)||^{2}$
(3)
$\displaystyle-2\alpha\sum_{t,r}\mathrm{RSS}^{(i)}_{t,r}b_{r}(\theta)b_{t}(\phi)\Big{\\}}$

where $\alpha$ is the received signal amplitude for a pair of given beam
patterns. Taking derivatives we find the $\alpha^{*}$ that minimizes the norm:

$D_{1}(\alpha)=2\alpha\Big{\\{}\sum_{t,r}||b_{r}(\theta)b_{t}(\phi)||^{2}\Big{\\}}-2\sum_{t,r}\mathrm{RSS}^{(i)}_{t,r}b_{r}(\theta)b_{t}(\phi)$

$D_{2}(\alpha)=2\Big{\\{}\sum_{t,r}||b_{r}(\theta)b_{t}(\phi)||^{2}\Big{\\}}\geq
0\ .$

Here $D_{1}(\alpha)$ and $D_{2}(\alpha)$ are the first and second-order
derivatives taken w.r.t $\alpha$. Setting $D_{1}(\alpha)=0$ and substituting
$\alpha^{*}$ in (3), we obtain

$[\hat{\theta},\hat{\phi}]=\operatorname*{arg\,min}_{\theta,\phi}\left\\{-\frac{|\sum_{t,r}\mathrm{RSS}^{(i)}_{t,r}b_{r}(\theta)b_{t}(\phi)|^{2}}{\sum_{t,r}||b_{r}(\theta)b_{t}(\phi)||^{2}}\right\\}$
(4)

$[\hat{\theta},\hat{\phi}]=\operatorname*{arg\,max}_{\theta,\phi}\left\\{\frac{|\sum_{t,r}\mathrm{RSS}^{(i)}_{t,r}b_{r}(\theta)b_{t}(\phi)|^{2}}{\sum_{t,r}||b_{r}(\theta)b_{t}(\phi)||^{2}}\right\\}\
.$ (5)

To solve this two-dimensional problem, we run a grid search over all
$\hat{\theta}$ and $\hat{\phi}$. For later discussion, with some abuse of
notation we refer to
$\mathbf{f(\theta,\phi)}=\frac{|\sum_{t,r}\mathrm{RSS}^{(i)}_{t,r}b_{r}(\theta)b_{t}(\phi)|^{2}}{\sum_{t,r}||b_{r}(\theta)b_{t}(\phi)||^{2}}$.

###  IV-B RSS-LS1D Algorithm

This algorithm is a variation of the algorithm mentioned in [8], where the
angle estimate is obtained by finding the strongest TX-RX beampattern pair
using BRP and minimizing over a fixed beam. Since we do not use BRP, we go
over each TX/RX beam pattern and minimize for $\theta$ or $\phi$ over the
summation of all such correlations, respectively.

$\hat{\theta}=\operatorname*{arg\,min}_{\theta}\min_{\alpha}\sum_{t,r}{}||{\mathrm{RSS}^{(i)}_{t,r}-\alpha
b_{r}(\theta)}||^{2}$ (6)

$\hat{\phi}=\operatorname*{arg\,min}_{\phi}\min_{\alpha}\sum_{t,r}||{\mathrm{RSS}^{(i)}_{t,r}-\alpha
b_{t}(\phi)}||^{2}$ (7)

Similarly, solving for $\alpha$ in (6) and (7) we obtain:

$\hat{\theta}=\operatorname*{arg\,max}_{\theta}\left\\{{|\sum_{t}\mathrm{R}^{(i)}_{t}B_{r}(\theta)|^{2}}\right\\}$
(8)

$\hat{\phi}=\operatorname*{arg\,max}_{\phi}\left\\{{|\sum_{r}\mathrm{R}^{(i)}_{r}B_{t}(\phi)|^{2}}\right\\}$
(9)

where $\mathrm{R}^{(i)}_{t}$ and $\mathrm{R}^{(i)}_{r}$ are the normalized
$\mathrm{RSS}^{(i)}_{t,r}$ values, such that
$\mathrm{R}^{(i)}_{t}=\frac{\mathrm{RSS}^{(i)}_{t,r}}{||\mathrm{RSS}^{(i)}_{t}||^{2}}$
and
$\mathrm{R}^{(i)}_{r}=\frac{\mathrm{RSS}^{(i)}_{t,r}}{||\mathrm{RSS}^{(i)}_{r}||^{2}}$.
Analogously, $B_{r}(\theta)$ = $\frac{b_{r}(\theta)}{||b_{r}(\theta)||^{2}}$
and $B_{t}(\phi)$ = $\frac{b_{t}(\phi)}{||b_{t}(\phi)||^{2}}$ denote the
normalized beampatterns. Since this is a one-dimensional projection of the
problem, the algorithm often locks to some local maxima instead of a global
one as we will see in the next section. For later discussion, with some abuse
of notation, we refer to
$\mathbf{\zeta(\theta)}={|\sum_{t}\mathrm{R}^{(i)}_{t}B_{r}(\theta)|^{2}}$ and
$\mathbf{\kappa(\phi)}={|\sum_{r}\mathrm{R}^{(i)}_{r}B_{t}(\phi)|^{2}}$.

##  V Results and Discussion

The proposed algorithm is validated using measured and simulated data. The
room at IMDEA is as shown in Fig. 1. We created a 2D ray-tracing-based model
assuming walls on all four sides. The windows, doors, and furniture were not
modeled. First and second-order MPCs were obtained and applied to the signal
model (1). The model is used only for verification of measured results and
therefore simplified. This generates I/Q sample data similar to what is
obtained by the measurement setup. Each signal component (LoS and NLoS)
experiences free space path loss and an additional 3 dB of attenuation per
reflection for NLoS. The transmitted pulse is the same Zadoff-Chu sequence
with $N_{\mathrm{K}}\times N_{\mathrm{\mathrm{BP},{TX}}}$ samples. For each
component, the AWGN is generated having power spectral density
$\mathrm{N_{0}}$, with $\mathrm{SNR}=10\log_{10}\frac{\lvert\lvert\alpha
b_{r}(\theta_{l})b_{t}(\phi_{l})\rvert\rvert^{2}}{\mathrm{N_{0}}}$, where
$\alpha$ is the magnitude of the free space path loss, and $b_{r}(\theta_{l})$
and $b_{t}(\phi_{l})$ are the receive and transmit beampatterns, respectively.

The synthetic I/Q samples use the beampattern matrix obtained from
measurements in an anechoic chamber (generated using the same codebook and
antenna array as for measurements) and are of dimension
$N_{\mathrm{\mathrm{BP},{TX}}}\times N_{\mathrm{\phi}}$ and
$N_{\mathrm{\mathrm{BP},{RX}}}\times N_{\mathrm{\theta}}$. Here
$N_{\mathrm{\mathrm{BP},{TX}}}=N_{\mathrm{\mathrm{BP},{RX}}}=64$ and
$N_{\mathrm{\theta}}=N_{\mathrm{\phi}}$ are the angular samples ranging from
$-79.5^{\circ}$ to $80^{\circ}$ with steps of $0.5^{\circ}$. The same measured
beampattern data is used for angle estimation.

(a) AoD (b) AoA

Figure 8: Plots depicting the objective functions for $\mathbf{RX}_{2}$ when TX is at position $\mathbf{p}_{5}$. The two plots compare the performance of RSS-LS2D algorithm with the RSS-LS1D algorithm. TABLE I: Estimated AoA and AoD values for TX at $\mathbf{p}_{1}$ | AoA | AoD
---|---|---
OR | RSS-LS2D | RSS-LS1D | OR | RSS-LS2D | RSS-LS1D
$\mathrm{RX}_{1}$ | $4^{\circ}$ | $4^{\circ}$ | $22.5^{\circ}$ | $-2.5^{\circ}$ | $-1^{\circ}$ | $-3^{\circ}$
$\mathrm{RX}_{2}$ | $-41.5^{\circ}$ | $-45^{\circ}$ | $-29^{\circ}$ | $9.5^{\circ}$ | $10^{\circ}$ | $-19.5^{\circ}$
$\mathrm{RX}_{3}$ | $22.5^{\circ}$ | $21.5^{\circ}$ | $37^{\circ}$ | $2^{\circ}$ | $2.5^{\circ}$ | $7^{\circ}$
$\mathrm{RX}_{4}$ | $5^{\circ}$ | $5.5^{\circ}$ | $14.5^{\circ}$ | $-7^{\circ}$ | $-6^{\circ}$ | $-3^{\circ}$
TABLE II: Estimated AoA and AoD values for TX at $\mathbf{p}_{5}$ | AoA | AoD
---|---|---
OR | RSS-LS2D | RSS-LS1D | OR | RSS-LS2D | RSS-LS1D
$\mathrm{RX}_{1}$ | $-15^{\circ}$ | $-16.5^{\circ}$ | $37^{\circ}$ | $-23.5^{\circ}$ | $-22^{\circ}$ | $-3^{\circ}$
$\mathrm{RX}_{2}$ | $32^{\circ}$ | $31^{\circ}$ | $37^{\circ}$ | $-15^{\circ}$ | $-15.5^{\circ}$ | $-1.5^{\circ}$
$\mathrm{RX}_{3}$ | $15.5^{\circ}$ | $13.5^{\circ}$ | $35^{\circ}$ | $-8.5^{\circ}$ | $-8^{\circ}$ | $27.5^{\circ}$
$\mathrm{RX}_{4}$ | $45^{\circ}$ | $45^{\circ}$ | $-3^{\circ}$ | $-41.5^{\circ}$ | $-40.5^{\circ}$ | $-1.5^{\circ}$

First, we compute the RSS matrix at each receiver for all TX positions. As a
sample, we can see in Fig. 6 and Fig. 7 the RSS plots for positions
$\mathbf{p}_{1}$ and $\mathbf{p}_{5}$ for all four receivers in comparison
with their simulated counterparts. Looking at each position in the room layout
and the geometry of the room we can predict some of the strongest paths for
each RX. Some positions of the TX perform better than others depending on
which receiver we observe. $\mathrm{RX}_{4}$ is the one facing the windows and
we observe that the data obtained is unreliable, our simulations perform
better because we do not model windows in our synthetic setup. If we were to
estimate the angle by observing the strongest RSS, the angle estimate obtained
is not bad due to the resolution of the beampatterns themselves. Still, it is
prone to some errors especially for receivers not facing the TX, as can be
seen for $\mathrm{RX}_{4}$ in a lot of cases. The white lines in the measured
RSS plots signify the missing samples of certain beampatterns that weren’t
strong enough to activate the RX.

The obtained angle estimates from the two algorithms for two sample positions
along with a rough estimate based on the observed RSS(OR) value are as in
Table I and II. We see that our proposed RSS-LS2D algorithm provides a more
accurate angle estimate based on the $\theta,\phi$ grid minimization compared
to RSS-LS1D which tends to vary from the true value. This is because the RSS-
LS1D algorithm locks on local maxima, if we were to fix the strongest TX/RX
beampattern it would be able to obtain the correct estimate. We notice, that
the angle obtained from RSS-LS1D is close to AoA/AoD of one of the MPCs just
not the strongest one.

To elaborate on the performance of RSS-LS1D we observe the objective functions
$\mathbf{\zeta(\theta)}$ and $\mathbf{\kappa(\phi)}$ as from (6) and (7) for
$\mathrm{RX}_{2}$ in Fig 8. We compare the function to a slice taken from
$\mathbf{f(\theta,\phi)}$ of the RSS-LS2D objective function from (5), where
$\theta$ and $\phi$ are maximized. We are essentially observing the objective
functions over which the $\operatorname*{arg\,max}$ is computed. For instance,
looking at Fig. 8(a) we see that the peak maximized by the RSS-LS1D algorithm
is present as a local maximum in the RSS-LS2D algorithm objective function
$\mathbf{f}(\theta,\phi)$. But where RSS-LS2D maximizes at $\phi=-15^{\circ}$,
RSS-LS1D maximizes at $\phi=-1.5^{\circ}$ (also seen in Table II). For Fig.
8(b), on the other hand, we see a similar behavior of the objective functions
and maximized values relatively close to each other. This is to emphasize that
we can observe plenty of cases where the RSS-LS1D tends to pick an MPC instead
of the strongest path since it doesn’t perform a maximization over the
complete RSS grid.

To further compare the accuracy of our proposed algorithm, we compute the CDF
of the error in AoA/AoD from the true value over all 10 receiver positions and
8 different snapshots of the measured data. The results are shown in Fig 9,
where we can see that our proposed RSS-LS2D algorithm performs well, with less
than a $10^{\circ}$ error in the majority of cases.

##  VI Conclusions

In this paper, we show that even with mmWave directional signaling and only
RSS information it is possible to exploit multipath information and obtain a
good NLoS angle estimate. Our proposed algorithm RSS-LS2D can successfully
estimate AoD and AoA for both LoS and NLoS cases with some stipulations. In
particular, we outperform the RSS-LS1D which requires information about the
strongest beampattern on at least one side to successfully estimate the angle
on the opposite side. Although we only estimated one dominant path per RX, in
the future we plan to estimate multiple AoA per RX array and add delay
estimates to improve position accuracy.

(a) Receiver 1 (b) Receiver 2 (c) Receiver 3 (d) Receiver 4

Figure 9: CDF plots for error in AoA/AoD estimation using RSS-LS2D and RSS-
LS1D over all positions.

## References

  * [1] F. Zafari, A. Gkelias, and K. K. Leung, “A survey of indoor localization systems and technologies,” _IEEE Communications Surveys & Tutorials_, vol. 21, no. 3, pp. 2568–2599, 2019.
  * [2] M. Shafi, J. Zhang, H. Tataria, A. F. Molisch, S. Sun, T. S. Rappaport, F. Tufvesson, S. Wu, and K. Kitao, “Microwave vs. millimeter-wave propagation channels: Key differences and impact on 5g cellular systems,” _IEEE Communications Magazine_ , vol. 56, no. 12, pp. 14–20, 2018.
  * [3] A. Shastri, N. Valecha, E. Bashirov, H. Tataria, M. Lentmaier, F. Tufvesson, M. Rossi, and P. Casari, “A review of millimeter wave device-based localization and device-free sensing technologies and applications,” _IEEE Communications Surveys & Tutorials_, 2022.
  * [4] P. Meissner, “Multipath-assisted indoor positioning,” PhD dissertation, Graz University of Technology, 2014. [Online]. Available: https://www2.spsc.tugraz.at/www-archive/downloads/Meissner
PhDThesis2014web.pdf

  * [5] P. Meissner and K. Witrisal, “Multipath-assisted single-anchor indoor localization in an office environment,” in _2012 19th International Conference on Systems, Signals and Image Processing (IWSSIP)_. IEEE, 2012, pp. 22–25.
  * [6] E. Leitinger, P. Meissner, and K. Witrisal, “Cognitive indoor positioning and tracking using multipath channel information,” _arXiv preprint arXiv:1610.05882_ , 2016.
  * [7] X. Li, E. Leitinger, and F. Tufvesson, “Rss-based localization of low-power iot devices exploiting aoa and range information,” in _2020 54th Asilomar Conference on Signals, Systems, and Computers_. IEEE, 2020, pp. 651–656.
  * [8] J. O. Lacruz, D. Garcia, P. J. Mateo, J. Palacios, and J. Widmer, “mm-flex: an open platform for millimeter-wave mobile full-bandwidth experimentation,” in _Proceedings of the 18th International Conference on Mobile Systems, Applications, and Services_ , 2020, pp. 1–13.
  * [9] Vadatech 2018 amc599 board. [Online]. Available: https://www.vadatech.com/media/AMC599_AMC599_Datasheet.pdf
  * [10] Vadatech 2018 amc726 board. [Online]. Available: https://www.vadatech.com/media/AMC726_AMC726_Datasheet.pdf
  * [11] Siversima 2018 evk06002 development kit. [Online]. Available: https://www.siversima.com/product/evk-600200


# Efficient Graph Reconstruction and Representation
Using Augmented Persistence Diagrams

Brittany Terese Fasy School of Computing and Department of Mathematical
Sciences, Montana State University<EMAIL_ADDRESS>Samuel Micka
Mathematics & Computer Science Department, Western Colorado University,
<EMAIL_ADDRESS>David L. Millman School of Computing, Montana State
University, Blocky, Inc<EMAIL_ADDRESS>Anna Schenfisch Department of
Mathematical Sciences, Montana State University<EMAIL_ADDRESS>Lucy Williams School of Computing, Montana State University,
<EMAIL_ADDRESS>
###### Abstract

Persistent homology is a tool that can be employed to summarize the _shape_ of
data by quantifying homological features. When the data is an object in
$\mathbb{R}^{d}$, the (augmented) persistent homology transform ((A)PHT) is a
family of persistence diagrams, parameterized by directions in the ambient
space. A recent advance in understanding the PHT used the framework of
reconstruction in order to find finite a set of directions to faithfully
represent the shape, a result that is of both theoretical and practical
interest. In this paper, we improve upon this result and present an improved
algorithm for graph—and, more generally one-skeleton—reconstruction. The
improvement comes in reconstructing the edges, where we use a radial binary
(multi-)search. The binary search employed takes advantage of the fact that
the edges can be ordered radially with respect to a reference plane, a feature
unique to graphs.

##  1 Introduction

At the heart of inverse problems in the field of topological data analysis is
the following question: _how many persistence diagrams are needed to
faithfully represent a shape?_ Since the introduction of persistence diagrams,
it has been known that many “shapes” can share the same persistence diagram.
With enough persistence diagrams, we arrive at a set of parameterized diagrams
(that is, diagrams labeled by direction) that uniquely, or _faithfully_
represents the underlying shape. Moreover, the set of parameterized diagrams
is _faithful_ if and only if it can be used to _reconstruct_ the underlying
shape.

The foundation for asking the question of how many diagrams are needed for a
faithful representation was first introduced in [16], where Turner et al.
defined the Persistent Homology Transform (PHT) and the Euler characteristic
curve transform (ECCT). These transforms map a simplicial complex embedded in
$\mathbb{R}^{d}$ (a shape) to sets of persistence diagrams (respectively,
Euler characteristic curves) parmeterized by $\mathbb{S}^{d-1}$, the set of
all directions in $\mathbb{R}^{d}$. They showed that, up to mild general
position assumptions, no two simplicial complexes can correspond to the same
PHT or ECCT, i.e., they showed that the uncountably infinite sets making up
the PHT and ECCT faithfully represent the shape.

However, the PHT and ECCT are uncountably infinite sets of diagrams (and Euler
characteristic curves). Thus, to bridge the gap between the theory and what
can be used in practice, a discretization of the PHT was needed; several
papers stepped up to the challenge and proved that there exists a finite
discretizations of topological transforms that are faithful for simplicial and
cubical complexes [1, 2, 4, 8, 5, 12, 7]. Beyond proving the existence of
these finite faithful sets, Belton et al. [2] explicitly give an algorithm for
using an oracle to reconstruct graphs embedded in $\mathbb{R}^{d}$ with $n$
vertices. Their reconstruction uses $n^{2}-n+d+1$ augmented persistence
diagrams in $O(dn^{d+1}+n^{4}+(d+n^{2})\Pi)$ time [2, Theorem 17], where
$\Theta(\Pi)$ is the time complexity it takes for an oracle to produce answer
a persistence diagram querry. Since the direction-labeled set of diagrams can
be used to reconstruct the underlying graph, the set is a faithful
discretization of the augmented PHT (APHT).

In the current work, we give a faithful discretization of the APHT using
$\Theta(\mathbb{V}_{\text{dgm}}+d+m\log n)$ diagrams and
$\Theta(\mathbb{V}_{\text{time}}+n^{2}+m\log n(\log n+d+\Pi))$ time. This
result is an improvement over [2], both in the size of the set and in the
speed of reconstruction. The crux is in the improvmeent to edge
reconstruction. While the method of [2] uses a linear scan of all possible
edges for each vertex, resulting in a quadratic number of diagrams needed,
here we show that we can detect edges with $\Theta(m\log n)$ diagrams using a
radial binary multi-search.

##  2 Background and Tools

In this section, we provide definitions of the tools used in the remainder of
the paper. We make use of standard notation such as using $e_{i}$ for the
$i$th standard basis vector in $\mathbb{R}^{d}$, where $1\leq i\leq d$. We use
the notation $(V,E)$ for a graph and its vertex and edge sets, and use $n=|V|$
and $m=|E|$. We assume the reader is familar with standard concepts in
topology (such as simplicial complexes, simplicial homology, Betti numbers),
and note that further information can be found in [6] for a general
introduction and in [7, Section 2.1] for a detailed definition of the
augmented persistence diagram.

First, we define our general postition assumption.

######  Assumption 1.

Let $V\subset\mathbb{R}^{d}$ be a finite set with $d\geq 2$. We say $V$ is in
_general position_ if the following properties are satisfied:

  1. (i)

Every set of $d+1$ points is affinely independent.

  2. (ii)

No three points are colinear after orthogonal projection into the space
$\pi(\mathbb{R}^{d})$, where $\pi:\mathbb{R}^{d}\to\mathbb{R}^{2}$ is the
orthogonal projection onto the plane spanned by the first two basis elements,
$e_{1}$ and $e_{2}$.

  3. (iii)

Every point has a unique height with respect to the direction $e_{2}$.

We call a graph _GP-immersed_ 111Here, we use _immersed_ rather than
_embedded_ in order to allow intersections of edges. Note, however, that this
can only happen when $d=2$. iff its vertex set is in general position in
$\mathbb{R}^{d}$.

We note that iii is not strictly necessary, however, it is convenient for
simplicity of exposition. How to handle this degeneracy is discussed in
Appendix B.

Given a graph GP-immersed in $\mathbb{R}^{d}$, we can filter the graph based
on the height in any direction $s$ in the sphere of directions
$\mathbb{S}^{d-1}$. To do so, we assign each simplex a height. A vertex $v\in
V$ is assigned the height $s\cdot v$, and an edge $[v_{0},v_{1}]\in E$ is
assigned the height $\max\\{s\cdot v_{0},s\cdot v_{1}\\}$. This function,
mapping vertices and edges to heights, is known as a _filter_ function, which
we use to compute persistent homology.

###  2.1 Persistence and the Oracle Framework

Given a filtered simplicial complex (that is, a simplicial complex with each
simplex assigned a “height”), the corresponding _augmented persistence diagram
(APD)_ is a record of of all homological events throughout the filtration. A
_birth event_ is the introduction or appearance of a new homological feature,
and a _death event_ is the merging of two features. A death is paired with the
most recent of the birth-labeled features that it merges together, creating a
birth-death pairing, leaving the remaining feature labeled by the elder birth
height. In an APD, every simplex corresponds to exactly one event (resulting
in some pairings where the birth and death heights are equal). As a result, by
construction, APDs contain at least one event at the height of each vertex.

For a simplicial complex (e.g., a graph) GP-immersed in $\mathbb{R}^{d}$ and a
direction in $S^{d-1}$, we use the _lower-star filtration_ : the nested
sequence of graphs that arise by looking at all simplices at or below a given
height and allowing that height to grow from $-\infty$ to $\infty$. Throughout
this paper, we denote the $i$-dimensional APD by
$\widehat{\mathcal{D}}_{i}(s)$ and write
$\widehat{\mathcal{D}}(s)=\sqcup_{i}\widehat{\mathcal{D}}_{i}(s)$, omitting
the graph from the notation (as it is always clear from context).222 When
calculating diagrams, we count $\widehat{\mathcal{D}}(s)$ as one diagram, not
multiple.

The _(augmented) persistence homology transform ((A)PHT)_ is the set of
(augmented) diagrams of lower-star filtrations in all possible directions,
parameterized by direction. That is, the set
$X=\\{(\widehat{\mathcal{D}}(s),s)\\}_{s\in\mathbb{S}^{d-1}}$. A _faithful
discretization_ is a finite subset of $X$ from which all other elements of $X$
can be deduced (and, by [16], corresponds to a unique simplicial complex). The
introduction of the (augmented) persistent homology transform has sparked
related research in applications of shape comparison [9, 11, 14, 10, 17, 18,
15, 3]. As such, finding a minimal faithful discretization is important for
the applicability of the (A)PHT. In what follows, we will only consider APDs,
and we may shorten notation and refer to an APD by the word _diagram_.

In this work, we assume an oracle framework. That is, we assume that we have
no knowledge of the shape itself, but we have access to an oracle from which
we can query directional diagrams.

######  Definition 1 (Oracle).

For a graph $(V,E)$ GP-immersed in $\mathbb{R}^{d}$ and a direction
$s\in\mathbb{S}^{d-1}$, the operation ${{\texttt{Oracle}}(s)}$ returns the
diagram $\widehat{\mathcal{D}}(s)$. We define $\Theta(\Pi)$ to be the time
complexity of this oracle query and note that the space complexity of
$\widehat{\mathcal{D}}(s)$ is $\Theta(n+m)$. We assume that the data structure
returned by the oracle allows queries for specific birth or death values in
$\Theta(\log n)$ time (for example, the we could have two arrays of
persistences points, one sorted by birth values and one sorted by death
values).

###  2.2 Constructions and Data Structures

In this subsection, we introduce the edge arc object and other definitions
useful for computing properties of immersed graphs. Throughout this paper, we
project points in $\mathbb{R}^{d}$ to the $(e_{1},e_{2})$-plane. As a result,
we use “above (below)” without stating with respect to which direction as
shorthand for “above (below) with respect to the direction $e_{2}$.” This
direction is intentionally chosen (and used in our GP assumption), as it
corresponds to our intuition of above (below) in the figures. When we measure
an angle of a vector $x$, denoted $\measuredangle x$, we mean the angle that
$\pi(x)$ makes with the positive $e_{1}$ axis.

Given a direction $s$ and a vertex in a graph immersed in $\mathbb{R}^{d}$, we
classify each edge $(v,v^{\prime})$ as either an “incoming” edge, when
$v^{\prime}$ is below $v$ with respect to $s$, or an “outgoing” edge, when
$v^{\prime}$ is above $v$ with respect to $s$. Note that all incoming edges
have the same height as the vertex with respect to the $e_{2}$ direction.

######  Definition 2 (Indegree).

Let $(V,E)$ be a graph GP-immersed in $\mathbb{R}^{d}$. Let $v\in V$ and
$s\in\mathbb{S}^{d-1}$. The _indegree of $v$ in direction $s$_, denoted
${\texttt{Indeg}}(v,s)$, is the number of edges incident to $v$ with height
$s\cdot v$.

The following lemma relates the number of edges at a given height to points in
the APD.

######  Lemma 1 (Edge Count).

Let $(V,E)$ be a graph and let $c\in\mathbb{R}$. Let $f\colon V\sqcup
E\to\mathbb{R}$ be a filter function. Then, the edges in $E$ with a function
value of $c$ are in one-to-one correspondence with the following multiset of
points in $\widehat{\mathcal{D}}(f)$, the diagram corresponding to $f$:

$\displaystyle\begin{split}\big{\\{}(b,d)&\in\widehat{\mathcal{D}}_{1}(f)\text{
s.t. }b=c\big{\\}}\\\ &\cup\big{\\{}(b,d)\in\widehat{\mathcal{D}}_{0}(f)\text{
s.t. }d=c\big{\\}}.\end{split}$ (1)

In other words, each edge corresponds to either a birth of a one-dimensional
homological feature or a death of a zero-dimensional feature in
$\widehat{\mathcal{D}}(f)$. For more details and a generalized proof, see [7,
Appendix A].

If $f$ is a lower-star filtration in direction $s\in\mathbb{S}^{d-1}$, we note
that whenever a vertex $v$ has a unique height with respect to a direction
$s$, the cardinality of the multiset above is exactly ${\texttt{Indeg}}(v,s)$.

######  Lemma 2 (Indegree Computation).

Let $(V,E)$ be a graph GP-immersed in $\mathbb{R}^{d}$. Let $v\in V$ and let
$s\in\mathbb{S}^{d-1}$ such that $s\cdot v\neq s\cdot v^{\prime}$ for any
$v^{\prime}\neq v\in V$. Then, ${\texttt{Indeg}}(v,s)$ can be computed via the
oracle using one diagram and $\Theta(\log n+\Pi)$ time.

###### Proof.

Let $\widehat{\mathcal{D}}={{\texttt{Oracle}}(s)}$. By the assumption on $s$,
the height of $v$ with respect to the direction $s$ is unique. Hence, we know
that any edge at height $c=s\cdot v$ must be incident to $v$. Thus, by the
definition of indegree, an edge has the height $c$ if and only if it
contributes to the indegree of $v$ in direction $s$. Using Lemma 1, we count
these edges by counting one-dimensional births and zero-dimensional deaths at
height $c$. Since $\widehat{\mathcal{D}}_{0}$ and $\widehat{\mathcal{D}}_{1}$
are sorted by both birth and death values (see Definition 1) and since
$\widehat{\mathcal{D}}$ has $\Theta(n+m)$ points, searching for these events
takes $\Theta(\log n+\log m)$. Adding $\Theta(\Pi)$ for the oracle query and
recalling that $m=O(n)$, the total runtime is $\Theta(\log n+\Pi)$. ∎

Table 1: Attributes of the edge arc object. lightgray ${\mathit{EA}}$ |  Edge Arc
---|---
${\mathit{v}}$ |  Vertex around which the edge arc is centered
$(\alpha_{1},\alpha_{2})$ |  Start and stop angles of the arc, with respect to the $e_{1}$ direction
${\mathit{verts}}$ |  Array of vertices in arc radially ordered clockwise in $(e_{1},e_{2})$-plane
${\mathit{count}}$ |  Number of edges incident to $v$ within the arc

We conclude this section by introducing a data structure, the _edge arc
object_ ; see Table 1 for a summary of the attributes of an edge arc and
Figure 1 for an example. An edge arc represents the region in the
$(e_{1},e_{2})$-plane centered at $v$ that is swept out between the two angles
$\alpha_{1}$ and $\alpha_{2}$ (the word ‘arc’ is referring to the arc of
angles between $\alpha_{1}$ and $\alpha_{2}$, where the angle is measured with
respect to the postive $e_{1}$ axis). We only consider edge arcs in the upper
half-space, with respect to the $e_{2}$ direction, so the maximal edge arc is
the upper half-plane and the start and stop angles always satisfy
$0\leq\alpha_{1}\leq\alpha_{2}\leq\pi$. An edge arc stores an array of
vertices sorted radially clockwise about $\pi(v)$ in the $(e_{1},e_{2})$-plane
in decreasing angle with the $e_{1}$-direction. By construction, the first
vertex in the array must be closest to $\alpha_{2}$ and the last closest to
$\alpha_{1}$. The edge arc also stores the count of edges of $E$ that have
vertices from ${\mathit{verts}}$ as endpoints. In implementation, angles
$\alpha_{1}$ and $\alpha_{2}$ do not need to be stored directly, but we
include them in psuedocode and discussions for clarity.

Figure 1: An edge arc ${\mathit{EA}}$ centered at vertex
${\mathit{EA}}.{\mathit{v}}=v$. Other attributes of the edge arc include its
start and stop angles, ${\mathit{EA}}.\alpha_{1}=1.75$ radians and
${\mathit{EA}}.\alpha_{2}=\pi\approx 3.14$ radians, the array of vertices
${\mathit{EA}}.{\mathit{verts}}=\\{v_{1},v_{2}\\}$, and the count of edges
${\mathit{EA}}.{\mathit{count}}=1$. Here, we also see that
${\texttt{Indeg}}(v,e_{1})=1$ and ${\texttt{Indeg}}(v,e_{2})=2$.

Given some arc ${\mathit{EA}}$ centered at vertex $v\in V$, we need to be able
to compute ${\mathit{EA}}.{\mathit{count}}$, the number of edges contained
${\mathit{EA}}$ that are adjacent to $v$. The following lemma provides such a
computation. We omit a proof because it is a straightforward adaptation of [7,
Theorem $16$] and [2, Lemma $13$].

######  Lemma 3 (Arc Count).

Let $(V,E)$ be a graph GP-immersed in $\mathbb{R}^{d}$. Let ${\mathit{EA}}$ be
an edge arc object, and let $v={\mathit{EA}}.{\mathit{v}}$. Let
$s\in\mathbb{S}^{d-1}$ be the direction perpendicular to $\alpha_{2}$ so that
the arc is entirely below $s\cdot v$. Let $E_{*}$ denote the edges with height
$s\cdot v$ that are not contained in ${\mathit{EA}}$. If no vertex in $V$ is
at the same height as $v$ in direction $s$, then

$\displaystyle{\mathit{EA}}.{\mathit{count}}={\texttt{Indeg}}(v,s)-|E_{*}|.$

As an illustration, again, consider Figure 1. Consider the direction
$s=e^{i(\alpha_{2}-\pi/2)}$, which is perpendicular to $e^{i\alpha}$ by
construction. In addition, the edge arc is below $s\cdot v$ (specifically, all
vertices in ${\mathit{EA}}.{\mathit{verts}}$ are below $s\cdot v$). Then,
${\texttt{Indeg}}(v,s)=2$ and $E_{*}=\\{[v,v_{5}]\\}$. By Lemma 3,
${\mathit{EA}}.{\mathit{count}}={\texttt{Indeg}}(v,s)-|E_{*}|=2-1=1$. When we
say that a list of vertices or edges is sorted clockwise around $v$, we mean
that the list is sorted clockwise (cw) around $\pi(v)$ once projected into the
$(e_{1},e_{2})$-plane with the largest angle first.

##  3 Fast Reconstruction

In this section, we provide an algorithm to reconstruct a graph (and, more
generally, a one-skeleton of a simplicial complex) using the oracle. We start
with an algorithm to find the edges, provided the vertex locations are known.
We end with describing the complete graph reconstruction method.

###  3.1 Fast Edge Reconstruction

In this subsection, we assume we have a graph $(V,E)$, where the vertex set
$V$ is known, but $E$ is unknown. Using the oracle and the known vertex
locations, we provide a reconstruction algorithm to find all edges in $E$
(Algorithm 3). This algorithm is a sweepline algorithm in direction $e_{2}$
that, for each vertex processed in the sweep, performs a radial binary multi-
search (Algorithm 2). This search is enabled by an algorithm that splits an
edge arc object into two edge arcs, each containing half of the vertices
(Algorithm 1). We provide the algorithms and relevant theorem statements here,
but defer the proofs to Appendix A.

Algorithm 1 ${\texttt{SplitArc}}({\mathit{EA}},\mathit{bigedges},\theta)$

0: ${\mathit{EA}}$, an edge arc; $\mathit{bigedges}$, an array of all edges
$({\mathit{EA}}.{\mathit{v}},v^{\prime})\in E$ such that
$\measuredangle\pi(v^{\prime}-{\mathit{EA}}.{\mathit{v}})<{\mathit{EA}}.\alpha_{1}$;
$\theta$, the minimum angle defined by any three vertices in $\pi(V)$

0: ${\mathit{EA}}_{\ell}$ and ${\mathit{EA}}_{r}$, edge arcs satisfying the
properties in Theorem 1

1: $n_{v}\leftarrow|{\mathit{EA}}.{\mathit{verts}}|$

2: $mid\leftarrow\lceil\frac{n_{v}}{2}\rceil$

3:
$\alpha\leftarrow\measuredangle\pi({\mathit{EA}}.{\mathit{verts}}[mid]-{\mathit{EA}}.{\mathit{v}})-\theta/2$

4: $s\leftarrow e^{i(\alpha-\frac{\pi}{2})}$

5: $m_{\ell}\leftarrow{\texttt{Indeg}}({\mathit{EA}}.{\mathit{v}},s)-|\\{b\in
bigedges~{}|~{}\measuredangle b<\pi+\alpha\\}|$

6: $m_{r}\leftarrow{\mathit{EA}}.{\mathit{count}}-m_{\ell}$

7: ${\mathit{EA}}_{\ell}\leftarrow$ edge arc where
${\mathit{EA}}_{\ell}.{\mathit{v}}={\mathit{EA}}.{\mathit{v}}$,
${\mathit{EA}}_{\ell}.\alpha_{1}=\alpha$,
${\mathit{EA}}_{\ell}.\alpha_{2}={\mathit{EA}}.\alpha_{2}$,
${\mathit{EA}}_{\ell}.{\mathit{verts}}={\mathit{EA}}.{\mathit{verts}}[:mid]$,
and ${\mathit{EA}}_{\ell}.{\mathit{count}}=m_{\ell}$

8: ${\mathit{EA}}_{r}\leftarrow$ edge arc where
${\mathit{EA}}_{r}.{\mathit{v}}={\mathit{EA}}.{\mathit{v}}$,
${\mathit{EA}}_{r}.\alpha_{1}={\mathit{EA}}.\alpha_{1}$,
${\mathit{EA}}_{r}.\alpha_{2}=\alpha$,
${\mathit{EA}}_{r}.{\mathit{verts}}={\mathit{EA}}.{\mathit{verts}}[mid+1:]$,
and ${\mathit{EA}}_{r}.{\mathit{count}}=m_{r}$

9: return $({\mathit{EA}}_{\ell}$, ${\mathit{EA}}_{r})$

In Algorithm 1, we find a direction $s$ in the $(e_{1},e_{2})$-plane such that
half of the vertices in ${\mathit{EA}}.{\mathit{verts}}$ are above $v$ and
half are below $v$ with respect to the direction $s$. This allows us to create
a new edge arcs corresponding to each half; see Figure 2. The properties of
Algorithm 1 are described in the following theorem.

Figure 2: The splitting of edge arc ${\mathit{EA}}$ into
${\mathit{EA}}_{\ell}$ and ${\mathit{EA}}_{r}$, as in Algorithm 1. The large
gray region is the region containing all edges of $bigedges$. That is, all
edges whose angle with the positive $e_{1}$-axis is at least
${\mathit{EA}}.\alpha_{1}$. On Line 5 of the algorithm, we compute the number
of edges in ${\mathit{EA}}_{\ell}$ by first computing the indegree of
${\mathit{EA}}.{\mathit{v}}$ in direction $s$ from the diagram in direction
$s$, then we subtract the number of edges in $bigedges$ that are below the
height of ${\mathit{EA}}.{\mathit{v}}$ in direction $s$ (i.e., below the blue
line). By the pigeonhole principal, we find
${\mathit{EA}}_{r}.{\mathit{count}}={\mathit{EA}}.{\mathit{count}}-{\mathit{EA}}_{\ell}.{\mathit{count}}$.

######  Theorem 1 (Arc Splitting).

Algorithm 1 uses one diagram and $\Theta(\log n+d+\Pi)$ time to split
${\mathit{EA}}$ into two new edge arcs ${\mathit{EA}}_{\ell}$ and
${\mathit{EA}}_{r}$ with the properties:

  1. (i)

The sets ${\mathit{EA}}_{r}.{\mathit{verts}}$ and
${\mathit{EA}}_{\ell}.{\mathit{verts}}$ partition the vertex set
${\mathit{EA}}.{\mathit{verts}}$ such that the vertices in
${\mathit{EA}}_{\ell}.{\mathit{verts}}$ come before those in
${\mathit{EA}}_{r}.{\mathit{verts}}$, with respect to the clockwise ordering
around ${\mathit{EA}}.{\mathit{v}}$.

  2. (ii)

$|{\mathit{EA}}_{\ell}.{\mathit{verts}}|=\lceil{\frac{1}{2}|{\mathit{EA}}.{\mathit{verts}}|}\rceil$.

  3. (iii)

$|{\mathit{EA}}_{r}.{\mathit{verts}}|=\lfloor{\frac{1}{2}|{\mathit{EA}}.{\mathit{verts}}|}\rfloor$.

Algorithm 2
${\texttt{UpEdges}}({v,V_{v},in_{v},\theta,\widehat{\mathcal{D}}})$

0: $v\in V$; $V_{v}$, array of all vertices in $V$ above $v$, ordered
clockwise; $in_{v}$, array of all incoming edges of $v$, sorted radially
clockwise; $\theta$, the minimum angle formed by any three vertices in
$\pi(V)$; and $\widehat{\mathcal{D}}$, the APD in direction $e_{2}$

0: array of all outgoing edges of $v$

1: $\mathit{indeg}\leftarrow$ indegree of $v$ in direction $-e_{2}$.

2: $\mathit{eastack}\leftarrow$ a stack of edge arc objects, initialized with
a single edge arc $A$, where $A.{\mathit{v}}=v$, $A.\alpha_{1}=0$,
$A.\alpha_{2}=\pi$, $A.{\mathit{count}}=\mathit{indeg}$, and
$A.{\mathit{verts}}=V_{v}$

3: $E_{v}\leftarrow\emptyset$

4: while $\mathit{eastack}$ is not empty do

5: ${\mathit{EA}}\leftarrow\mathit{eastack}.pop()$

6: if ${\mathit{EA}}.{\mathit{count}}=0$ then

7: Continue to top of while loop

8: else if ${\mathit{EA}}.{\mathit{count}}=|{\mathit{EA}}.{\mathit{verts}}|$
then

9: Append $v\times{\mathit{EA}}.{\mathit{verts}}$ to ${E}_{v}$, in order

10: Continue to top of while loop

11: else

12:
$({\mathit{EA}}_{\ell},{\mathit{EA}}_{r})\leftarrow{\texttt{SplitArc}}({\mathit{EA}},in_{v}\cup{E}_{v},\theta)$

13: Push ${\mathit{EA}}_{r}$ onto $\mathit{eastack}$

14: Push ${\mathit{EA}}_{\ell}$ onto $\mathit{eastack}$

15: end if

16: end while

17: return $E_{v}$

In Algorithm 2, we use Algorithm 1 to find all outgoing edges from a given
vertex. In particular, the algorithm maintains a stack of edge arc objects.
When processing an edge arc (the while loop in Lines 4–16), we are determining
which of the vertices in ${\mathit{verts}}$ form edges with ${\mathit{v}}$. If
an edge arc has ${\mathit{count}}=0$, it contains no edges, and it can be
ignored (Lines 6–8). If it has ${\mathit{count}}$ exactly equal to the number
of vertices in ${\mathit{verts}}$, each vertex in ${\mathit{verts}}$ must form
an edge with $v$ (Lines 8–11). Otherwise, as demonstrated in Figure 3, the
edge arc is split in half using Algorithm 1 and each half is put on the stack
to be processed.

(a) At vertex $v$.

(b) Split edge arc.

(c) Split again.

Figure 3: We demonstrate one step of Algorithm 2. (3(a)) By assumption, we
initially know $[v_{5},v]\in E$. From Line 1 of Algorithm 2 we also know that
two of the four vertices above $v$ are adjacent to $v$. Thus, we create an
edge arc object ${\mathit{EA}}$ with ${\mathit{EA}}.{\mathit{count}}=2$, and
${\mathit{EA}}.{\mathit{verts}}=(v_{1},v_{2},v_{3},v_{4})$. (3(b)) In
Algorithm 1, we choose a direction $s$ such that half of the vertices in
${\mathit{EA}}$ are below $v$. We use this split to create two edge arcs,
${\mathit{EA}}_{r}$ and ${\mathit{EA}}_{\ell}$, corresponding to the pink
shaded regions on the right and left of the blue line defined by $s$. We push
${\mathit{EA}}_{r}$ onto a stack to be processed later and focus on the arc
${\mathit{EA}}_{\ell}$. Since two edges contribute to $v$’s indegree in
direction $s$ and one is the known edge $[v_{5},v]$, we have
${\mathit{EA}}.{\mathit{count}}=2-1=1$. (3(c)) Next, we find a new direction
$s$ that splits ${\mathit{EA}}_{\ell}.{\mathit{verts}}$ into two sets of size
one. We push the set above $s$ onto our stack. The edge arc containing only
$v_{1}$ also has ${\mathit{EA}}.{\mathit{count}}=2-1=1$, so $[v_{1},v]\in E$.
After all steps of Algorithm 2 are applied to find the edges above a
particular vertex, Algorithm 2 is then applied to the next highest vertex,
eventually processing every vertex in $V$ in a sweep (Algorithm 3).

######  Theorem 2 (Finding Edges Above a Vertex).

Algorithm 2 finds the sorted array of edges above $v$ using
$\Theta(\deg(v)\log n)$ augmented persistence diagrams in $\Theta((\deg(v)\log
n)(\log n+d+\Pi))$ time.

Finally, our main algorithm (Algorithm 3) is a sweepline algorithm, where we
consider the vertices in increasing order of their $e_{2}$-coordinates and
find the outgoing edges of the vertex being considered.

Algorithm 3 ${\texttt{FindEdges}}(V)$

0: $V$, array of all vertices in the unknown graph

0: $E$, array of all edges in the unknown graph

1: $\widehat{\mathcal{D}}\leftarrow{{\texttt{Oracle}}(-e_{2})}$

2: $E\leftarrow\\{\\}$

3: $\mathit{vertsabove}\leftarrow$ for each $v\in V$, an array clockwise
ordering all vertices in $V$ that are above $v$

4: $\theta\leftarrow$ min angle defined by any three vertices of $\pi(V)$

5: for $v$ in $V$, in increasing height in direction $e_{2}$ do

6: $\mathit{inedges}\leftarrow$ clockwise sorted array of edges in $E$
incident to $v$

7:
$E+={\texttt{UpEdges}}({v,\mathit{vertsabove}[v],\mathit{inedges},\theta,\widehat{\mathcal{D}}})$

8: end for

9: return $E$

######  Theorem 3 (Edge Reconstruction).

Let $(V,E)$ be a graph GP-immersed in $\mathbb{R}^{d}$. Given $V$, Algorithm 3
reconstructs $E$ using $\Theta(m\log n)$ augmented persistence diagrams in
$\Theta(n^{2}+m\log n(\log n+d+\Pi))$ time.

###  3.2 Putting it Together: Full Reconstruction

The results of Section 3.1 are related to just part of the full process of
reconstruction, since reconstruction begins with no knowledge of the
underlying simplicial complex. Identifying the location of all vertices is the
first step, and is one that has been previously examined in detail. In
particular, Belton et al. provide an algorithm to reconstruct $V$ in
$\Theta(dn^{d+1}+d\Pi)$ time and $d+1$ oracle queries, assuming stricter
general position assumptions; see [2, Algorithm 1 & Theorem 9]. An appropriate
vertex reconstruction method taken together with Theorem 3, gives us the
following for a full reconstruction process.

######  Theorem 4 (Graph Reconstruction).

Suppose we use a vertex reconstruction method that requires
$\Theta(\mathbb{V}_{\text{time}})$ time and $\Theta(\mathbb{V}_{\text{dgm}})$
diagrams. Then we can reconstruct an unknown graph immersed in
$\mathbb{R}^{d}$ using $\Theta(\mathbb{V}_{\text{dgm}}+d+m\log n)$ diagrams in
$\Theta(\mathbb{V}_{\text{time}}+n^{2}+m\log n(\log n+d+\Pi))$ time.

We omit a proof of Theorem 4, as it simply combines the results of vertex
reconstruction methods (such as [2, Theorem 9]) and Theorem 3 of the current
paper. Note that, taking the stricter general position assumptions of [2],
Theorem 4 gives us a runtime of $\Theta(dn^{d+1}+d\Pi+n^{2}+m\log n(\log
n+d+\Pi))$ with $\Theta(d+m\log n)$ oracle queries.

Observing that the methods presented here are immediately applicable in the
reconstruction of one-skeletons of general simplicial complexes, we have the
following corollary:

######  Corollary 1 (One-Skeleton Reconstruction).

Let $K$ be an unknown simplicial complex GP-immersed in $\mathbb{R}^{d}$. A
vertex reconstruction algorithm that requires
$\Theta(\mathbb{V}_{\text{time}})$ time and $\Theta(\mathbb{V}_{\text{dgm}})$
diagrams taken together with Algorithm 3 of the current paper reconstruct the
one-skeleton of $K$ using $\Theta(\mathbb{V}_{\text{dgm}}+d+m\log n)$
augmented persistence diagrams in $\Theta(\mathbb{V}_{\text{time}}+n^{2}+m\log
n(\log n+d+\Pi))$ time.

Finally, we note that embedding a graph (or simplicial complex) in
$\mathbb{R}^{2}$ is a special case, as $m=O(n)$ and $d$ is constant. This
speeds up edge reconstruction as well as potentially, vertex reconstruction.
For instance, with the stricter general position assumptions, the method of
[2, Theorem 6] gives a vertex reconstruction of a graph embedded in
$\mathbb{R}^{2}$ that uses three diagrams and $\Theta(n\log n+\Pi)$ time.
Broadening our view again to general vertex reconstruction algorithms, we
obtain a result for plane graph reconstruction:

######  Corollary 2 (Reconstruction in $\mathbb{R}^{2}$).

Suppose we use a vertex reconstruction method that requires
$\Theta(\mathbb{V}_{\text{time}})$ time and $\Theta(\mathbb{V}_{\text{dgm}})$
diagrams. We can use an oracle to reconstruct the one-skeleton of an unknown
simplicial complex embedded in $\mathbb{R}^{2}$ using
$O(\mathbb{V}_{\text{dgm}}+n\log n)$ diagrams and
$O(\mathbb{V}_{\text{time}}+n^{2}+n\log n(d+\Pi))$ time.

##  4 Discussion

One way of proving that a discretization of the APHT is faithful is through
the method of _reconstructing_ the underlying simplicial complex. That is, by
showing that the underlying simplicial complex can be recovered with the data
of the discretization alone. In this paper, we take that approach and provide
an algorithm for reconstructing a graph immersed in $\mathbb{R}^{d}$. We use
fewer persistence diagrams than presented in alternate approaches. For
example, the algorithm that we present for edge reconstruction (when the
vertex locations are known) uses $\Theta(m\log n)$ diagrams. In contrast, [2,
Theorem 16] uses $n^{2}-n$ diagrams. Note that, for a very dense edge set,
that is, when $m=\Theta(n^{2})$, the method in [2, Theorem 16] uses fewer
diagrams. However, if $m=O(n)$, as is common in many complexes, the
represntation computed in this paper has fewer diagrams. Moreover, we
emphasize that the number of diagrams is not exponential in the ambient
dimension.

One might hope to use binary search strategies to reconstruct a simplicial
complexe, but the methods presented here are unique to one-skeletons. Radially
ordering higher dimensional simplices is not well-defined, and this issue
prevents the methods presented here from being immediately transferrable. On
the other hand, with the representation in this paper being output-sensitive
(as opposed to testing if every pair of vertices is a simplex), we have hope
for the discretization of the (A)PHT of a simplicial complex immersed in
$\mathbb{R}^{d}$ being proportional to the size of the complex itself.

We also observe that not all diagrams used in our reconstruction algorithms
were strictly necessary (i.e., the set of diagrams used were not a minimal
faithful set). One straightforward way to reduce the number of diagrams used
without altering the method much would be to split the region above a vertex
in the sweep into arcs that contain exactly the same number of edges as
vertices, or no edges. This property can then be validated by a simple
difference of indegrees. In ongoing work, we hope to make these claims
precise. We also hope to extend our methods to use topological descriptors
that are _not_ dimension-returning (such as augmented Euler Characteristic
curves).

## References

  * [1] R. L. Belton, B. T. Fasy, R. Mertz, S. Micka, D. L. Millman, D. Salinas, A. Schenfisch, J. Schupbach, and L. Williams.  Learning simplicial complexes from persistence diagrams.  In Canadian Conference on Computational Geometry, August 2018.  Also available at arXiv:1805.10716.
  * [2] R. L. Belton, B. T. Fasy, R. Mertz, S. Micka, D. L. Millman, D. Salinas, A. Schenfisch, J. Schupbach, and L. Williams.  Reconstructing embedded graphs from persistence diagrams.  Computational Geometry: Theory and Applications, 2020.
  * [3] P. Bendich, J. S. Marron, E. Miller, A. Pieloch, and S. Skwerer.  Persistent homology analysis of brain artery trees.  The Annals of Applied Statistics, 10(1):198, 2016.
  * [4] L. M. Betthauser.  Topological Reconstruction of Grayscale Images.  PhD thesis, University of Florida, 2018.
  * [5] J. Curry, S. Mukherjee, and K. Turner.  How many directions determine a shape and other sufficiency results for two topological transforms.  arXiv:1805.09782, 2018.
  * [6] H. Edelsbrunner and J. Harer.  Computational Topology: An Introduction.  American Mathematical Society, 2010.
  * [7] B. T. Fasy, S. Micka, D. L. Millman, A. Schenfisch, and L. Williams.  A faithful discretization of the augmented persistent homology transform.  2022\.  arXiv:1912.12759v4.
  * [8] R. Ghrist, R. Levanger, and H. Mai.  Persistent homology and Euler integral transforms.  Journal of Applied and Computational Topology, 2(1-2):55–60, 2018\.
  * [9] C. Giusti, E. Pastalkova, C. Curto, and V. Itskov.  Clique topology reveals intrinsic geometric structure in neural correlations.  Proceedings of the National Academy of Sciences, 112(44):13455–13460, 2015.
  * [10] P. Lawson, A. B. Sholl, J. Q. Brown, B. T. Fasy, and C. Wenk.  Persistent homology for the quantitative evaluation of architectural features in prostate cancer histology.  Scientific Reports, 9, 2019.
  * [11] Y. Lee, S. D. Barthel, P. Dłotko, S. M. Moosavi, K. Hess, and B. Smit.  Quantifying similarity of pore-geometry in nanoporous materials.  Nature Communications, 8:15396, 2017.
  * [12] S. A. Micka.  Searching and Reconstruction: Algorithms with Topological Descriptors.  PhD thesis, Montana State University, 2020.
  * [13] D. L. Millman and V. Verma.  A slow algorithm for computing the Gabriel graph with double precision.  Proceedings of the 23rd Annual Canadian Conference on Computational Geometry, 2011.
  * [14] A. H. Rizvi, P. G. Camara, E. K. Kandror, T. J. Roberts, I. Schieren, T. Maniatis, and R. Rabadan.  Single-cell topological RNA-seq analysis reveals insights into cellular differentiation and development.  Nature Biotechnology, 35(6):551, 2017.
  * [15] G. Singh, F. Mémoli, and G. E. Carlsson.  Topological methods for the analysis of high dimensional data sets and 3d object recognition.  SPBG, 91:100, 2007.
  * [16] K. Turner, S. Mukherjee, and D. M. Boyer.  Persistent homology transform for modeling shapes and surfaces.  Information and Inference: A Journal of the IMA, 3(4):310–344, 2014\.
  * [17] S. Tymochko, E. Munch, J. Dunion, K. Corbosiero, and R. Torn.  Using persistent homology to quantify a diurnal cycle in hurricanes.  Pattern Recognition Letters, 2020.
  * [18] Y. Wang, H. Ombao, and M. K. Chung.  Statistical persistent homology of brain signals.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1125–1129. IEEE, 2019.

##  Appendix A Algorithmic Proofs

In this appendix, we provide the proofs omitted from Section 3. These proofs
provide justification for the runtimes, diagram compelxity, and correctness of
the algorithms presented in this paper.

###  A.1 Proof of Theorem 1

See 1

###### Proof.

For the runtime, we walk through the algorithm and analyze the time and
diagram complexity of each line. In Lines 1–3, we find the angle $\alpha$ that
splits ${\mathit{EA}}.{\mathit{verts}}$ into two equal sets, then in Line 4
compute a direction $s$ orthogonal to $\alpha$. See Figure 2. Lines 1–4 use no
diagrams and can be done in constant time when restricting our attention to
the $(e_{1},e_{2})$-plane. However, we need $s$ to be a direction in
$\mathbb{R}^{d}$ (as opposed to only in the $(e_{1},e_{2})$-plane), so the
computation takes $\Theta(d)$ time.333With some clever data structures, this
$\Theta(d)$ can be reduced to constant time. For example, we could require
vectors in $\mathbb{R}^{2}$ are automatically padded with $0$’s to become
vectors in $\mathbb{R}^{d}$ when needed. However, this is out of the scope of
the real RAM model of computation. Specifically, $s$ is the vector

$\displaystyle s$ $\displaystyle=e^{\frac{1}{2}i(2\alpha-\pi-\theta)}$ (2)
$\displaystyle=\Bigg{(}\cos{\left(\alpha-\frac{1}{2}\pi-\frac{1}{2}\theta\right)},\sin{\left(\alpha-\frac{1}{2}\pi-\frac{1}{2}\theta\right)},0,0,\ldots,0\Bigg{)}.$

To compute $m_{\ell}$ in Line 5, we compute ${\texttt{Indeg}}(v,s)$ then
subtract the cardinality of the set $S:=\\{b\in bigedges~{}|~{}\measuredangle
b<\pi+\alpha\\}$, where $\measuredangle b$ is taken to mean the angle $b$
makes with the $e_{1}$-axis, when viewed as a vector with
${\mathit{EA}}.{\mathit{v}}$ as the origin. By Lemma 2, we compute
${\texttt{Indeg}}(v,s)$ via the oracle using one diagram and $\Theta(\log
n+\Pi)$ time. Since $\mathit{bigedges}$ is sorted and since $s$ lies in the
$(e_{1},e_{2})$-plane, we can find the set $S$ in
$\Theta(\log(|\mathit{bigedges}|))$ time. The subtraction in Line 5 takes
constant time, as does Line 6.

In Lines 7 and 8, we create two edge arc objects. The time complexity of
creating them is proportional to the size of the obejcts themselves. All
attributes of edge arc objects, except the array of vertices
(${\mathit{verts}}$), are constant size. By construction,
${\mathit{EA}}_{\ell}.{\mathit{verts}}$ and
${\mathit{EA}}_{r}.{\mathit{verts}}$ split ${\mathit{EA}}.{\mathit{verts}}$
into two sets, which can be done naïvely in
$\Theta(d|{\mathit{EA}}.{\mathit{verts}}|)$ time by walking through
${\mathit{EA}}.{\mathit{verts}}$ and storing each one explicitly. However, we
improve this to $\Theta(\log|{\mathit{EA}}.{\mathit{verts}}|)$ time if we have
a globally accessible array of vertices (sorted cw around $v$) and just
computes the pointers to the beginning and end of the sub-arrays corresponding
to the ${\mathit{verts}}$ attributes of the new edge arc objects. In total,
Algorithm 1 and takes $\Theta(d+\log
n+\Pi+\log(|\mathit{bigedges}|)+1+\log(|{\mathit{EA}}.{\mathit{verts}}|))=\Theta(\log
n+d+\Pi)$ time and uses uses one diagram.

Now that we have walked through the algorithm and established the runtime and
diagram complexity, we prove correctness. To do so, we first show that
${\mathit{EA}}_{\ell}$ and ${\mathit{EA}}_{r}$ are edge arc objects. In
particular, this means showing that they have the correct values for
${\mathit{count}}$ and ${\mathit{verts}}$. We prove this for
${\mathit{EA}}_{\ell}$; the proof for ${\mathit{EA}}_{r}$ follows a similar
argument.

_${\mathit{EA}}_{\ell}.{\mathit{count}}$ :_ We must show that
${\mathit{EA}}_{\ell}.{\mathit{count}}$ is the number of edges in
${\mathit{EA}}_{\ell}$ incident to ${\mathit{EA}}_{\ell}.{\mathit{v}}$. By
Lemma 1, the value returned from
${\texttt{Indeg}}({\mathit{EA}}.{\mathit{v}},s)$ counts all edges incident to
${\mathit{EA}}.{\mathit{v}}$ and below $s\cdot{\mathit{EA}}.{\mathit{v}}$ in
direction $s$. By Lemma 3, this is exactly the total number of edges in
${\mathit{EA}}_{\ell}$ plus edges
$({\mathit{EA}}.{\mathit{v}},v^{\prime})\in{E}_{v}$ for which $s\cdot
v^{\prime}<s\cdot{\mathit{EA}}.{\mathit{v}}$. Thus, by subtracting $|\\{b\in
bigedges~{}|~{}\measuredangle b<\pi+\alpha\\}|$ from
${\texttt{Indeg}}({\mathit{EA}}.{\mathit{v}},s)$ on Line 5, we are left with
$m_{\ell}$, the number of edges incident to ${\mathit{EA}}.{\mathit{v}}$
contained in ${\mathit{EA}}_{\ell}$. Setting
${\mathit{EA}}_{\ell}.{\mathit{count}}=m_{\ell}$ on Line 7, we see that
${\mathit{EA}}_{\ell}.{\mathit{count}}$ is correct.

_${\mathit{EA}}_{\ell}.{\mathit{verts}}$ :_ We must show that
${\mathit{EA}}_{\ell}.{\mathit{verts}}$ contains an array of all verices
contained in ${\mathit{EA}}_{\ell}$ radially ordered clockwise. This follows
from the fact that ${\mathit{EA}}.{\mathit{verts}}$ is all vertices contained
in ${\mathit{EA}}$ ordered clockwise, so when we restrict
${\mathit{EA}}.{\mathit{verts}}$ to ${\mathit{EA}}.{\mathit{verts}}[:mid]$ on
Line 7, we are eliminating vertices not contained in ${\mathit{EA}}_{\ell}$,
so ${\mathit{EA}}_{\ell}.{\mathit{verts}}$ is correct.

Next, we prove Statement i. Recall that ${\mathit{EA}}.{\mathit{verts}}$
orders the vertices in decreasing angle with $e_{1}$. In Line 3, the angle
$\measuredangle\pi({\mathit{EA}}.{\mathit{verts}}[mid]-{\mathit{EA}}.{\mathit{v}})$
is the angle made by ${\mathit{EA}}.{\mathit{v}}$ with the middle vertex. We
tilt this angle by $\theta/2$ on Line 3 to obtain angle $\alpha$. By
construction of $\alpha$,

$\measuredangle\pi({\mathit{EA}}.{\mathit{verts}}[mid]-{\mathit{EA}}.{\mathit{v}})>\alpha.$

By definition of $\theta$, the angle $\alpha$ satisfies:

$\alpha>\measuredangle\pi({\mathit{EA}}.{\mathit{verts}}[mid+1]-{\mathit{EA}}.{\mathit{v}}).$

Since the array ${\mathit{EA}}.{\mathit{verts}}$ is sorted, all vectors in the
set $\pi({\mathit{EA}}.{\mathit{verts}}[:mid]-{\mathit{EA}}.{\mathit{v}})$
have an angle of at least $\alpha$ with $e_{1}$ and all vectors in
$\pi({\mathit{EA}}.{\mathit{verts}}[:mid]-{\mathit{EA}}.{\mathit{v}})$ have an
angle of at most $\alpha$.

By Lines 1–2 and Line 5, we know that ${\mathit{EA}}.{\mathit{verts}}$
contains the first
$m=\lceil{\frac{1}{2}|{\mathit{EA}}.{\mathit{verts}}|}\rceil$ vertices in
${\mathit{EA}}.{\mathit{verts}}$. Hence, Statement ii holds. Statement iii
follows from Statements i and ii. ∎

###  A.2 Proof of Theorem 2

See 2

###### Proof.

First, we analyze the time complexity of the algorithm and the number of
diagrams it requires. By Lemma 2, Line 1 can be computed in $\theta(\log n)$
time (since we are given the diagram and do not need an additional oracle
query). Storing $A.{\mathit{verts}}$ by storing a pointer to $V_{v}$, we
initialize $\mathit{eastack}$ and $E_{v}$ in Lines 2 and 3 in constant time.

To analyze the complexity of the loop in Lines 4–16, we first note that this
is a radial binary multi-search. When processing an edge arc, we decide
whether all edges have been found or if we need to split the edge arc. If
there is only one edge in the arc (i.e., ${\mathit{EA}}.{\mathit{count}}=1$),
then this loop is a binary search for an edge, using the angle with $e_{1}$ in
the $(e_{1},e_{2})$-plane as the search key. When
${\mathit{EA}}.{\mathit{count}}>1$, we search for all edges, finding them in
decreasing angle order (since arcs with larger angles are added after arcs of
smaller angles). The else-if statement in Lines 8–11 is where the edges are
added to ${E}_{v}$. Note that this shortcuts additional edge arc splitting by
stopping the process once we find that the number of edges in the arc is equal
to the number of potential vertices that can form the edges. As a result, each
edge above $v$ contributes to $O(\log n)$ edge arcs being added to
$\mathit{eastack}$ and, in the case that every other vertex is incident to an
edge with $v$, we have $\Theta(\log n)$ edge arcs added to the stack. All
operations in the while loop are constant time, except splitting the edge arc
object in Line 12, which uses one diagram and takes $\Theta(\log n+d+\Pi)$
time.

The complexity of Algorithm 2 is dominated by the complexity of the while
loop: the algorithm uses $\Theta(\deg(v)\log n)$ augmented persistence
diagrams and takes takes $\Theta((\deg(v)\log n)(\log n+d+\Pi))$ time.

To prove correctness of this algorithm, we track how edges are stored. In this
algorithm, $\mathit{eastack}$ is a stack of edge arc objects. We abuse
notation slightly and say that an edge $(v^{\prime},v)\in E$ is in
$\mathit{eastack}$ if there exists an edge arc $A\in\mathit{eastack}$ such
that $\measuredangle(v^{\prime},v)\in(A.\alpha_{1},A.\alpha_{2})$. In this
case, we also say $(v^{\prime},v)$ is in $A$. We claim the while loop has the
following loop invariant:

  1. (i)

every outgoing edge from $v$ is either in $\mathit{eastack}$ or is in the list
of edges $E_{v}$;

  2. (ii)

the edges of $E_{v}$ are radially sorted clockwise order;

  3. (iii)

the edge arcs in $\mathit{eastack}$ are non-overlapping and in radially
clockwise order; and,

  4. (iv)

for all $(v^{\prime},v)\in E_{v}$ and
$(v^{\prime\prime},v)\in\mathit{eastack}$, the angle
$\measuredangle(v^{\prime},v)$ is larger than the angle
$\measuredangle(v^{\prime\prime},v)$.

Entering the loop, $E_{v}$ is empty and $\mathit{eastack}$ is a single edge
arc accounting for the entire half space above $v$. Thus, all outgoing edges
are in $\mathit{eastack}$ and $E_{v}$ is empty, so the loop invariant is
initially true. Denote the edge arc that is processed in iteration $j$ by
${\mathit{EA}}^{j}$. Suppose entering iteration $j$, the loop invariant holds.
In Line 5, we pop ${\mathit{EA}}^{j}$ from $\mathit{eastack}$. There are three
cases to consider. First, if ${\mathit{EA}}^{j}.{\mathit{count}}=0$, there are
no edges in ${\mathit{EA}}^{j}$, so we enter the if clause in Line 6 and
continue to the top of the while loop. Since no edges were removed from
$\mathit{eastack}$ when ${\mathit{EA}}^{j}$ was popped, we do not add any
edges to $E_{v}$ nor do we change any ordering of the previously sorted edges
or edge arcs in $\mathit{eastack}$. Thus, the loop invariant is maintained.
Second, if
${\mathit{EA}}^{j}.{\mathit{count}}=|{\mathit{EA}}^{j}.{\mathit{verts}}|$, we
enter the else-if clause in Line 8. Here, we know the vertices in the current
edge arc are in a one-to-one correspondence with the edges in
${\mathit{EA}}^{j}$, so we add these edges in sorted order to $E_{v}$ in Line
9; since the loop iteration was true upon entering the loop, these edges are
between the edges already in $E_{v}$ and the remaining edges in
$\mathit{eastack}$. Thus, iii and iv are maintained. Furthermore, since we add
the edges of ${\mathit{EA}}^{j}$ in radially clockwise order, $E_{v}$ is still
radially sorted in clockwise order, satisfying ii. Since all edges of
${\mathit{EA}}^{j}$ are now in $E_{v}$, i is satisfied. So again, the loop
invariant is maintained. The third and final case is that
${\mathit{EA}}.{\mathit{count}}\neq 0$ and
${\mathit{EA}}^{j}.{\mathit{count}}\neq|{\mathit{EA}}^{j}.{\mathit{verts}}|$.
Then, we split the current edge arc on Line 12, pushing the two new halves
onto the stack in a radially clockwise order, maintaining iii. In this case,
$E_{v}$ remains as it was going into the loop and the edges in
$\mathit{eastack}$ remain the same, so the loop invariant is maintained. To
finalize the proof of partial correctness, suppose the loop ends and that the
loop invariant is true. Since the loop ends, $\mathit{eastack}$ is empty. By i
and ii, all outgoing edges are in $E_{v}$ in radially clockwise order. The
runtime analysis tells us that the algorithm terminates, and thus, the
algorithm is correct. ∎

###  A.3 Proof of Theorem 3

See 3

###### Proof.

We first analyze the runtime and diagram count for Algorithm 3 by walking
through the algorithm line-by-line. In Line 1, we ask the oracle for the
diagram in direction $-e_{2}$, which takes $\Theta(\Pi)$ time. In [2, Theorem
14 (Edge Reconstruction)], simultaneously find the cyclic ordering around all
vertices in $\Theta(n^{2})$ time by Lemmas $1$ and $2$ of [13]. In Line 3, we
do that as well; however, we do not store vertices that are above $v$ in the
array $\mathit{vertsabove}[v]$, and thus this line takes $\Theta(n^{2})$ time.
We note that such a cyclic ordering exists around each vertex by Assumption
1projectedindep. Once we have $\mathit{vertsabove}$, to find the minimum angle
defined by any three vertices of $V$, we check all angles between vectors
$\mathit{vertsabove}[v][i]-v$ and $\mathit{vertsabove}[v][i+1]-v$ in Line 4 in
$\Theta(n+m)$ time.

The for loop in Lines 5–8 is repeated $n$ times, once for each vertex in $V$.
To determine the order of processing the vertices in $V$, we follow the births
in $\widehat{\mathcal{D}}_{0}$, in decreasing order (since
$\widehat{\mathcal{D}}_{0}$ is the lower-star filtration in direction
$-e_{2}$). Thus, finding the order takes $\Theta(n)$ time. In each iteration,
we compute the incoming edges (those whose other vertex is below $v$) in Line
6 followed by all outgoing edges (those whose other vertex is above $v$) in
Line 7. By Assumption 1iii, every edge is either incoming or outgoing with
respect to direction $e_{2}$. Thus, all edges incident to $v$ are in $E$ once
$E$ is updated in Line 6. By Theorem 2, when processing vertex $v$, the call
to Algorithm 2 on Line 7 takes $\Theta((\deg(v)\log n)(\log n+d+\Pi))$ time
and uses $\Theta(\deg(v)\log n)$ diagrams. Summing over all vertices, we see
that the loop in Lines 5–8 takes

$\displaystyle\sum_{v\in V}\Theta((\deg(v)\log n)(\log n+d+\Pi))$
$\displaystyle=\Theta(m\log n(\log n+d+\Pi))$

time and uses $\sum_{v\in V}\Theta(\deg(v)\log n)=\Theta(m\log n)$ augmented
persistence diagrams.

In total, Algorithm 3 takes
$\Theta(\Pi)+\Theta(n^{2})+\Theta(n+m)+\Theta(n)+\Theta(m\log n(\log
n+d+\Pi))=\Theta(n^{2}+m\log n(\log n+d+\Pi))$ time and uses
$\Theta(1)+\Theta(m\log n)=\Theta(m\log n)$ diagrams.

Next, we prove the correctness of Algorithm 3 (i.e., that all edges are
found). In order to process vertices in order of their heights in the $e_{2}$
direction, we first sort them in Line 5. For $1\leq j\leq n$, let $v_{j}$ be
the $j$th vertex in this ordering. To show that Algorithm 3 finds all edges in
$E$, we consider the loop invariant (LI): when we process $v_{j}$, all edges
with maximum vertex height equal or less than the height of $v_{j}$ are known.
The LI is trivially true for $v_{1}$. We now assume that it is true for
iteration $j$, and show that it must be true for iteration $j+1$. By
assumption, all edges $(v_{i},v_{j})$ with $1\leq i<j$ are known, and so by
Theorem 2, Algorithm 2 finds all edges $(v_{k},v_{j})$, where $k>j$, and we
add them to the edge set $E$. Note that, by assumption, all edges
$(v_{x},v_{i})$ for $1\leq i\leq j$ are also already known, and so the
invariant is maintained. Thus, after the loop terminates, all edges are found.
∎

##  Appendix B Basis

In Assumption 1iii, we assume all vertices of the underlying graph are unique
with respect to the first basis direction $e_{2}$. In this appendix, we
provide details of how to _find_ a basis where all vertices have a unique
height with respect to the second basis direction. I.e., this appendix allows
us to remove one general position assumption by showing it can be satisfied
deterministically, at an added cost of $\Theta(|P|\log|P|+d+\Pi)$ time.

######  Lemma 4 (Creation of Orthonormal Basis).

Given a point set $P\subset\mathbb{R}^{d}$ satisfying Assumption 1i and
Assumption 1ii, we can use two diagrams and $\Theta(|P|\log|P|+d+\Pi)$ time to
create the orthonormal basis $\\{b_{1},b_{2},e_{3},e_{4},\ldots,e_{d}\\}$ so
that all points of $P$ have a unique height in direction $b_{2}$.

###### Proof.

Algorithm 6 (Tilt) of [7] takes diagrams from two linearly independent
directions $s,s^{\prime}\in\mathbb{S}^{d-1}$, the point set $P$, and returns a
direction $s_{*}$ in $\Theta(|P|\log|P|+d+\Pi)$ time444While [7] does not
account for diagram computation time, there are two diagrams used in this
process, hence our addition of $\Theta(\Pi)$ to the total runtime. so that the
following properties holds for all $p_{1},p_{2}\in P$:

  1. (i)

If $p_{1}$ is strictly above (below) $p_{2}$ with respect to direction $s$,
then $p_{1}$ is strictly above (below, respectively) $p_{2}$ with respect to
direction $s_{*}$.

  2. (ii)

If $p_{1}$ and $p_{2}$ are at the same height with respect to direction $s$
and $p_{1}$ is strictly above (below) $p_{2}$ with respect to direction
$s^{\prime}$, then $p_{1}$ is strictly above (respectively, below) $p_{2}$
with respect to direction $s_{*}$.

  3. (iii)

If $p_{1}$ is is at the same height as $p_{2}$ with respect to both directions
$s$ and $s^{\prime}$, then $p_{1}$ and $p_{2}$ are at the same height with
respect to direction $s_{*}$.

A proof of correctness is given in [7, Lemma 32 (Tilt)].

We start with the standard basis for $\mathbb{R}^{d}$,
$\\{e_{1},e_{2},\ldots,e_{d}\\}$, and we replace the first two basis elements
as follows. Let $b_{2}$ be the direction obtain by using Tilt with $s=e_{1}$,
$s^{\prime}=e_{2}$, and $P=P$.

By Assumption 1ii, no three points of $P$ are colinear when projected onto the
first two coordinates. In particular, this means no two points share the same
heights in both the $e_{1}$ and $e_{2}$ directions. Then, by Statements i-ii
above, the direction $b_{2}$ must order all vertices of $P$ uniquely. Using
only the first two coordinates of $b_{2}$ and $e_{1}$, we then perform Gram
Schmidt orthanormalization to compute the first two coordinates of $b_{1}$.
More precisely, letting $b_{i}^{(j)}$ denote the $j$th coordinate of $b_{i}$,
we compute

$\begin{pmatrix}b_{1}^{(1)}\\\ b_{1}^{(2)}\end{pmatrix}=\begin{pmatrix}1\\\
0\end{pmatrix}-\frac{\bigg{\langle}\begin{pmatrix}b_{2}^{(1)}&b_{2}^{(2)}\end{pmatrix}^{T},\begin{pmatrix}1&0\end{pmatrix}^{T}\bigg{\rangle}}{\bigg{|}\bigg{|}\begin{pmatrix}b_{2}^{(1)}&b_{2}^{(2)}\end{pmatrix}^{T}\bigg{|}\bigg{|}^{2}}\begin{pmatrix}b_{2}^{(1)}\\\
b_{2}^{(2)}\end{pmatrix}$ (3)

We then set $b_{1}^{(j)}=0$ for $2<j\leq d$, so that
$b_{1}\in\text{span}\\{e_{1},e_{2}\\}$, $b_{2}\perp b_{1}$, and $||b_{1}||=1$.
Only considering the first two coordinates of $b_{2}$ and $e_{1}$ means this
process takes constant time. The remaining $e_{i}$ for $2\leq i\leq d$ can be
used to fill the basis. Finally, we have a basis satisfying all assumptions of
Assumption 1, namely, $\\{b_{1},b_{2},e_{3},e_{4},\ldots,e_{d}\\}$. ∎