# Deep Surface Reconstruction from Point Clouds with Visibility Information

Raphael Sulzer1, Loïc Landrieu1, Alexandre Boulch3, Renaud Marlet23 and Bruno
Vallet1  1LASTIG, Univ Gustave Eiffel, IGN-ENSG, F-94160 Saint-Mandé, France
Email<EMAIL_ADDRESS>2LIGM, Ecole des Ponts, Univ Gustave Eiffel,
CNRS, Marne-la-Vallée, France  3Valeo.ai, Paris, France

###### Abstract

Most current neural networks for reconstructing surfaces from point clouds
ignore sensor poses and only operate on raw point locations. Sensor
visibility, however, holds meaningful information regarding space occupancy
and surface orientation. In this paper, we present two simple ways to augment
raw point clouds with visibility information, so it can directly be leveraged
by surface reconstruction networks with minimal adaptation. Our proposed
modifications consistently improve the accuracy of generated surfaces as well
as the generalization ability of the networks to unseen shape domains.

##  I Introduction

The problem of reconstructing a watertight surface from a point cloud has
recently been addressed by a variety of deep learning based methods. Compared
to traditional approaches, deep surface reconstruction (DSR) can learn shape
priors [1, 2] and leverage shape similarities [3] to complete missing parts
[4], filter outliers, or smoothen noise in defect-laden point clouds. DSR
methods, however, often derive priors from training datasets with few shape
classes, generalizing poorly to unseen categories or datasets. Learning more
local priors improves consistency across different objects or scenes [5, 6]
but may result in higher sensitivity to noise or other defects. Besides, lack
of global context complicates surface orientation.

For real world point clouds, usually acquired via active or passive methods
such as LiDAR scanning or multi-view stereo (MVS), the sensor position can be
known and used to relate each observed point with a line of sight. Such
visibility information can then help to orient surface normals [7] or predict
occupancy [8, 9, 10]. While visibility is key for MVS, it has largely been
ignored by DSR methods. In fact, sensor positions are usually not given in
reconstruction benchmarks from point clouds.

|
---|---
(a) Reconstruction using only the points position.

|
---|---
(b) Reconstruction with visibility augmented point cloud.

Figure 1: Surface Reconstruction with Visibility Information. We augment each
3D point  with a sightline vector  pointing towards the sensor observing it.
Additionally, two auxiliary points are placed before  and after  the observed
point along the sightline. This allows DSR networks, with very little
modification, to reconstruct a significantly more accurate surface.

To remedy this, we consider virtual scanning rather than uniform sampling, and
we show that many DSR methods can easily be adapted to benefit from visibility
(cf. Figure 1). Our main contributions are as follows:

  * •

We propose two simple ways to add visibility information to 3D point clouds,
and we detail how to adapt DSR methods to utilize them, with very little
changes.

  * •

Using synthetic and real data, at object and scene level, we show for a wide
range of state-of-the-art DSR methods that models leveraging visibility
reconstruct higher-quality surfaces and are more robust to domain shifts.

Incidentally, our benchmarks also allow us to compare a range of recent state-
of-the-art DSR methods on the same ground.

##  II Related Work

Many traditional surface reconstruction methods use visibility information [8,
9, 10, 11, 12, 13, 14]. They are usually based on a 3D Delaunay
tetrahedralization, which is intersected with lines of sight to attribute
visibility features to Delaunay cells. While such methods can scale to
billions of points [15] and are robust to moderate levels of noise and
outliers, they do not incorporate learned shape priors.

In contrast, recent DSR methods have shown to produce more accurate surfaces
than traditional approaches, especially for shape categories encountered
during training. Many DSR methods use an implicit surface representation,
either based on occupancy [2, 16, 17], or on the distance to the surface,
whether it is signed [1, 17, 18, 19] or unsigned [20, 21, 22, 23]. To
integrate local information, different forms of convolutions are used, either
on regular grids [4, 5, 24, 25, 26], directly on points [27, 28] or via an MLP
instead [29]. Other methods rather use an explicit surface representation such
as a mesh, which is deformed [3] or whose elements classified [6, 30].

Figure 2:  Visibility-Augmented Point Cloud. Each observed point is associated
to a sightline unit vector pointing towards its sensor. Two new points before
and after each point are added. They help to disambiguate occupancy.

A key issue is to get a sense of surface point orientation, to choose between
reconstructing a thin volume (two main opposite orientations) or a thicker one
(one main orientation at void-matter interface). Some methods dismiss the
orientation issue by requiring oriented normals as input [5, 27, 31, 32],
albeit producing such normals is a challenging task in itself [7, 33, 34]. We
show that oriented normals can be advantageously replaced by visibility
information.

Only few deep-learning methods make use of visibility information, typically
from multiple views with camera pose information. RayNet [35] aggregates
features from pixels of different views that intersect in the same voxel, but
it outputs a dense point cloud, not a watertight surface mesh. Neural radiance
fields [36, 37] somehow also model the free space between a point and its
sensor. They, however, generally assume numerous and dense views (i.e.,
images), and leverage little or no shape priors. We argue that DGNN [6], that
classifies Delaunay cells with a graph neural network, currently is the only
general DSR method from point clouds with visibility. However, DGNN relies on
handcrafted visibility features requiring substantial geometry processing,
while, we propose to directly augment the input point clouds. For point clouds
for which visibility information is not available, Vis2Mesh [38] shows that
rendering virtual views and learning point sensor visibility can significantly
improve the reconstruction quality of a traditional method.

##  III Method

We consider a 3D point cloud $P$ where each point $p\,{\in}\,P$ has some
coordinates $X_{p}\,{\in}\,\mathbb{R}^{3}$ and knows the position
$S_{p}\,{\in}\ \mathbb{R}^{3}$ of a sensor observing it. Instead of only using
the raw point coordinates $(X_{p})_{p\in P}$ as the input $(I_{p})_{p\in P}$
of a DSR network, we propose two simple ways to augment point cloud $P$ with
visibility information, and adapt DSR methods accordingly.

###  III-A Sightline Vector (SV)

For each $p\in P$, we define a unit vector $\mathbf{v}_{p}$ pointing from the
observation $X_{p}$ to the sensor $S_{p}$:
$\mathbf{v}_{p}\,{=}\,{(S_{p}\,{-}\,X_{p})}/{\|S_{p}\,{-}\,X_{p}\|}$. This
contains useful information for surface orientation. We normalize the vector
as the distance to the sensor is not as relevant as the viewing angle.

###  III-B Auxiliary Points (AP)

To help the network predict empty and full space immediately in front of and
behind the observed surface, we consider two auxiliary points to each point
$p$: a _before-point_ $p_{\mathrm{b}}$ and an _after-point_ $p_{\mathrm{a}}$,
located along the sightline on each side of $p$:
$X_{p_{\mathrm{b}}}=X_{p}+d\mathbf{v}_{p},X_{p_{\mathrm{a}}}=X_{p}-d\mathbf{v}_{p}$,
where $d$ is a characteristic distance in the point cloud $P$, e.g., the
average distance from a point to its nearest neighbor. By construction,
$p_{\mathrm{b}}$ is likely outside the scanned object or scene (modulo sensing
noise and outliers), and ${p_{\mathrm{a}}}$, likely inside (modulo object
thickness too).

###  III-C Visibility-Augmented Point Cloud

We use sightline vectors and auxiliary points to add visibility information to
an input point cloud, separately or together.

  * (SV) To use sightline information only, we simply concatenate the sightline vector channelwise to the point coordinates to form the network input: $I_{p}=(X_{p}\oplus\mathbf{v}_{p})\in\mathbb{R}^{6}$.

  * (AP) To use auxiliary points only, we add before-points $p_{\mathrm{b}}$ and after-points $p_{\mathrm{a}}$ to $P$, with tags $\mathbf{t}\in\mathbb{R}^{2}$ concatenated to point coordinates to identify the point type, i.e., $I_{q}\,{=}\,(X_{q}\,{\oplus}\,\mathbf{t}_{q})\in\mathbb{R}^{5}$ with $q\in\\{p,p_{\mathrm{b}},p_{\mathrm{a}}\\}$, where $\mathbf{t}_{p}\,{=}\,[0\,\,0]$ (observed point), $\mathbf{t}_{p_{\mathrm{b}}}\,{=}\,[1\,\,0]$ (before-point), or $\mathbf{t}_{p_{\mathrm{a}}}\,{=}\,[0\,\,1]$ (after-point).

  * (SV+AP) When combining both kinds of visibility information, before-points $p_{\mathrm{b}}$ and after-points $p_{\mathrm{a}}$ are given the same sightline vector as their reference point, i.e., $\mathbf{v}_{p_{\mathrm{b}}}=\mathbf{v}_{p_{\mathrm{a}}}=\mathbf{v}_{p}$, and we take as input $I_{p}=(X_{p}\oplus\mathbf{v}_{p}\oplus\mathbf{t}_{p})\in\mathbb{R}^{8}$.

While holding a similar kind of information, no augmentation can be reduced to
the other one. SVs alone are not enough to place APs, and APs alone, as they
are not associated to their observed point in $P$, cannot determine SVs (cf.
Figure 2).

###  III-D Modifying an Existing Architecture

We can adapt most DSR networks to handle visibility-augmented point clouds
with only few modifications:

  * •

We change the input size (number of channels) of the first layer of the
network (generally an encoder), increasing it by 2, 3 or 5, depending on the
augmentation.

  * •

We directly add auxiliary points to the point cloud, thus tripling the number
of input points. For methods based on neighboring point sampling, we add
auxiliary points after sampling for more efficiency.

The batch size may need to be adjusted to fit a larger point cloud in memory,
but the rest of the network stays unchanged. Its size is mostly unaltered
(e.g., +0.005% for ConvONet [25]).

Bare |  |  |  |  |  |
---|---|---|---|---|---|---
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
| Input |  ConvONet-2D [25] |  Points2Surf [29] |  Shape As Points [39] |  POCO [28] | Ground Truth
Figure 3:  Object-Level Reconstruction. Reconstructed shapes from the
ModelNet10 test set using four different DSR methods trained on ModelNet10.
Top rows of each object use the bare point cloud as input, and bottom rows use
the point cloud augmented with visibility information.

##  IV Experiments

To assess our proposal, we first detail our simple adaptation of six different
DSR baseline networks to leverage our visibility information, then compare the
quality of the reconstructed surfaces and analyze the generalization ability
of the networks.

###  IV-A DSR Baselines

####  IV-A1 ConvONet [25]

This method first extracts point features and projects them on three 2D grids,
or one 3D grid (variant). 2D or 3D grid convolutions then create features
capturing local occupancy. Last, the occupancy of a query-point is estimated
after interpolating grid features. We consider the $3\times 64^{2}$ 2D-plane
encoder and the $64^{3}$ 3D-volume variant. To adapt them, we change the input
size of the point encoder’s first layer.

####  IV-A2 Points2Surf [29]

This method predicts both the occupancy of a query point and its unsigned
distance to the surface. It uses both a local query-point neighborhood and a
global point-cloud sampling. We use the best-performing variant (uniform
global sampling, no spatial transformer). To adapt it, we increase the input
size of the first layer of both the local and global encoders, and when a
point is sampled, locally or globally, we add its two auxiliary points on the
fly.

####  IV-A3 Shape As Points [39]

For each input point, the method estimates its normal as well as $k$ point
offsets that are used to correct and densify the point cloud. The resulting
point cloud of size $k|P|$ is then fed to a differentiable Poisson solver
[31]. To adapt the method, we change the input size of the first layer of the
encoder, and of the normal and offset decoders as they also input the point
cloud. We directly add auxiliary points as input, whose normal and offsets
will thus be computed too.

####  IV-A4 Local Implicit Grids (LIG) [5]

This method trains an auto-encoder from dense point cloud patches. For
inference, a given sparse patch with oriented normals is first augmented,
close to our idea, with 10 new points along each normal; then reconstruction
uses latent vectors minimizing a decoder-based training loss, and a post-
processing removes falsely-enclosed volumes. As training code is unavailable,
we use the model pretrained on ShapeNet (without noise). For oriented normals,
we use Jets [40] oriented with a minimum spanning tree [7], as in [23]. To
exploit visibility, we replace normals with sightline vectors; we do not add
(more) auxiliary points.

####  IV-A5 POCO [28]

This method extracts point features using point cloud convolution [41], then
estimates the occupancy of a query point with a learning-based interpolation
on nearest neighbors. To adapt it, we increase the input size of the first
layer and add auxiliary points on the fly only in the first layer.

####  IV-A6 DGNN [6]

This method uses a graph neural network to estimate the occupancy of Delaunay
cells in a point cloud tetrahedralization. A graph-cut-based optimization then
reinforces global consistency. The method, which already uses visibility,
outperforms other traditional reconstruction methods that use visibility
information. As it already exploits visibility, we do not alter it, but use it
as baseline for comparison.

For all methods, unless otherwise stated, training and evaluation are
unchanged; we keep the value of the hyperparameters used in the original
papers. When Marching cubes [42] are needed for surface extraction, we use a
grid resolution of $128^{3}$.

###  IV-B Datasets

We consider a variety of object and scene datasets, both synthetic and real,
to show the versatility of our approach.

####  IV-B1 ModelNet10

We use the official train/test splits of all 10 object classes of ModelNet10
[43]. We hold out 10% of the train set for validation. We synthetically scan
the models by placing $10$ random range scanners in two bounding spheres
around the objects and shooting random rays to a sphere inscribed within the
convex hull of the object. We sample $3\,000$ points per object and add
Gaussian noise with zero mean and standard deviation $0.005$ as in [25].

####  IV-B2 ShapeNet

We study the generalizability of models trained on ModelNet10 by testing on
100 shapes per class from the ShapeNet [44] test set of Choy et al. [45] ($9$
out of $13$ classes are not represented in ModelNet10). We use the same
scanning procedure as for ModelNet10.

####  IV-B3 Synthetic Room

We use the train/val/test splits of Synthetic Rooms [25]. For virtual
scanning, we only place sensors in the upper hemispheres, and scan $10\,000$
points as in [25].

####  IV-B4 SceneNet

We test on a few synthetic scenes of SceneNet [46] using the given virtual
scans, voxel-decimated to 1 cm.

####  IV-B5 ScanNet

We test on a few real scenes of ScanNet [47] using the provided real RGB-D
scans, voxel-decimated to 2 cm.

####  IV-B6 Tanks and Temples

We use the real LiDAR point cloud of the _Ignatius_ statue from the Tanks and
Temples dataset [48].

####  IV-B7 Middlebury

We use an MVS point cloud of the _TempleRing_ from Middlebury [49], made with
OpenMVS [50].

###  IV-C Metrics

We report volumetric intersection over union (IoU), mean Chamfer distance
$\times 100$ (CD) and normal consistency (NC).

###  IV-D Ablation Study

TABLE I:  Ablation Study.

The vanilla model of ConvONet trained and tested on ModelNet10 with different
ways to add visibility or normal information.

Model | SV | AP | IoU $\uparrow$
---|---|---|---
ConvONet-2D ($3\times 64^{2}$) [25] |  |  | 0.853
\+ sightline vectors (SV) only | ✓ |  | 0.871
\+ auxiliairy points (AP) only |  | ✓ | 0.881
\+ both SV and AP | ✓ | ✓ | 0.886
\+ sensor position | $S_{p}$ |  | 0.870
\+ unnormalized SV | $S_{p}-X_{p}$ |  | 0.870
\+ estim. normals / estim. orientation | Jets [40] / MST [7] |  | 0.853
\+ estim. normals / sensor orientation | Jets [40] / sensor-based [7] |  | 0.868
\+ true normals | GT normals |  | 0.879

To validate our design, we compare in Table I various ways to add visibility
information to the vanilla model of ConvONet.

Independently, SVs and APs significantly improve performance (+$1.8$ and
+$2.8$ IoU pts). A reason why APs are more profitable could be that the
network is tailored for points, not points with sightline features. While SVs
and APs capture a similar kind of information, they are, however,
complementary: combining them is even more beneficial (+$3.5$ IoU pts). Our
general interpretation is that SVs help to decide whether a locally “thin”
point cloud is to be considered as a noisy scan of a single surface, or as a
(less noisy) scan on both sides of a thin surface. They thus have an impact on
local shape topology, which can bring a notable gain. Auxiliary points convey
similar information, but also contribute more directly to refine the surface
position. Replacing SVs by the sensor position or by the unnormalized point-
sensor vector gives essentially the same performance than our unit vector.
This can be explained by the fact that our scanning procedure does not
introduce significant variation in terms of distance to the sensor. Yet, for
real world acquistions, with a larger range of sensor distances, normalizing
the SV ensures more stability.

Adding SVs outperforms estimated normals [40] with estimated orientation [7],
and even estimated normals with sensor-based orientation. While using ground-
truth normals is slightly more beneficial than SVs, combining SV+AP yields the
best overall performance, which highlights the richness of our visibility
information.

We also experiment with adding more than two auxiliary points: (i) at distance
$0.5d$ or $2d$, (ii) at the midpoint between sensor and point, or (iii) as
grazing points, estimated by densely sampling the sightlines with auxiliary
points and keeping the ones close to an input point. None of these strategies
brought significant improvements over adding two points at distance $d$ on
both sides of the real point.

TABLE II:  Object-Level Reconstruction.

DSR methods trained and tested on ModelNet10, with and without sightline
vectors (SV) or auxiliary points (AP). † Trained on ShapeNet.

Model | SV | AP | IoU $\uparrow$ | CD $\downarrow$ | NC $\uparrow$
---|---|---|---|---|---
ConvONet-2D [25] |  |  | 0.853 | 0.618 | 0.934
ConvONet-2D [25] | ✓ |  | 0.871 | 0.557 | 0.936
ConvONet-2D [25] | ✓ | ✓ | 0.886 | 0.518 | 0.943
ConvONet-3D [25] |  |  | 0.885 | 0.493 | 0.949
ConvONet-3D [25] | ✓ |  | 0.911 | 0.424 | 0.956
ConvONet-3D [25] | ✓ | ✓ | 0.913 | 0.423 | 0.957
Points2Surf [29] |  |  | 0.842 | 0.590 | 0.890
Points2Surf [29] | ✓ |  | 0.859 | 0.544 | 0.896
Points2Surf [29] | ✓ | ✓ | 0.856 | 0.548 | 0.897
Shape As Points [39] |  |  | 0.903 | 0.438 | 0.948
Shape As Points [39] | ✓ |  | 0.907 | 0.430 | 0.950
Shape As Points [39] | ✓ | ✓ | 0.913 | 0.414 | 0.953
POCO [28] |  |  | 0.907 | 0.422 | 0.945
POCO [28] | ✓ |  | 0.915 | 0.408 | 0.950
POCO [28] | ✓ | ✓ | 0.917 | 0.406 | 0.950
† LIG [5] |  |  | – | 0.974 | 0.849
† LIG [5] | ✓ |  | – | 0.880 | 0.882
DGNN [6] | ✓ |  | 0.866 | 0.543 | 0.884

###  IV-E Object-Level Reconstruction

Table II reports the performance on ModelNet10 of various models, with and
without sightline vectors or auxiliary points.

ConvONet (both planar and volumetric) gains about +$3$ IoU pts with visibility
information. The resulting surface is more accurate, especially in concave
parts, as illustrated in Figure 3.

Points2Surf improves with sightline vectors, but auxiliary points do not
improve further: the sensor vectors are enough to resolve ambiguities for the
occupancy estimation, but distance estimation does not further benefit from
auxiliary points.

Shape As Points benefits from sightline vectors, although not as much as other
methods, probably because the model also estimates normals which provide a
similar information as visibility. Still, adding auxiliary points further
gains +$0.6$ IoU pts, yielding more complete and smoother surfaces.

POCO similarly benefits +$1$ IoU pts from sightline vectors but not much from
the further addition of auxiliary points. While sightline vectors help for
surface orientation, POCO is already accurate enough for APs to bring little
refinement.

LIG produces poor results, likely because the only available model is trained
on ShapeNet, with uniform sampling, little or no noise, and because oriented
normals are only estimated. We cannot report IoU because LIG’s post-processing
creates holes in some objects. Yet, replacing the estimated normals by
sightline vectors improves the predicted surface.

DGNN, which already exploits visibility and outperforms ConvONet-2D and
Points2Surf, is outdistanced on this dataset by methods that use our augmented
point clouds.

Running time is mostly unaffected when adding sightline vectors. The effect of
auxiliary points depends on the method. For ConvONet and POCO, it is
negligible ($<$ +2%); running time is $\times$ 2.25 for Points2Surf and
$\times$ 1.25 for Shape As Points.

TABLE III:  Scene-Level Reconstruction.

ConvONet trained and tested in sliding-window mode on Synthetic Rooms.

Model | SV | AP | IoU $\uparrow$ | CD $\downarrow$ | NC $\uparrow$
---|---|---|---|---|---
ConvONet-3D [25] |  |  | 0.805 | 0.598 | 0.906
ConvONet-3D [25] | ✓ | ✓ | 0.832 | 0.569 | 0.911

TABLE IV:  Out-of-Domain Object-Level Reconstruction.

DSR methods trained on ModelNet10 and tested on ShapeNet, with and without
sightline vectors (SV) or auxiliary points (AP). † Trained on ShapeNet.

Model | SV | AP | IoU $\uparrow$ | CD $\downarrow$ | NC $\uparrow$
---|---|---|---|---|---
† ConvONet-2D [25] |  |  | 0.852 | 0.560 | 0.929
ConvONet-2D [25] |  |  | 0.685 | 0.979 | 0.878
ConvONet-2D [25] | ✓ |  | 0.667 | 1.042 | 0.833
ConvONet-2D [25] | ✓ | ✓ | 0.750 | 0.891 | 0.878
ConvONet-3D [25] |  |  | 0.628 | 0.972 | 0.885
ConvONet-3D [25] | ✓ |  | 0.759 | 0.724 | 0.905
ConvONet-3D [25] | ✓ | ✓ | 0.854 | 0.554 | 0.925
Points2Surf [29] |  |  | 0.807 | 0.561 | 0.876
Points2Surf [29] | ✓ |  | 0.836 | 0.516 | 0.886
Points2Surf [29] | ✓ | ✓ | 0.833 | 0.522 | 0.887
† Shape As Points [39] |  |  | 0.838 | 0.577 | 0.923
Shape As Points [39] |  |  | 0.494 | 0.997 | 0.859
Shape As Points [39] | ✓ |  | 0.749 | 0.843 | 0.881
Shape As Points [39] | ✓ | ✓ | 0.821 | 0.617 | 0.919
POCO [28] |  |  | 0.391 | 1.119 | 0.839
POCO [28] | ✓ |  | 0.832 | 0.618 | 0.901
POCO [28] | ✓ | ✓ | 0.815 | 0.635 | 0.887
DGNN [6] | ✓ |  | 0.844 | 0.549 | 0.854

###  IV-F Scene-Level Reconstruction

To study the impact of visibility information at scene level, we train and
test ConvONet on Synthetic Rooms, in sliding-window mode [25]. We report
quantitative results in Table III and qualitative results in the supplementary
material. The model gains almost +$3$ IoU pts with visibility information,
showing that benefits scale to scenes, not just objects.

###  IV-G Generalization to New Domains

| SceneNet | ScanNet
---|---|---
Input |  |
Bare |  |
Augmented |  |
Figure 4: Out-of-Domain Scene-Level Reconstruction. POCO trained on
ModelNet10, with and without visibility information, is run on scenes from
SceneNet (synthetic RGB-D scan) and ScanNet (real RGB-D scan).

To evaluate the impact of added visibility on the generalization capability of
DSR methods, we train on ModelNet10 and test on ShapeNet (Table IV).

We observe that ConvONet, Shape As Points and POCO trained with visibility
information generalize much better on the new objects and classes, with a gain
up to +44 IoU pts. For comparison, we also show the scores of official models
trained on ShapeNet, although trained on uniformly sampled points rather than
virtual scans, which explains the drop of performance compared to the numbers
in the papers [25, 39]. Points2Surf also improves by up to +3 IoU pts with
added sightline vectors, but not further with APs.

Bare |  |  |  |  |  |
---|---|---|---|---|---|---
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
| Input |  ConvONet-3D [25] |  Points2Surf [29] |  Shape As Points [39] |  POCO [28] | HD Scan
Figure 5:  Out-of-Domain Object-Level Reconstruction. Reconstructed shapes from a LiDAR point cloud (top, _Ignatius_ from Tanks And Temples) and a MVS point cloud (bottom, _TempleRing_ from Middlebury) using four different DSR methods trained on ModelNet10. Top rows of each object use the bare point cloud as input, and bottom rows use the point cloud augmented with visibility information. _HD Scan_ is a high-density point cloud.  |  ConvONet-3D [25] |  POCO [28]
---|---|---
Bare |  |
Augmented |  |
Figure 6: Cut of Out-of-Domain Object Reconstruction. A cut (along the green
curve) of the reconstructed surface of _Ignatius_. The reconstructions from
the bare point cloud (top row) include empty space enclosed inside the object
with backfaces, leading to a poor volumetric IoU. The reconstructions from the
point cloud augmented with visibility information (bottom row) include only
one surface, close to the input points.

The increased generalization capability of the models is also validated when
reconstructing surfaces from real-world scans obtained with LiDAR or MVS. In
Figures 4 and 5, we show that networks using visibility information can
reconstruct a more accurate and more complete surface. The reason for the
largely improved volumetric IoU when using visibility information is
illustrated in Figure 6. For out-of-domain reconstructions, the baseline
methods often predict hollow shapes, i.e., empty space enclosed inside an
object. This leads to backfaces behind the real surface and a poor volumetric
IoU. On the contrary, our models, trained on visibility-augmented point
clouds, learn to distinguish between empty and full space more reliably and do
not produce such artifacts.

##  V Limitations and Perspectives

The position of auxiliary points depends on parameter $d$, which is the
average distance, across the whole scene, from a point to its nearest
neighbor. To better handle point density variations, it could be set locally
rather than globally. Besides, as this positioning is also sensitive to
sampling noise, $d$ could also be directly adjusted after noise estimation.

Our current approach only associates each point with a single sensor, while
MVS points typically have several. A more efficient and versatile approach
than simply duplicating sightlines is still an open issue.

Last, we resort to virtual scans because current 3D reconstruction benchmarks
do not provide sensor positions. While we show that using our augmented point
clouds allows common architectures to successfully generalize from virtual to
real scenes, our training set may fail to replicate some challenging
configurations encountered using actual sensors.

##  VI Conclusion

The sensor poses are often ignored in point cloud processing, even though
available with most acquisition technologies. We present two straightforward
ways to exploit sensor positions to augment point clouds with visibility
information. Our experiments show that various deep surface reconstruction
methods can be adapted with minimal effort to exploit these visibility-
augmented point clouds, resulting in improved accuracy and completeness of
reconstructed surfaces, as well as a substantial increase in generalization
capability.

##### Acknowledgments

This work was partially funded by the ANR-17-CE23-0003 BIOM grant.

## References

  * [1] J. J. Park, P. Florence, J. Straub, R. Newcombe, and S. Lovegrove, “DeepSDF: Learning continuous signed distance functions for shape representation,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * [2] L. Mescheder, M. Oechsle, M. Niemeyer, S. Nowozin, and A. Geiger, “Occupancy networks: Learning 3D reconstruction in function space,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * [3] R. Hanocka, G. Metzer, R. Giryes, and D. Cohen-Or, “Point2Mesh: A self-prior for deformable meshes,” _ACM Transaction on Graphics_ , 2020.
  * [4] A. Dai, C. Diller, and M. Nießner, “SG-NN: Sparse generative neural networks for self-supervised scene completion of RGB-D scans,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * [5] C. M. Jiang, A. Sud, A. Makadia, J. Huang, M. Nießner, and T. Funkhouser, “Local implicit grid representations for 3D scenes,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * [6] R. Sulzer, L. Landrieu, R. Marlet, and B. Vallet, “Scalable surface reconstruction with delaunay‐graph neural networks,” _Eurographics Symposium on Geometry Processing (SGP)_ , 2021.
  * [7] N. Schertler, B. Savchynskyy, and S. Gumhold, “Towards globally optimal normal orientations for large point clouds,” _Computer Graphics Forum (CFG)_ , 2017\.
  * [8] P. Labatut, J. P. Pons, and R. Keriven, “Robust and efficient surface reconstruction from range data,” _Computer Graphics Forum (CGF)_ , 2009.
  * [9] H. H. Vu, P. Labatut, J. P. Pons, and R. Keriven, “High accuracy and visibility-consistent dense multiview stereo,” _IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI)_ , 2012.
  * [10] M. Jancosek and T. Pajdla, “Multi-view reconstruction preserving weakly-supported surfaces,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2011.
  * [11] L. Caraffa, M. Brédif, and B. Vallet, “3D watertight mesh generation with uncertainties from ubiquitous data,” in _Asian Conference on Computer Vision (ACCV)_ , 2017.
  * [12] A. Bódis-Szomorú, H. Riemenschneider, and L. Van Gool, “Efficient volumetric fusion of airborne and street-side data for urban reconstruction,” in _International Conference on Pattern Recognition (ICPR)_ , 2016.
  * [13] M. Jancosek and T. Pajdla, “Exploiting visibility information in surface reconstruction to preserve weakly supported surfaces,” _International Scholarly Research Notices_ , 2014.
  * [14] Y. Zhou, S. Shen, and Z. Hu, “Detail preserved surface reconstruction from point cloud,” _Sensors_ , 2019.
  * [15] L. Caraffa, Y. Marchand, M. Brédif, and B. Vallet, “Efficiently distributed watertight surface reconstruction,” in _3DV_ , 2021.
  * [16] Z. Chen and H. Zhang, “Learning implicit fields for generative shape modeling,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * [17] R. Chabra, J. E. Lenssen, E. Ilg, T. Schmidt, J. Straub, S. Lovegrove, and R. Newcombe, “Deep local shapes: Learning local SDF priors for detailed 3D reconstruction,” in _European Conference on Computer Vision (ECCV)_ , 2020.
  * [18] M. Michalkiewicz, J. Pontes, D. Jack, M. Baktashmotlagh, and A. Eriksson, “Implicit surface representations as layers in neural networks,” in _International Conference on Computer Vision (ICCV)_ , 2019.
  * [19] A. Gropp, L. Yariv, N. Haim, M. Atzmon, and Y. Lipman, “Implicit geometric regularization for learning shapes,” in _International Conference on Machine Learning (ICML)_ , 2020.
  * [20] M. Atzmon and Y. Lipman, “SAL: Sign agnostic learning of shapes from raw data,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * [21] J. Chibane, T. Alldieck, and G. Pons-Moll, “Implicit functions in feature space for 3D shape reconstruction and completion,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020.
  * [22] M. Atzmon and Y. Lipman, “SALD: sign agnostic learning with derivatives,” in _International Conference on Learning Representations (ICLR)_ , 2021.
  * [23] W. Zhao, J. Lei, Y. Wen, J. Zhang, and K. Jia, “Sign-agnostic implicit learning of surface self-similarities for shape modeling and reconstruction from raw point clouds,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2021.
  * [24] J. Chibane, A. Mir, and G. Pons-Moll, “Neural unsigned distance fields for implicit function learning,” in _Conference on Neural Information Processing Systems (NeurIPS)_ , 2020.
  * [25] S. Peng, M. Niemeyer, L. Mescheder, M. Pollefeys, and A. Geiger, “Convolutional occupancy networks,” in _European Conference on Computer Vision (ECCV)_ , 2020.
  * [26] J. Tang, J. Lei, D. Xu, F. Ma, K. Jia, and L. Zhang, “SA-ConvONet: Sign-agnostic optimization of convolutional occupancy networks,” in _International Conference on Computer Vision (ICCV)_ , 2021.
  * [27] B. Ummenhofer and V. Koltun, “Adaptive surface reconstruction with multiscale convolutional kernels,” in _International Conference on Computer Vision (ICCV)_ , 2021.
  * [28] A. Boulch and R. Marlet, “POCO: Point convolution for surface reconstruction,” _arXiv preprint arXiv:2201.01831_ , 2022.
  * [29] P. Erler, S. Ohrhallinger, N. Mitra, and M. Wimmer, “Points2Surf: Learning implicit surfaces from point clouds,” in _European Conference on Computer Vision (ECCV)_ , 2020.
  * [30] N. Sharp and M. Ovsjanikov, “”pointtrinet: Learned triangulation of 3d point sets”,” in _European Conference on Computer Vision (ECCV)_ , 2020.
  * [31] M. Kazhdan and H. Hoppe, “Screened Poisson surface reconstruction,” _ACM Transaction on Graphics._ , 2013.
  * [32] F. Williams, M. Trager, J. Bruna, and D. Zorin, “Neural splines: Fitting 3d surfaces with infinitely-wide neural networks,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2021.
  * [33] S. König and S. Gumhold, “Consistent propagation of normal orientations in point clouds.” in _International Workshop on Vision, Modeling, and Visualization (VMV)_ , 2009.
  * [34] G. Metzer, R. Hanocka, D. Zorin, R. Giryes, D. Panozzo, and D. Cohen-Or, “Orienting point clouds with dipole propagation,” _ACM Transactions on Graphics (TOG)_ , 2021.
  * [35] D. Paschalidou, O. Ulusoy, C. Schmitt, L. Van Gool, and A. Geiger, “Raynet: Learning volumetric 3d reconstruction with ray potentials,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2018.
  * [36] B. Mildenhall, P. P. Srinivasan, M. Tancik, J. T. Barron, R. Ramamoorthi, and R. Ng, “NeRF: Representing scenes as neural radiance fields for view synthesis,” in _European Conference on Computer Vision (ECCV)_ , 2020.
  * [37] M. Oechsle, S. Peng, and A. Geiger, “Unisurf: Unifying neural implicit surfaces and radiance fields for multi-view reconstruction,” in _International Conference on Computer Vision (ICCV)_ , 2021.
  * [38] S. Song, Z. Cui, and R. Qin, “Vis2mesh: Efficient mesh reconstruction from unstructured point clouds of large scenes with learned virtual view visibility,” in _International Conference on Computer Vision (ICCV)_ , 2021\.
  * [39] S. Peng, C. M. Jiang, Y. Liao, M. Niemeyer, M. Pollefeys, and A. Geiger, “Shape as points: A differentiable poisson solver,” in _Conference on Neural Information Processing Systems (NeurIPS)_ , 2021.
  * [40] F. Cazals and M. Pouget, “Estimating differential quantities using polynomial fitting of osculating jets,” _Computer Aided Geometric Design_ , 2005.
  * [41] A. Boulch, “Convpoint: Continuous convolutions for point cloud processing,” _Computers & Graphics_, 2020.
  * [42] W. E. Lorensen and H. E. Cline, “Marching cubes: A high resolution 3D surface construction algorithm,” _ACM SIGGRAPH Computer Graphics_ , vol. 21, no. 4, 1987.
  * [43] Z. Wu, S. Song, A. Khosla, F. Yu, L. Zhang, X. Tang, and J. Xiao, “3d shapenets: A deep representation for volumetric shapes,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015.
  * [44] A. X. Chang, T. Funkhouser, L. Guibas, P. Hanrahan, Q. Huang, Z. Li, S. Savarese, M. Savva, S. Song, H. Su, J. Xiao, L. Yi, and F. Yu, “ShapeNet: an information-rich 3D model repository,” Stanford University, Princeton University, Toyota Technological Institute at Chicago, Tech. Rep., 2015.
  * [45] C. B. Choy, D. Xu, J. Gwak, K. Chen, and S. Savarese, “3D-R2N2: A unified approach for single and multi-view 3D object reconstruction,” in _European Conference on Computer Vision (ECCV)_ , 2016.
  * [46] A. Handa, V. Patraucean, S. Stent, and R. Cipolla, “SceneNet: An annotated model generator for indoor scene understanding,” in _International Conference on Robotics and Automation (ICRA)_ , 2016.
  * [47] A. Dai, A. Chang, M. Savva, M. Halber, T. Funkhouser, and M. Nießner, “ScanNet: Richly-annotated 3D reconstructions of indoor scenes,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017.
  * [48] A. Knapitsch, J. Park, Q.-Y. Zhou, and V. Koltun, “Tanks and Temples,” _ACM Transactions on Graphics (TOG)_ , 2017.
  * [49] S. Seitz, B. Curless, J. Diebel, D. Scharstein, and R. Szeliski, “A comparison and evaluation of multi-view stereo reconstruction algorithms,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2006.
  * [50] D. Cernea, “OpenMVS: Multi-view stereo reconstruction library,” 2020. [Online]. Available: https://cdcseacave.github.io/openMVS
  * [51] J. Huang, Y. Zhou, and L. Guibas, “Manifoldplus: A robust and scalable watertight manifold surface generation method for triangle soups,” _arXiv preprint arXiv:2005.11621_ , 2020.
  * [52] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017.
  * [53] R. Jensen, A. Dahl, G. Vogiatzis, E. Tola, and H. Aanæs, “Large scale multi-view stereopsis evaluation,” in _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2014.

## Supplementary Material

In this supplementary document, we first provide additional information about
the datasets that we use (Section VII), formal definitions of the evaluation
metrics (Section VIII) and additional quantitative and qualitative results
(Section IX). Our code, data and pretrained models can be found online:
https://github.com/raphaelsulzer/dsrv-data .

##  VII Datasets

###  VII-A Scanning procedure

In Figure 7, we represent a visualization of our scanning procedure.
Viewpoints and ray target points are uniformly sampled on the surface of the
spheres.

Figure 7: Scanning Procedure. We randomly place sensors on two spheres (red
and blue) around the object, and consider rays aiming at uniformly sampled
points on a sphere inscribed in the convex hull of the object (green).

###  VII-B ModelNet

We use the official ModelNet10 dataset and make the models watertight using
ManifoldPlus [51]. We scan the watertight models using the procedure described
above.

###  VII-C ShapeNet

We use the watertight models provided111https://s3.eu-
central-1.amazonaws.com/avg-projects/occupancy_networks/data/watertight.zip by
the authors of Occupancy Networks [2] and scan the models using the procedure
described above. We apply a transformation to the models (and scans) to match
their orientation to the orientation of the ModelNet10 objects (except for
networks marked with $\dagger$ in Table IV, which were trained with the
original orientation).

###  VII-D Synthetic Rooms Dataset

We use the watertight scenes provided222https://s3.eu-
central-1.amazonaws.com/avg-
projects/convolutional_occupancy_networks/data/room_watertight_mesh.zip by the
authors of ConvONet [25]. We scan the scenes using the procedure described
above, limiting the sensors to lie in the upper hemisphere only.

###  VII-E Tanks And Temples, Middlebury and DTU

For the reconstructions in Figure 5, 6 and Figure 8, we downsample the input
point clouds to 10,000 points.

Bare |  |  |  |  |  |
---|---|---|---|---|---|---
Augmented |  |  |  |  |
| Input |  ConvONet-3D [25] |  Points2Surf [29] |  Shape As Points [39] |  POCO [28] | HD Scan
Figure 8:  Out-of-Domain Object-Level Reconstruction. Reconstructed shape from
a MVS point cloud of _scan1_ from DTU, using four different DSR methods
trained on ModelNet10. The top row uses the bare point cloud as input, and the
bottom row uses the point cloud augmented with visibility information. _HD
Scan_ is a high-density point cloud.

##  VIII Metrics

We evaluate the quality of reconstructions with the volumetric IoU (IoU),
symmetric Chamfer distance (CD) and normal consistency (NC).

Let $\mathcal{M_{G}}$ be the ground truth mesh and $\mathcal{M_{P}}$ be the
reconstructed mesh. The volumetric IoU is defined as:

$\displaystyle\text{IoU}\left(\mathcal{M_{G}},\mathcal{M_{P}}\right)=$
$\displaystyle\frac{\operatorname{volume}(\mathcal{M_{G}}\cap\mathcal{M_{P}})}{\operatorname{volume}(\mathcal{M_{G}}\cup\mathcal{M_{P}})},$

We approximate volumetric IoU by sampling $100,000$ points in the union of the
bounding boxes of $\mathcal{M_{G}}$ and $\mathcal{M_{P}}$.

To compute the Chamfer distance and normal consistency, we sample a set of
points $S_{G}$ on the ground-truth mesh and a set of points $S_{P}$ on the
reconstructed mesh with $|S_{G}|=|S_{P}|=100,000$. We approximate the two-
sided Chamfer distance between $\mathcal{M_{G}}$ and $\mathcal{M_{P}}$ as
follows:

$\displaystyle\text{CD}(\mathcal{M_{G}},\mathcal{M_{P}})=$
$\displaystyle\frac{1}{2|S_{{G}}|}\sum_{x\in S_{G}}\min_{y\in
S_{P}}{||x-y||}_{2}$ $\displaystyle+$
$\displaystyle\frac{1}{2|S_{{P}}|}\sum_{y\in S_{P}}\min_{x\in
S_{G}}{||y-x||}_{2}$

Let $n(x)$ be the unit normal associated to a point $x$ taken on a mesh, and
$\langle\cdot{,}\cdot\rangle$ the Euclidean scalar product in
$\mathbb{R}^{3}$. Normal consistency is defined as:

$\displaystyle\text{NC}(\mathcal{M_{G}},\mathcal{M_{P}})=$
$\displaystyle\frac{1}{2|S_{{G}}|}\sum_{x\in S_{G}}\left\langle
n(x),n\left(\operatorname*{argmin}_{y\in
S_{P}}{||x-y||}_{2}\right)\right\rangle$ $\displaystyle+$
$\displaystyle\frac{1}{2|S_{{P}}|}\sum_{y\in S_{P}}\left\langle
n(y),n\left(\operatorname*{argmin}_{x\in
S_{G}}{||y-x||}_{2}\right)\right\rangle$

##  IX Additional Results

###  IX-A Runtimes

In Table V, we report detailed runtimes for the tested methods, with and
without visibility information.

Adding sightline vectors does not significantly increase the runtime for any
of the tested methods. The effect of auxiliary points depends on the method.
For ConvONet, most of the processing time is spent computing grid features.
The encoding of 3d points is performed by a small PointNet network [52], whose
runtime is only a small fraction of the total time. As a consequence, adding
APs does not incur significant changes in computation time. In contrast,
Points2Surf uses a large point encoding network and is $2.2$ times slower with
APs. Shape As Points is $1.3$ times slower due to the fact that we decode $3$
times as many points as the baseline method. POCO is also essentially
unaffected by the addition of auxiliary points as they only impact the first
(small) layer of the point-convolution backbone.

TABLE V:  Runtimes for Object-Level Reconstruction.

Average times (in seconds) for reconstructing one object from a point cloud of
3 000 points with and without sightline vectors (SV) or auxiliary points (AP).
MC is marching cubes. Times are averaged over the ModelNet10 test set.

Model | SV | AP | Encoding | Decoding | MC | Total
---|---|---|---|---|---|---
ConvONet-2D [25] |  |  | 0.016 | 0.25 | 0.17 | 0.44
ConvONet-2D [25] | ✓ |  | 0.016 | 0.27 | 0.17 | 0.47
ConvONet-2D [25] | ✓ | ✓ | 0.016 | 0.26 | 0.17 | 0.45
Points2Surf [29] |  |  | 69.06 | 11.51 | 80.57
Points2Surf [29] | ✓ |  | 71.92 | 11.35 | 83.27
Points2Surf [29] | ✓ | ✓ | 173.2 | 11.41 | 184.7
Shape As Points [39] |  |  | 0.022 | 0.017 | 0.047 | 0.088
Shape As Points [39] | ✓ |  | 0.023 | 0.017 | 0.046 | 0.086
Shape As Points [39] | ✓ | ✓ | 0.024 | 0.041 | 0.047 | 0.114
POCO [28] |  |  | 0.088 | 13.72 | 0.33 | 15.74
POCO [28] | ✓ |  | 0.091 | 13.68 | 0.33 | 15.66
POCO [28] | ✓ | ✓ | 0.093 | 13.70 | 0.33 | 15.67

###  IX-B DTU Dataset

In Figure 8, we show the reconstruction of an MVS point cloud, generated with
OpenMVS, of _scan1_ from the DTU dataset [53]. The point cloud represents an
open scene, while all methods were trained on the closed ModelNet10 objects.
The methods using our augmented point clouds with visibility information cope
much better with this domain shift.

###  IX-C ShapeNet

Bare |  |  |  |  |  |
---|---|---|---|---|---|---
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
| Input |  ConvONet-3D [25] |  Points2Surf [29] |  Shape As Points [39] |  POCO [28] | Ground Truth
Figure 9:  Object-Level Reconstruction on ShapeNet. Reconstructed shapes from
the ShapeNet test set using four different DSR methods trained on ModelNet10.
Top rows of each object use the bare point cloud as input, and bottom rows use
the point cloud augmented with visibility information. The last two rows show
a cut of the reconstructions that are shown on the two other rows immediately
above.

We show the results of object-level reconstruction on ShapeNet in Figure 9
from methods trained on ModelNet10. All methods benefit from added visibility
information. In particular, ConvONet produces very accurate and complete
surfaces of the unseen shape classes.

A common problem for most baseline methods is the false reconstruction of
hollow shapes with enclosed outside space. Using visibility information can
address this issue and the artifacts then do not occur.

###  IX-D ModelNet10

Bare |  |  |  |  |  |
---|---|---|---|---|---|---
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
Bare |  |  |  |  |  |
Augmented |  |  |  |  |
| Input |  ConvONet-2D [25] |  Points2Surf [29] |  Shape As Points [39] |  POCO [28] | Ground Truth
Figure 10:  Object-Level Reconstruction On ModelNet10. Reconstructed shapes
from the ModelNet10 test set using four different DSR methods trained on
ModelNet10. Top rows of each object use the bare point cloud as input, and
bottom rows use the point cloud augmented with visibility information.

In Figure 10, we represent additional results of object-level reconstruction
on ModelNet10. Concave parts of the objects are frequently reconstructed more
accurately by methods using point clouds with visibility information. Surfaces
also tend to be more complete when visibility information is used.

###  IX-E Synthetic Rooms Dataset

Bare |  |  |  |  |
---|---|---|---|---|---
Augmented |  |  |  |  |
ConvONet [25] |  |  |  |  |
ConvONet [25] (+SV+AP) |  |  |  |  |
Ground Truth |  |  |  |  |
Figure 11:  Scene-Level Reconstruction on Synthetic Rooms. Reconstructed
scenes of the synthetic room dataset using ConvONet [25] in sliding-window
mode, with and without visibility information.

In Figure 11, we show the reconstruction results of ConvONet on Synthetic
Rooms, with and without visibility information. In contrast to [25], we
evaluate in sliding window mode, which explains the difference with figures in
[25]. Improvements here are visually not as obvious as with other datasets.


# A More Scalable Sparse Dynamic Data Exchange

1st Andrew Geyko
Department of Computer Science
University of New Mexico
Albuquerque, USA
<EMAIL_ADDRESS>2nd Gerald Collom
Department of Computer Science
University of New Mexico
Albuquerque, USA
<EMAIL_ADDRESS>3rd Derek Schafer
Department of Computer Science
University of New Mexico
Albuquerque, USA
<EMAIL_ADDRESS>4th Patrick Bridges
Department of Computer Science
University of New Mexico
Albuquerque, USA
<EMAIL_ADDRESS>5th Amanda Bienz
Department of Computer Science
University of New Mexico
Albuquerque, USA
<EMAIL_ADDRESS>
###### Abstract

Parallel architectures are continually increasing in performance and scale,
while underlying algorithmic infrastructure often fail to take full advantage
of available compute power. Within the context of MPI, irregular communication
patterns create bottlenecks in parallel applications. One common bottleneck is
the sparse dynamic data exchange, often required when forming communication
patterns within applications. There are a large variety of approaches for
these dynamic exchanges, with optimizations implemented directly in parallel
applications. This paper proposes a novel API within an MPI extension library,
allowing for applications to utilize the variety of provided optimizations for
sparse dynamic data exchange methods. Further, the paper presents novel
locality-aware sparse dynamic data exchange algorithms. Finally, performance
results show significant speedups up to 20x with the novel locality-aware
algorithms.

###### Index Terms:

sparse dynamic data exchange, MPI, HPC, parallel performance, locality-aware,
extension library

##  I Introduction

Parallel architectures are continuously improving in scale, computational
power, and efficiency, allowing for increasingly large and efficient high
performance computing (HPC) applications. Despite improvements in hardware
capabilities, parallel applications fail to fully utilize emerging hardware
due to high costs associated with data movement and inter-process
communication. Standard communication algorithms treat all communication as
equal. However, the costs of various paths of communication on current and
emerging computers vary greatly with relative locations of the sending and
receiving processes, with notable differences between intra-socket, inter-
socket, and inter-node communication.

A large number of parallel applications, including sparse solvers and
simulations, are bottlenecked by the performance of irregular communication.
This type of communication consists of each process sending messages to a
subset of other processes, with the pattern of communication dependent on the
specific problem being solved such as the sparsity pattern of a given matrix.
One bottleneck in irregular communication is forming a communication pattern,
determining the set of processes to which each process must send as well as
the set from which it receives. Often, such as in the case of row-wise
partitioned sparse matrix operations, each processes can use the local
sparsity pattern to determine which data should be received from other
processes, but has insufficient information regarding processes to which data
should be sent. As a result, a crucial first step in every communication phase
is determining this pattern of communication. The creation of this
communication pattern requires a _sparse dynamic data exchange_ phase
(_SDDE_). The cost of creating the communication pattern through the SDDE
often outweighs the cost of subsequent communication unless a number of
iterations of communication are performed.

There are two main use cases for sparse dynamic data exchanges: exchanges with
constant data sizes, such as determining only the communication pattern and
per-message byte count, and variable-sized data exchanges, in which each
process must dynamically send a variable amount of data to a subset of other
processes. Cell adaptive mesh refinement applications, such as CELLAR. require
an SDDE with constant size messages each time the mesh is altered. Each
process can locally determine the subset of processes to which it must send
and what data should be sent to each. However, an SDDE is needed to determine
from which processes each receives, along with the size of each receive. This
SDDE requires constant message sizes as only message sizes must be exchanged.
A large majority of SDDEs, however, require variable-sized exchanges. Linear
solvers typically consist of distributing sparse matrices across processes
row-wise, allowing each process to locally determine the receive portion of
the communication pattern. However, in this case, an SDDE is needed for each
process to find the subset of processes to which it must send along with which
data to send to each. As a result, indices corresponding to these data values
must be sent, requiring each message to be equal to the size of subsequent
iterations of communication, and therefore variable.

There are many existing approaches for the sparse dynamic data exchange, each
hand implemented within applications that rely on the SDDE. There are two
typical implementations for the sparse dynamic data exchange, which can apply
to both variable and constant data sizes. The personalized method consists of
all processes performing an MPI_Allreduce to determine data exchange sizes
before indices are dynamically communicated. Alternatively, the non-blocking
implementation proposed by Torsten et al. dynamically exchanges indices
without the need for an initial reduction. The methods each outperform the
other in a subset of sparse dynamic data exchange problems, with performance
dependent on communication patterns, data exchange sizes, and process count.
The personalized method outperforms the non-blocking approach in the case of
SDDEs with very large message counts, as the allreduce has minimal overhead.
However, SDDEs among a large set of process in which there are few messages
are more optimally implemented through the non-blocking algorithm, avoiding
the significant overhead of the reduction.

Applications of constant-sided SDDEs, such as Cellar, have further optimized
the algorithm with one-sided MPI remote-memory access (RMA) communication. As
all messages are of size count, each process $p$ can use MPI_Put to directly
put its contribution into the destination process’s receive buffer at position
$p\times\texttt{count}$. This implementation, however, cannot be naturally
extended to the variable-sized SDDE as variable sized data would need to be
exchanged after the message sizes are exchanged with RMA. While the initial
RMA step would allow for a standard neighbor collective or point-to-point
exchange thereafter, rather than dynamic communication involving probes, very
large scales are required for RMA to outperform the personalized and non-
blocking approaches. At these scales, the significant queue search costs, or
overhead of matching received data with posted MPI_Irecvs, cause standard
point-to-point exchanges to be more costly than dynamically receiving the
first available method with MPI_Probe. As a result, the RMA SDDE presented in
this paper is only applied to SDDEs with constant-sized messages.

This paper presents optimizations to existing sparse dynamic data exchanges
through message aggregation. Locality-aware optimizations greatly improve the
performance and scalability of irregular data exchanges within the context of
iterative methods. Communication within a region, such as intra-socket,
greatly outperforms inter-region messages, such as inter-socket or inter-node.
Locality-aware aggregation techniques gather data within a region and minimize
the number of messages communicated between regions. This paper presents two
novel locality-aware sparse dynamic data exchange algorithms, which optimize
point-to-point communication constraints for large message counts through
message aggregation. Locality-aware variants are presented for both the
personalized and the non-blocking SDDE approaches, _achieving up to 20x
speedup over existing personalized and non-blocking SDDE methods._ While the
paper does not investigate locality-aware aggregation within the RMA SDDE, the
presented aggregation techniques could be naturally extended to one-sided
communication.

All SDDE algorithms presented in this paper are implemented within the open-
source MPI eXtension library MPI Advance [1]. Application programmers can
replace current hand-implemented SDDE methods within widely used solvers with
the standard API presented within this paper to utilize the presented
optimizations. The authors present scaling studies of the various algorithms
as implemented within MPI Advance and analyze the trade-offs between the
existing and locality-aware approaches.

The remainder of the paper is outlined as follows. An overview of the sparse
dynamic data exchange problem and related works are presented in Section II. A
common MPI eXtension (MPIX) API for SDDE algorithms is detailed in Section 4.
Existing algorithms for the SDDE are detailed in Section IV, while novel
locality-aware extensions are described in Section IV-D. Scaling studies are
presented in Section V, along with analysis of when each algorithm performs
best. Finally, concluding remarks and future work are discussed in Section VI.

##  II Background

The sparse dynamic data exchange problem consists of a set of processes, each
of which contains local data and is looking to perform an operation on the
global collection of data. Each process must receive a subset of the global
data before it can complete its local portion of the operation. Initially,
every process knows which global data it must receive from other processes.
For instance, consider a sparse matrix operation, in which the sparse matrix
and corresponding vectors are partitioned row-wise across all available
processes. Each process must receive data corresponding to each non-zero
column within its subset of rows. However, before data can be exchanged, each
process must discover the subset of processes to which it must send all or
part of its local data. The problem is formalized as follows.

######  Definition 1 (Sparse Dynamic Data Exchange Problem)

Let $\Sigma:=\\{\sigma_{1},\dots,\sigma_{n}\\}$ be a collection of processes.
For each $\sigma_{i}\in\Sigma$, let there be an associated local data segment
$D_{i}$ and a subset $A_{i}$ of $\Sigma$ of processes from which $\sigma_{i}$
needs to receive data from. The solution to the problem is a communication
pattern such that after execution, each $\sigma_{i}\in\Sigma$ holds a copy of
the set $\\{j\;|\;\sigma_{i}\in A_{j}\\}$.

Often, parallel applications rely on the SDDE problem when forming a
communication pattern consisting of all processes to which each rank sends and
receives data, the sizes of these messages, and indices of the data to be
communicated. There are multiple existing approaches to forming these
communication patterns, including the personalized, non-blocking, and RMA
approaches, detailed in Section IV. The performance of existing SDDE
approaches consists of overhead costs, such as window synchronization or
reduction costs, as well as the cost of point-to-point or one-sided dynamic
communication of data. The trade-offs between existing implementations depend
on whether the bottleneck of an SDDE is due to overhead or transfer of data.
Factors that determine this overhead include the number of messages each
process communicates, the size of these messages, and the total number of
processes. Furthermore, the performance of the required point-to-point or one-
sided communication is dependent on the locality of the messages, MPI
implementation, and system architecture.

###  II-A Parallel Matrix Partitioning

Figure 1: An example of a $4\times 4$ matrix and a vector distributed across 4
processes.

Sparse matrices are typically partitioned row-wise, so that each process holds
a contiguous subset of rows of the matrix along with corresponding vector
values, as shown in Figure 1. Consider a system with $n$ rows split row-wise
across $p$ processes. Each process holds $\frac{n}{p}$ consecutive rows of the
matrix assuming the number of processes divides the number of rows; otherwise,
some processes will hold one additional row. Furthermore each process holds
the vector values corresponding to its local columns. There are other
partitioning strategies, such as column-wise or 2.5d partitions [2, 3]. All
SDDE strategies discussed in this paper are applied to row-wise partitions of
sparse matrices and extend naturally to column-wise partitioning. While multi-
dimensional partitioning strategies are outside the scope of this paper, the
required SDDE problems are similar regardless of the partitioning strategy.

Send Variable | Input/Output | Definition | Receive Variable | Input/Output | Definition
---|---|---|---|---|---
send_nnz | Input | Number of dynamic sends | recv_nnz | Input/Output | Number of dynamic sends
send_size | Input | Size of dynamic sends | recv_size | Input/Output | Size of dynamic sends
dest | Input | Destinations of messages | src | Output | Destinations of messages
sendcounts | Input | Send per-message count | recvcounts | Output | Receive per-message count
sdispls | Input | Send displacements | rdispls | Input | Receive displacements
sendtype | Input | Datatype of sends | recvtype | Output | Datatype of receives
sendvals | Input | Data to send | recvvals | Output | Data to receive
TABLE I: Variables required in MPIX sparse dynamic data exchange APIs. Red
text indicates variables only required in the variable SDDE, which black
variables are required in both APIs.

###  II-B Data Dependencies in Matrix-Vector Multiplications

Consider the matrix displayed in Figure 1. When performing a matrix-vector
product $A\vec{v}$ in parallel, each process must retrieve vector entries
corresponding to off-process non-zeroes to perform the operation. For example,
$P_{2}$ needs vector entries $v_{0}\text{ and }v_{3}$ in order to performs its
local portion of the matrix-vector product. A complete data dependency chart
is given by the following figure:

Figure 2: Data dependency chart when performing the matrix vector product
illustrated in Figure 1

As shown in this example, a process does not know in advance which other
processes need portions of its data as this communication pattern is dependent
on the associated sparsity pattern. Therefore, before data can be exchanged,
SDDE is required to determine the subset of processes with which each must
communicate, along with which data must be communicated. Communication
patterns for irregular data exchanges are typically sparse, with each process
communicating with a small subset of the total number of processes, typically
$\mathcal{O}(\log(n))$ neighbors. For more dense communication patterns, all-
to-all exchanges are more efficient, removing the needs for the SDDE phase and
allowing for optimizations such as pairwise exchange.

###  II-C Related Works

Irregular communication bottlenecks many parallel applications, including
sparse solvers and simulations. As a result, there has been a significant
amount of research into irregular communication optimizations. However, the
majority of research into irregular communication optimizes messages after the
communication pattern has been formed. Iterative point-to-point communication
has been optimized within MPI through persistent communication [4, 5], which
performs initialization once before reusing the communicator, amortizing setup
costs over all iterations. MPI 4.0 provides partition communication [6, 7], an
optimization for multithreaded applications that allows each thread to send a
portion of a message when ready.

Many architecture-aware optimizations for irregular and collective
communication utilize locality to reduce communication costs. Hierarchical
communication reduces inter-node communication costs by utilizing only one or
a small number of MPI processes per node during inter-node communication [8,
9]. Locality-aware aggregation minimizes the number and size of non-local
messages, such as inter-node, through local, such as intra-node, aggregation
of data [10, 11, 12]. Similar node-aware aggregation strategies, such as YGM,
have been applied to graph operations [13]. Multilane algorithms minimize the
size of inter-node messages, having each process per node send an equal and
minimal amount of data [14].

Communication constraints can be minimized through optimal partitions. Rather
than a standard row-wise partition, graph or hypergraph partitioners minimize
edge cuts between partitions [15]. Each partition is then assigned to a single
process or node, reducing the amount of data needed to be exchanged throughout
applications. While graph partitioners can be used in combination with sparse
dynamic data exchanges, the use of these methods is outside the scope of this
paper due to their large overheads.

##  III MPI Extension APIs

A large variety of sparse dynamic data exchange implementations are used
throughout existing parallel applications, with each application required to
determine the optimal method for their program. This section presents the API
that the authors have added to an open-source MPI Extension library, MPI
Advance. Existing methods and novel optimizations are added within the
provided APIs, allowing an existing application to access any of the available
optimizations rather than hand optimize. Further, the API allows for future
research in which the optimal algorithm can be selected dynamically. Note that
sparse data exchanges exist within MPI through neighbor collectives, but these
methods do not allow for sparse _dynamic_ data exchanges.

Sparse dynamic data exchange algorithms can be represented with two separate
APIs: MPIX_Alltoall_crs and MPIX_Alltoallv_crs. All variables required for the
APIs are defined in Table I.

The MPIX_Alltoall_crs method allows for dynamically exchanging constant data
sizes with a sparse set of neighbors. One use case of this method occurs when
processes know all data that must be sent to other processes, but have no
knowledge of the receive portion of the communication pattern. This type of
SDDE algorithm occurs within cell-based adaptive mesh refinement solvers, such
as CELLAR [16].

The MPIX_Alltoall_crs API is defined in Figure 3.

⬇

int MPIX_Alltoall_crs(int send_nnz,

int* dest, int sendcount,

MPI_Datatype sendtype,

void* sendvals, int* recv_nnz,

int* src, int recvcount,

MPI_Datatype recvtype,

void* recvvals, MPIX_Info* xinfo,

MPIX_Comm* xcomm)

Figure 3: MPIX_Alltoall_crs API

This method takes the send portion of the SDDE as input and returns received
data. The variable recv_nnz is both input and output, allowing the user to
input this information if it is already known rather than requiring the method
to redetermine number of receives. Note the API is following current MPI
standard requirements, and as a result recvvals must be allocated, potentially
to some upper-bound, before the method is instantiated.

The MPIX_Alltoallv_crs method provides an interface for dynamically exchanging
variable sized data with a sparse set of neighbors. This method is required
for any distributed sparse matrix operation, in which processes know all data
that must be received, but have no knowledge of processes to which they must
send or what data must be send to each. As a result, a variable sized message
must be exchanged containing this boundary information. Typically, this
variable SDDE is used to initially form a communication package to then be
used during each future sparse matrix operation, during which communication is
typically implemented either through point-to-point messages or neighborhood
collectives. Variable SDDEs exist within many widely used solvers, including
Hypre BoomerAMG.

The MPIX_Alltoallv_crs API is defined in Figure 4.

⬇

int MPIX_Alltoallv_crs(int send_nnz,

int send_size, int* dest,

int* sendcounts, int* sdispls,

MPI_Datatype sendtype,

void* sendvals, int* recv_nnz,

int* recv_size, int* src,

int* recvcounts, int* rdispls,

MPI_Datatype recvtype,

void* recvvals, MPIX_Info* xinfo,

MPIX_Comm* comm)

Figure 4: MPIX_Alltoallv_crs API

Similar to the non-variable API, this method takes the send portion of the
variable SDDE as input and returns received data. The parameters recv_nnz and
recv_size are both input and output, allowing the user to input both receive
count and/or receive size if either are known. Due to MPI standard
requirements, this method also requires all pointers to be returned by the
method to be allocated before they are passed to the method. This includes
recvcounts and rdispls, which each need to be allocated to hold recv_nnz
integers at mininum, as well as recvvals, which needs to be allocated to hold
at least recv_size variables of size recvtype.

##  IV Implementations

There are a number of existing methods for determining the communication
pattern of a given irregular data exchange. The personalized and non-blocking
methods, described in detail below, are widely used throughout existing
parallel applications. Furthermore, an RMA-based optimization, as implemented
within Cellar, allows for MPIX_Alltoall_crs communication to be performed with
one-sided communication. All algorithms presented in this section, except for
the RMA method, can be applied to both the MPIX_Alltoall_crs and
MPIX_Alltoallv_crs methods, with minor changes required for the variable-sized
SDDE, as highlighted in red.

###  IV-A Personalized Method

The standard personalized method, utilized in a number of irregular codebases
and described in [17], is detailed in Algorithm 1.

Input: rank

{Process ID}

args

{MPIX_Alltoallv_crs Arguments}

$\texttt{ctr}\leftarrow 0$

for  _$i\leftarrow 0$ to send_nnz_ do

proc = dest[$i$]

size = sendcounts[$i$]

sizes[proc] = size

MPI_Isend size at sendvals[ctr] to proc

$\texttt{ctr}\leftarrow\texttt{ctr}+\texttt{size}$

MPI_Allreduce(sizes) for  _$i\leftarrow 0$ to sizes[rank]_ do

Probe for message and receive dynamically

Wait for all sends to complete

Algorithm 1 Personalized

The personalized protocol completes a sparse dynamic data exchange by
performing an MPI_Allreduce to gather either the total number of size of
messages that every process will receive during the dynamic communication
phase. Data to be exchanged is then communicated with non-blocking sends.
Finally, each process dynamically receives all messages sent to it using
MPI_Probe.

When process counts are large, the overhead of the MPI_Allreduce increases.
However, as message count increases, the reduction overhead is quickly
outweighed by the dynamic exchange of data. The use of probes requires
received data to be added to the unexpected message queue, further reducing
performance for large messages. However, the MPI_Allreduce can yield large
benefits as it allows all communication structures to be allocated before the
dynamic communication step.

###  IV-B Non-blocking Method

The non-blocking method, presented in Torsten et al. [17] and detailed in
Algorithm 2, is a modification of the standard algorithm which avoids the
overhead and collective synchronization required by the MPI_Allreduce.

Input: rank

{Process ID}

args

{MPIX_Alltoallv_crs Arguments}

$\texttt{ctr}\leftarrow 0$

for  _$i\leftarrow 0$ to send_nnz_ do

proc = dest[$i$]

size = sendcounts[$i$]

MPI_Isend size at sendvals[ctr] to proc

$\texttt{ctr}\leftarrow\texttt{ctr}+\texttt{size}$

while  _All sends have not completed_

if  _Non-blocking probe finds a message_

Dynamically receive message

Non-blocking barrier

while  _Non-blocking barrier has not completed_

if  _Non-blocking probe finds a message_

Dynamically receive message

Algorithm 2 Non-blocking

In this algorithm, each process $p$ first sends non-blocking synchronous sends
to every process to which it must communicate. A synchronous send is only
considered complete when the destination process posts the associated receive.
Each process dynamically receives messages until the synchronous sends have
been completed on all processes. This is accomplished through MPI_Iprobe,
which checks if a message is available, and only then is data received. While
checking for messages, each process also tests its synchronous sends for
completion Once all message sent from a process have been received, that
process calls a non-blocking barrier to indicate that it is finished sending.
This process continues probing for messages until all processes have reached
the barrier, at which point all processes’ sends are received and dynamic
communication is completed.

While the non-blocking algorithm improves over the personalized method for
large process counts by reducing synchronization and cost associated with the
MPI_Allreduce, communication pattern structures must now be dynamically
allocated. Furthermore, dynamic communication and unexpected messages remain a
dominant cost of the method.

###  IV-C RMA Constant-Size SDDE

One-sided optimizations exist for sparse dynamic data exchanges in which all
data is of a constant size, as implemented in Cellar and detailed in Algorithm
3.

Input: rank

{Process ID}

args

{MPIX_Alltoall_crs Arguments}

Create window

Synchronize on window

$\texttt{proc\\_vals}[\texttt{num\\_procs}\times\texttt{sendcount}]$

$\texttt{ctr}\leftarrow 0$

for  _$i\leftarrow 0$ to send_nnz_ do

$\texttt{proc}\leftarrow\texttt{dest}[i]$

$\texttt{orig\\_address}\leftarrow\texttt{sendvals}[\texttt{ctr}]$

$\texttt{dest\\_address}\leftarrow\texttt{proc\\_vals}[\texttt{rank}]\times\texttt{sendcount}$

MPI_Put sendcount values from orig_address to dest_address

$\texttt{ctr}\leftarrow\texttt{ctr}+\texttt{sendcount}$

Synchronize on window

for  _$i\leftarrow 0$ to num_procs_ do

$\texttt{}{val}\leftarrow\texttt{proc\\_vals}[i\times\texttt{sendcount}]$

if  _val_

$\texttt{recvvals}[\texttt{recv\\_nnz}++]\leftarrow\texttt{val}$

Algorithm 3 RMA (MPIX_Alltoall_crs only)

Each process firsts allocates a shared memory window with enough space for all
processes to add sendcount values. Then, each process can independently use
MPI_Put to move the required data to position
$\texttt{rank}\cdot\texttt{sendcount}$ on each corresponding process. Finally,
each process can move received data from the window into the recvvals array.

The one-sided approach requires window creation, but this can be amortized
over the cost of the application, and can potentially be reused within any
subsequent one-sided communication. There is also a synchronization required
to use a shared memory window. The author’s currently use MPI_Fence for this,
but more optimized forms of synchronization, such as locks, could be explored.
RMA sparse dynamic data exchanges allow for all data to be exchanged without
any dynamic MPI communication. However, the method does not naturally extend
to variable-sized dynamic data exchanges. It can be used to determine
communication sizes, followed a standard data exchange to avoid dynamic
messages, but the authors were unable to find a case where this would
outperform other methods due to large matching costs. It is possible, however,
locality-aware data exchanges could improve the performance of the subsequent
data exchange, but this optimization is outside the scope of this paper.

###  IV-D Locality-Aware Method

At large scales, the cost of both the personalized and non-blocking methods
are bottlenecked by the point-to-point communication requirements. Irregular
communication has been extensively modeled and analyzed, showing large costs
associated with high inter-socket message counts due to lower bandwidth than
intra-socket messages, injection bandwidth limits, queue search (or matching)
costs, and network contention. This irregular communication can be optimized
through locality-awareness, aggregating messages within a region, such as a
socket, to minimize the number of messages communicated between regions. This
greatly reduces the number of times inter-region latency is incurred along
with associated queue search costs. While other locality-aware techniques
remove duplicate values to reduce inter-region message sizes, this paper only
concatenates messages as the overhead of determining duplicate values would
outweigh benefits of reduced message size within a single iteration of
communication, such as required during the SDDE problem.

(a) dielFilterV2clx

(b) germany_osm

(c) human_gene1

(d) NLR

Figure 5: The cost of the various MPIX_alltoall_crs methods using Mvapich2.

(a) dielFilterV2clx

(b) germany_osm

(c) human_gene1

(d) NLR

Figure 6: The cost of the various MPIX_alltoall_crs methods using OpenMPI.

(a) dielFilterV2clx

(b) germany_osm

(c) human_gene1

(d) NLR

Figure 7: The cost of the various MPIX_alltoallv_crs methods using Mvapich2.

(a) dielFilterV2clx

(b) germany_osm

(c) human_gene1

(d) NLR

Figure 8: The cost of the various MPIX_alltoallv_crs methods using OpenMPI.

This paper presents locality-aware aggregation for both the personalized and
non-blocking sparse dynamic data exchanges, as shown in Algorithms 4 and 5,
respectively.

Input: rank

{Process ID}

args

{MPIX_Alltoallv_crs Arguments}

local_rank

{Local rank of process within region}

region_size

{Number of processes per region}

$\texttt{ctr}\leftarrow 0$

for  _$i\leftarrow 0$ to send_nnz_ do

proc = dest[$i$]

region = get_region(proc)

size = sendcounts[$i$]

Copy proc, size, and sendvals[ctr:ctr+size] into buf[region]

$\texttt{ctr}\leftarrow\texttt{ctr}+\texttt{size}$

for each  _region to which rank sends_

$\texttt{proc}=\texttt{region}\cdot\texttt{region\\_size}+\texttt{local\\_rank}$

$\texttt{sizes}[proc]=\texttt{size}$

Non-blocking send of buf[region] to proc

MPI_Allreduce(sizes)

for  _$i\leftarrow 0$ to sizes[rank]_ do

Probe for message and receive dynamically

for  _$\texttt{proc}\leftarrow 0$ to region_size_ do

sizes[proc] = size of buf[proc]

Non-blocking send of buf[proc] to proc

MPI_Allreduce(sizes)

for  _$i\leftarrow 0$ to sizes[rank]_ do

Probe for message and receive dynamically

Wait for all sends to complete

Algorithm 4 Locality-Aware Personalized

Input: rank

{Process ID}

args

{MPIX_Alltoallv_crs Arguments}

local_rank

{Local rank of process within region}

region_size

{Number of processes per region}

for  _$i\leftarrow 0$ to send_nnz_ do

proc = dest[$i$]

region = get_region(proc)

size = sendcounts[i]

Copy proc, size, and sendvals[ctr:ctr+size] into buf[region]

$\texttt{ctr}\leftarrow\texttt{ctr}+\texttt{size}$

for each  _region to which rank sends_

$\texttt{proc}=\texttt{region}\cdot\texttt{region\\_size}+\texttt{local\\_rank}$
Non-blocking send of buf[region] to proc

while  _All sends have not completed_

if  _Non-blocking probe finds a message_

Dynamically receive buf from origin_proc

for each  _proc , size, indices in buf_

Copy origin_proc, size, indices into local_buf[proc]

Non-blocking barrier

while  _Non-blocking barrier has not completed_

if  _Non-blocking probe finds a message_

Dynamically receive buf from origin_proc

for each  _proc , size, indices in buf_

Copy origin_proc, size, indices into local_buf[proc]

for  _$\texttt{proc}\leftarrow 0$ to region_size_ do

sizes[proc] = size of buf[proc]

Non-blocking send of buf[proc] to proc

MPI_Allreduce(sizes)

for  _$i\leftarrow 0$ to sizes[rank]_ do

Probe for message and receive dynamically

Wait for all sends to complete

Algorithm 5 Locality-Aware Nonblocking

Both locality-aware optimizations consist of concatenating all messages that a
process $p$ must send to any process within a single region. For example, the
node-aware version of this method would aggregate all messages that process
$p$ sends to any process on node $n$, and send all data as a single message to
a corresponding process. The corresponding process is determined to be the
process in the region of destination with the same local rank as the sending
process $p$, where local rank is determined as the rank of $p$ within its
region. For instance, if there are PPN processes per region and ranks are laid
out sequentially across the regions, each process $p$ has local rank
$\frac{p}{\texttt{PPN}}$. A process with local rank $r$ communicates only with
other processes of local rank $r$ during the inter-region step.

After the inter-region communication completes with either the personalized or
non-blocking approach, all data is redistributed within the region. This is
currently implemented as a personalized algorithm due to the fact that there
are often a small number of processes within each socket or node, and they
typically will redistribute data to a large percentage of other processes
within the region. In recent years, the number of cores per node have largely
increased. On emerging systems, such as sapphire rapids, it is possible that
some sparsity patterns would be better optimized with the non-blocking
algorithm within each region. Additional possible optimizations include using
an MPI_Alltoallv for dense intra-region redistribution. However, as inter-
region communication bottlenecks the SDDE, intra-region communication
optimizations are beyond the scope of this paper.

##  V Results

This section analyzes all SDDE algorithms across a subset of matrices from the
Suitesparse Matrix Collective [18]. The subset of matrices is sequential in
order on the Suitesparse website when ordered by non-zeros and are four
matrices with non-zero counts closest to $25$ million. These matrices are
small enough to be stored on just two nodes but also large enough to be scaled
out to $64$ nodes.

The SDDE algorithms are analyzed on the Quartz supercomputer at Lawrence
Livermore National Lab. Each node of Quartz contains two Intel Xeon E5-2695 v4
processors, totalling 36 available CPU cores per node. All algorithms are
tested on both system versions of MPI available on Quartz, OpenMPI 4.1.2 and
Mvapich2 2.3.7. For simplicity, $32$ processes per node are used in each test,
allowing for power of $2$ process counts.

The MPIX_Alltoall_crs is analyzed across the subset of Suitesparse matrices
with OpenMPI in Figure 6 and Mvapich2 in Figure 5. For these tests, the sparse
dynamic data exchange communicates data sizes based on the sparsity pattern of
each matrix, so all messages contain only a single integer. For instance, if a
process $p$ holds non-zeros in columns $x$, $y$, and $z$, for which
corresponding rows are stored on process $q$, the SDDE will contain a message
with the integer $3$ sent from process $p$ to $q$, indicating that during
subsequent operations, $q$ should send $3$ elements to process $p$.

The costs of the various MPIX_Alltoall_crs algorithms are displayed as black
lines, while the maximum number of inter-node messages sent by any process are
displayed as red dots. The locality-aware SDDE methods communicate only the
aggregated number of inter-node messages while the personalized, non-blocking,
and RMA approaches all communicate the standard message count.

For all matrices, the maximum number of inter-node messages is greatly
decreased with aggregation. At large scales, the locality-aware non-blocking
SDDE typically outperforms the other methods due to the lack of reduction
along with the low message count. Furthermore, the scalability of this method
is greatly improved over the other approaches, allowing for increased
efficiency at even larger scales. Note, some node counts have no results on
OpenMPI due to a UCX error during RMA communication.

Figures 8 and 7 analyze the costs of the various MPIX_Alltoallv_crs algorithms
using OpenMPI and Mvapich2, respectively. The variable-sized SDDE methods form
the communication pattern that is required to perform sparse matrix operations
with the given matrices. As a result, the message sizes vary, with each
message containing indices that should be communicated during subsequent
sparse matrix operations. For example, if a process $p$ holds non-zeros in
columns $x$, $y$, and $z$, for which corresponding rows are stored on process
$q$, these SDDEs will consist of process $p$ sending a message of size $3$ to
$q$ containing the numbers $x$, $y$, and $z$. Standard and aggregated message
counts are equivalent to those within the MPIX_Alltoall_crs tests. However,
rather than containing only a single integer, each message contains a number
of integers equal to the message count, or the number of edges between the two
processes. Note, there are no RMA results as the RMA approach only applies to
MPIX_Alltoall_crs methods.

The locality-aware non-blocking SDDE method outperforms the others at scale
for all matrices other than dielFilterV2clx, which has the smallest message
count. For the other matrices, the locality-aware non-blocking method improves
the performance at $64$ nodes by up to 20x. Furthermore, the scalability of
the locality-aware non-blocking algorithm is best for all matrices, with the
SDDE algorithm scaling to at least $32$ or $64$ nodes for all matrices. As a
result, this algorithm improves not only performance at large scales, but also
allows for sparse matrix operations and the linear solvers that rely on them
to be efficiently performed at increasingly large scales.

##  VI Conclusions and Future Directions

The optimal sparse dynamic data exchange algorithm varies with communication
pattern, problem size, process count, and architecture. As per-region process
counts increase, the potential for locality-aware optimizations increases. The
impact of locality-aware aggregation increases with the number of processes
per message, and a large number of messages are expected to have a common
region of destination. Parallel applications with large irregular
communication constraints are expected to achieve significant benefits from
locality-aware SDDE algorithms at strong scaling limits. When tested across
$32$ nodes, with $32$ processes per node, the locality-aware approach achieves
up to 20x speedup over the second best approach for matrices with high message
counts, while incurring slowdown for matrices that require few messages. This
impact is seen for both the MPIX_Alltoall_crs and MPIX_Alltoallv_crs methods
across both system installations of MPI. While this paper did not explore
locality-aware aggregation for the RMA method, similar concatenation
strategies could be used within MPI_Puts to reduce the synchronization
overheads as well as communication costs.

Future performance models are needed to dynamically select the optimal SDDE
algorithm for given sparsity patterns on each existing and emerging
architecture. Furthermore, as heterogeneous architectures increase in
prevalence, the SDDE algorithm should be adapted for these systems. There are
an increasingly large number of possible paths of data movement on emerging
systems, including GPUDirect, copying to a single CPU, and copying portions of
data to all available CPU cores, which should be combined with existing
algorithms to optimize on emerging heterogeneous architectures.

## References

  * [1] A. Bienz, D. Schafer, and A. Skjellum, “Mpi advance : Open-source message passing optimizations,” in _EuroMPI’23: 30th European MPI Users’ Group Meeting_ , 2023. [Online]. Available: https://eurompi23.github.io/assets/papers/EuroMPI23_paper_33.pdf
  * [2] B. Vastenhouw and R. H. Bisseling, “A two-dimensional data distribution method for parallel sparse matrix-vector multiplication,” _SIAM Review_ , vol. 47, no. 1, pp. 67–95, 2005. [Online]. Available: https://doi.org/10.1137/S0036144502409019
  * [3] A. Lazzaro, J. VandeVondele, J. Hutter, and O. Schütt, “Increasing the efficiency of sparse matrix-matrix multiplication with a 2.5d algorithm and one-sided mpi,” in _Proceedings of the Platform for Advanced Scientific Computing Conference_ , ser. PASC ’17. New York, NY, USA: Association for Computing Machinery, 2017. [Online]. Available: https://doi.org/10.1145/3093172.3093228
  * [4] Y. Ishikawa, K. Nakajima, and A. Hori, “Revisiting persistent communication in mpi,” in _Recent Advances in the Message Passing Interface_ , J. L. Träff, S. Benkner, and J. J. Dongarra, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 296–297.
  * [5] M. Hatanaka, A. Hori, and Y. Ishikawa, “Optimization of mpi persistent communication,” in _Proceedings of the 20th European MPI Users’ Group Meeting_ , ser. EuroMPI ’13. New York, NY, USA: Association for Computing Machinery, 2013, p. 79–84. [Online]. Available: https://doi.org/10.1145/2488551.2488566
  * [6] R. E. Grant, M. G. F. Dosanjh, M. J. Levenhagen, R. Brightwell, and A. Skjellum, “Finepoints: Partitioned multithreaded mpi communication,” in _High Performance Computing_ , M. Weiland, G. Juckeland, C. Trinitis, and P. Sadayappan, Eds. Cham: Springer International Publishing, 2019, pp. 330–350.
  * [7] M. G. Dosanjh, A. Worley, D. Schafer, P. Soundararajan, S. Ghafoor, A. Skjellum, P. V. Bangalore, and R. E. Grant, “Implementation and evaluation of mpi 4.0 partitioned communication libraries,” _Parallel Computing_ , vol. 108, p. 102827, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0167819121000752
  * [8] J. L. Träff and A. Rougier, “Mpi collectives and datatypes for hierarchical all-to-all communication,” in _Proceedings of the 21st European MPI Users’ Group Meeting_ , ser. EuroMPI/ASIA ’14. New York, NY, USA: Association for Computing Machinery, 2014, p. 27–32. [Online]. Available: https://doi.org/10.1145/2642769.2642770
  * [9] X. Luo, W. Wu, G. Bosilca, Y. Pei, Q. Cao, T. Patinyasakdikul, D. Zhong, and J. Dongarra, “Han: a hierarchical autotuned collective communication framework,” in _2020 IEEE International Conference on Cluster Computing (CLUSTER)_ , 2020, pp. 23–34.
  * [10] A. Bienz, W. D. Gropp, and L. N. Olson, “Node aware sparse matrix–vector multiplication,” _Journal of Parallel and Distributed Computing_ , vol. 130, pp. 166–178, 2019.
  * [11] ——, “Reducing communication in algebraic multigrid with multi-step node aware communication,” _The International Journal of High Performance Computing Applications_ , vol. 34, no. 5, pp. 547–561, 2020.
  * [12] M. Hidayetoğlu, T. Bicer, S. G. de Gonzalo, B. Ren, V. De Andrade, D. Gursoy, R. Kettimuthu, I. T. Foster, and W.-m. W. Hwu, “Petascale xct: 3d image reconstruction with hierarchical communications on multi-gpu nodes,” in _SC20: International Conference for High Performance Computing, Networking, Storage and Analysis_ , 2020, pp. 1–13.
  * [13] B. Priest, T. Steil, G. Sanders, and R. Pearce, “You’ve got mail (ygm): Building missing asynchronous communication primitives,” in _2019 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW)_ , 2019, pp. 221–230.
  * [14] J. L. Träff and S. Hunold, “Decomposing mpi collectives for exploiting multi-lane communication,” in _2020 IEEE International Conference on Cluster Computing (CLUSTER)_ , 2020, pp. 270–280.
  * [15] A. Buluç, H. Meyerhenke, I. Safro, P. Sanders, and C. Schulz, _Recent Advances in Graph Partitioning_. Cham: Springer International Publishing, 2016, pp. 117–158. [Online]. Available: https://doi.org/10.1007/978-3-319-49487-6_4
  * [16] G. Shipman, “Cellar,” https://github.com/lanl/CELLAR, 2022.
  * [17] T. Hoefler, C. Siebert, and A. Lumsdaine, “Scalable communication protocols for dynamic sparse data exchange,” _SIGPLAN Not._ , vol. 45, no. 5, p. 159–168, jan 2010. [Online]. Available: https://doi.org/10.1145/1837853.1693476
  * [18] T. A. Davis and Y. Hu, “The university of florida sparse matrix collection,” _ACM Trans. Math. Softw._ , vol. 38, no. 1, dec 2011. [Online]. Available: https://doi.org/10.1145/2049662.2049663