Neural Splines: Fitting 3D Surfaces with Infinitely-Wide Neural Networks
Francis Williams, Matthew Trager, Joan Bruna, Denis Zorin
Introduction
Estimating a 3D surface from a scattered point cloud is a classical and important problem in computer vision and computer graphics. In this task, the input is a set of 3D points sampled from an unknown surface and, possibly, normals of the surface at those points. The goal is to estimate a representation of the complete surface from which the input samples were obtained, for example, a polygonal mesh or an implicit function. This problem is challenging in practice: It is inherently ill-posed, since an infinite number of surfaces can interpolate the data. Furthermore, the input 3D points are often incomplete and noisy, as they are acquired from range sensors such as LIDAR, structured light, and laser scans. Ideally, the recovered surface should not interpolate noise but preserve key features and surface details.
Many early surface reconstruction techniques consider a kernel formulation of the surface reconstruction problem, using translation-invariant kernels such as biharmonic RBFs , Gaussian kernels , or compactly supported RBFs . Currently, the most widely used method for surface reconstruction is Screened Poisson Surface Reconstruction , which solves a variant of the Poisson equation to find an implicit function whose zero-level set produces an approximating surface. We show in this paper that this method can also be viewed as a kernel method for a particular choice of kernel.
More recently, many papers have used neural networks to represent an implicit function or a local chart in a manifold atlas as a means of reconstructing a surface . These methods can be integrated into a data-driven learning pipeline or directly applied in the so called “overfitting” regime, where a massively overparameterized (\ie, more parameters than input points) neural network is fitted to a single input point cloud as a functional representation for a surface. Empirical evidence has shown that these methods enjoy some form of “implicit regularization” that biases the recovered surface towards smoothness. Moreover, employing early stopping in gradient descent can prevent these approaches from interpolating noise.
Under certain parameter initializations, infinitely-wide neural networks de facto behave as kernel machines , defining Reproducing Kernel Hilbert Spaces (RKHS), whose kernel is obtained by linearizing the neural network mapping around its initialization. While the kernel regime simplifies non-linear neural network learning into a convex program and provides a simple explanation for the successful optimization of overparameterized models, it cannot explain the good generalization properties observed for high-dimensional problems due to the inability of the RKHS to approximate non-smooth functions . However, the situation for low-dimensional problems such as surface reconstruction is entirely different, and RKHS can provide powerful characterizations of regularity. In this context, shows that in the univariate case the RKHS norm associated with a wide ReLU network is a weighted curvature, and leads to cubic spline interpolants. In higher dimensions, similar (albeit more complex) characterizations of the RKHS norm exist (see also Proposition 2). In order to assess the benefits of neural networks on such low-dimensional problems, it is thus important to first understand their linearized counterparts, given by their associated kernel machines.
In this work, we demonstrate that in fact kernels arising from shallow ReLU networks are extremely competitive for 3D surface reconstruction, achieving state-of-the-art results: outperforming classical methods as well as non-linear methods based on far more complex neural network optimization.
Kernels provide many advantages over neural networks in our context: (i) They are well understood theoretically. (ii) Kernel regression boils down to solving a linear system, and avoids gradient descent that suffers from slow convergence. (iii) Kernel-based interpolants are represented using a number of parameters that is linear in the size of the input, whereas overparameterized neural networks require many more parameters than points. (iv) The inductive bias of kernel methods can be characterized explicitly via the RKHS norm (Section 3.3). (v) Kernel methods lend to scalable and efficient implementations (Section 3.5). We provide explicit expressions for two kinds of infinite-width ReLU kernels, their derivatives, and their corresponding RKHS norms. We further argue that these kernels can be viewed as a multidimensional generalization of cubic spline interpolation in 1D. Moreover, we show that Poisson Surface Reconstruction can itself be viewed as a kernel method and give an expression for its RKHS norm, suggesting that kernels are a broad framework which enable the rigorous understanding both traditional and modern surface reconstruction techniques.
Methods for 3D surface reconstruction can mostly be divided by their choice of surface representation: These are (i) the zero level set of a volumetric scalar function , (ii) the fixed point of a projection operator onto locally fitted patches , (iii) a mesh connecting the input points , (iv) a collection of local parametric maps , or (v) the union of parametric shapes . For a recent up-to-date survey of surface reconstruction techniques, we refer the reader to .
The “random feature” kernels used this paper arise from training the top-layer weights in two-layer networks of infinite width and were described in . More recently, a lot of work has focused on a different kernel, known as the “neural tangent kernel” , that linearly approximates the training of both layers of a neural network. We chose to use random feature kernels since, when the input dimension is small, typical initializations in neural networks lead to mainly training the top layer weights. More broadly, the function spaces associated with shallow neural networks were studied in .
Neural Spline Formulation
We begin by introducing the Neural Spline family of functions through the lens of finite-width shallow neural networks in Section 2.1, and turn to the infinite-width limit formulation in Section 2.2. While our approach is simple and our derivations straightforward, we differ slightly from a standard kernel regression setting since we employ a multi-dimensional kernel that includes the gradient of the fitted function. This formulation allows us to fit points and normals simultaneously.
We first assume that the model is a two-layer ReLU neural network with neurons, but we keep the bottom layer weights fixed from initialization:
In typical initialization schemes for neural networks (\eg, Kaiming He initalization ), the weights of each layer are initialized with a variance that is inversely proportional to the number of inputs of that layer. Our choice of fixing is motivated by the fact that using a standard two-layer network, for a given target model and as , only the top layer weights tend to vary throughout training .
This minimizer is given explicitly by where denotes the Moore-Penrose pseudo-inverse of and
where for compactness we used , and .
2 Infinite-Width Kernel
As the number of neurons tends to infinity, the kernel (4) converges to defined by
The uniform distribution (8) corresponds to the default initialization of linear layers in PyTorch and, as we argue in Section 3.1, it leads to a direct generalization of cubic spline interpolation. However, the Gaussian initialization actually leads to simpler analytical expressions for and produces almost the same results. We refer to Appendix B for a discussion and comparison of the two distributions.
Discussion
The expression assuming is obtained by swapping and . For fixed , the map is piecewise cubic and twice continuously differentiable (). This implies that kernel regression with (10) yields cubic spline interpolation. Applying the Neural Spline objective (1) with derivative constraints at the samples in also yields a piecewice cubic interpolant, although this curve is in general only . For other distributions of and , the kernel is no longer cubic, but the norm in the RKHS is a weighed norm of curvature (see for details). In this sense, our approach with the initialization (8) can be viewed as a multi-dimensional version of spline interpolation.
2 Regularization and Robustness to Noise
To deal with noisy data, we can optionally add a simple regularizer term to our formulation that corresponds to penalizing the RKHS norm of the interpolant (“kernel ridge regression”). Concretely, we replace (6) with
for . The regularizer term affects the spectrum of the Gram matrix by smoothing its smallest eigenvalues, with a similar effect to early stopping in gradient descent (see \eg ). Figure 2 shows examples of applying this regularizer on 2D and 3D problems.
3 Inductive Bias of Interpolants
Our interpolant belongs to the Hilbert space given by
where is a measure over the weights (\eg, (8) or (9)), and can be viewed as a continuous analog of the outer-layer weights of the finite network given by (2). The inductive bias of Neural Splines is thus determined by the RKHS norm in (12) since our method outputs the interpolating function which minimizes that norm. As noted in , if and then, by differentiating twice, we note that the Laplacian of is given by
By comparing (12) and (13), we see that the RKHS norm and the Laplacian of are closely related. More precisely, the Laplacian is the Dual Radon Transform of . Under certain assumptions, (13) can be inverted, yielding an explicit expression for in terms of . Intuitively, this shows that bounding the RKHS norm imposes a constraint on the Laplacian of , and thus encourages to be smooth. We report the following statement from and we refer to Appendix D for more details.
Let , and the constant . If we assume that holds, then
where is the Radon Transform of .
4 Poisson Surface Reconstruction as a Kernel
We cast Screened Poisson Surface Reconstruction in kernel form to facilitate comparisons. In its simplest form, Poisson reconstruction, extracts the level set of a smoothed indicator function determined as the solution of
The vector field is obtained by interpolating the normals using a fixed-grid spline basis and barycentric coordinates of the sample points with respect to the grid cell containing it. This is equivalent to using a non-translation invariant non-symmetric locally-supported kernel :
where .
To study the qualitative properties of this kernel, we replace with a radial kernel (see Appendix C) which has qualitatively similar behavior (see Figure 3). Since both and the Laplace kernel are radial functions, their convolution is also radial, yielding a translation-invariant radial approximation of :
The RKHS norm of the corresponding to the approximate Poisson kernel is
where is the Fourier transform.
We discuss the kernel formulation of Poisson Reconstruction in more detail in Appendix C.
5 Fast and Scalable Implementation
We provide a fast and scalable implementation of Neural Spline kernels based on FALKON , a recently-proposed solver for kernel-ridge-regression which runs in parallel on the GPU. While naïve kernel ridge-regression with points requires solving and storing an dense linear system, FALKON uses conjugate gradient descent requiring only storage, and convergence (though in practice we find that even for very large inputs, we converge in fewer than 10 iterations). To speed up convergence, FALKON can optionally store an preconditioner matrix in CPU memory (where , see paragrph below). To maximize performance and reduce memory overhead, we rely on KeOps to evaluate kernel matrix-vector products symbolically on the GPU, which means our implementation uses only a small constant amount GPU memory and can be readily used on commodity hardware. Section 4.4 compares the performance of our implementation against other state of the art surface reconstruction techniques. We note that in principle, low-dimensional kernel methods can be accelerated using fast multipole-based approaches (in particular, in the context of 3D surface reconstruction this was used in ); this yields optimal time complexity, for dense matrix-vector multiplication.
Full kernel ridge regression predicts a function which is supported on every input point as in (6), requiring coefficients to store the resulting function. We rely on Nyström sampling to instead produce a kernel function which is supported on a small -sized subset of the input points (while still minimizing a loss on all the points). This is equivalent to approximating the kernel matrix with a low rank linear system. To choose Nyström samples, we leverage the geometric nature of our problem and select these by downsampling the input point cloud to have a blue-noise distribution using Bridson’s algorithm . We demonstrate the effect of varying the number of Nyström samples qualitatively in Figure 4.
Experiments and Results
We now demonstrate the effectiveness of Neural Splines on the task of surface reconstruction. For all the experiments in this section, we used the analytical form of the kernel (5) with Gaussian initialization (8). Appendix A.7 compares the uniform (8) and Gaussian (9) kernels, showing almost no measurable difference in the reconstructions produced by either. We compare the empirical and analytical kernels in detail in Appendix A.6. An implementation of Neural Splines is available at https://github.com/fwilliams/neural-splines.
We performed a quantitative evaluation on a subset of the Shapenet dataset to demonstrate that the inductive bias of Neural Splines is particularly effective for reconstructing surfaces from sparse points. We chose 1024 random points and normals sampled from the surface of 20 shapes per category across 13 categories (totalling 260 shapes). Using this dataset, we compared our method against Implicit Geometric Regularization (IGR) , SIREN , Fourier Feature Networks , Biharmonic RBF (Biharmonic) , SVM surface modelling (SVR) , and Screened Poisson Surface Reconstruction (Poisson) . The first three techniques are modern neural network based methods, while the latter three techniques are classical methods based on kernels or solving a PDE. As criteria for the benchmark, we consider the Intersection over Union (IoU) and Chamfer Distance between the reconstructed shapes and the ground truth shapes. The former metric captures the accuracy of the predicted occupancy function, while the latter metric captures the accuracy of the predicted surface. Under both metrics, our method outperforms all other methods by a large margin. Table 1 reports quantitative results for the experiment and Figure 5 shows visual results on a few models. We report per-category results in Appendix A.4 and show many more figures in Appendix A.3.
Several of the above methods have parameters which can be tuned to increase performance. To ensure a fair comparison, we ran parameters sweeps for each Shapenet model for these methods, reporting the maximimum of each metric under consideration (see Appendix A.1 for a detailed description). For our method we did not tune parameters. We used no regularization and 1024 as Nyström samples for all models.
2 Surface Reconstruction Benchmark
We evaluated our method on the Surface Reconstruction Benchmark which consists of simulated noisy range scans (points and normals) taken from 5 shapes with challenging properties such as complex topologies, sharp features, and small surface details. We evaluate our method against Deep Geometric Prior (DGP) , Implicit Geometric Regularization (IGR) , SIREN , and Fourier Feature Networks (FFN) . We remark that DGP establishes itself as superior to a dozen other classical methods on this benchmark and IGR furhter outperforms DGP. As in and , Table 2 reports the Hausdorff () and Chamfer () distances between the reconstruction and ground-truth. We also report the one sided Hausdorff () and Chamfer () distances between the scan and the reconstruction, which measures how much the reconstructions overfits noise in the input. Our reconstructions are quantitatively closer to ground truth on all but one model. Figure 6 shows visual examples of a few models from the benchmark. All models are shown in Appendix A.3. As with the Shapenet benchmark, we did parameter sweeps for those methods which have tunable parameters, choosing the best model for each metric under consideration. See Appendix A.2 for details.
3 Large Scale Reconstruction of Full Scenes
Figure 7 shows a full scene consisting of 9 million points reconstructed using our method. This is the same scene used in , and contains many thin features that are difficult to reconstruct (\eg, curtains, plant, and lampshade). To generate the input point cloud for this experiment, we densely sampled a mesh extracted from a occupancy grid of the scene. To perform the reconstruction, we subdivided the space into 8x8x8 cells containing between 10k and 500k samples and reconstructed each cell interdependently using up to 15k Nyström samples. The whole process takes 1.5 hours on a machine with an NVIDIA 1080Ti GPU.
4 Timing and Performance
Table 3 compares the average running time and GPU memory usage of our method and others when reconstructing point clouds from the Surface Reconstruction Benchmark described in Section 4.2. To ensure a fair comparison we ran all neural network based methods for 5000 iterations. We do not report exact CPU memory usage for the methods since it is hard to measure but we remark that, by observation, system memory usage never exceeded 3GiB for any of the methods. Appendix A.5 reports the running times and memory usages for individual models in the benchmark.
5 Our Method versus Neural Networks
Implicit Geometric Regularization demonstrates that ReLU networks have a natural inductive bias making them good at reconstructing shapes. However, methods based on ReLU networks suffer from slow convergence (see Figure 8). SIREN drastically improves convergence speed by replacing ReLU with sinusoidal activations and a clever initialization, however, SIREN occasionally underfits on sparse inputs and requires an additional loss computed on points in the volume around the shape to prevent reconstruction artefacts. This loss increase runtimes, and suggest that SIREN’s inductive bias may not be ideal for sparse reconstruction tasks. By projecting input points onto random Fourier features, Fourier Feature Networks (FFNs) present a principled approach rooted in Kernel methods to control the inductive bias of the solution as well as convergence speed. While FFNs are capable of producing high quality solutions, they are sensitive to the choice of feature distribution when data is sparse (see Figure 9) and thus require tuning to work well. In contrast, our technique converges in seconds (Section 4.4), has a well suited inductive bias for shape representation (Section 3.3), and requires minimal parameter tuning.
Conclusion and Future Work
We have shown that Neural Spline kernels arising from infinitely wide shallow ReLU networks are very effective tools for 3D surface reconstruction, outperforming state-of-the-art methods while being computationally efficient and conceptually simple. In a sense, our work bridges the gap between traditional reconstruction methods and modern neural network based methods by leveraging the deep connection between neural networks and kernels.
We remark that our kernel formulation is fully differentiable. In the future, we hope to integrate Neural Splines into deep learning pipelines and apply them to other 3D tasks such as shape completion and sparse reconstruction. On the theory side, we would like to investigate and compare the approximation properties of different kernels (in particular those arising from infinite width sinusoidal networks) in the context of 3D reconstruction.
This work is partially supported by the Alfred P. Sloan Foundation, NSF RI-1816753, NSF CAREER CIF 1845360, NSF CHS-1901091, and Samsung Electronics.
References
Appendix A Additional Experiments
Many of the methods we compared against on Shapenet have tunable parameters which can drastically alter the quality of reconstructed outputs. To ensure a proper comparison, we ran sweeps over these parameters where appropriate choosing the best reconstruction for each model under both metrics (Chamfer and IoU). For our method we did no parameter sweeps, using no regularization and 1024 Nyström samples for each model in the dataset. We describe the exerimental methodology for each method in the benchmark below.
We trained each model for 5k iterations with Adam and a learning rate of 0.001 using the same parameters and architecture as proposed in the original paper. We included the normals in the loss with the parameter set to 1. The Eikonal regularization term was set to 0.1. While IGR can slightly improve by using a very large number of iterations. doing so is prohibitively slow over many models. Figure 8 motivates our choice of iterations, demonstrating only a slight improvements between 5k and 100k iterations (the latter which required 2 hours of fitting on a NVIDIA-1080-Ti GPU).
Screened Poisson Surface Reconstruction [31]
We considered every possible combination of the following parameters: the octree depth in , the point weight (which controls the degree to which the method interpolates the input) in . Since all the shapes in the benchmark are watertight meshes, we used Dirichlet boundary constraints for the reconstruction.
SIREN [40]
We trained used a SIREN network with 5 hidden layers each containing 256 neurons for 5000 iterations, using a learning rate of 1e-4 with Adam. We used the same loss for shapes as in the original SIREN paper, sampling an additional (where is the number of input points) points in an axis aligned bounding box whose diagonal is 10% larger than the object’s bounding box. We minimize the same loss for shapes as the SIREN paper (see section 4.2 in the original SIREN paper).
Fourier Feature Networks [42]
We used an 8-layer ReLU MLP with 256 Fourier features sampled from a Gaussian distribution. This is the same architecture as the shape representation experiment in the original paper. For each model in the benchmark, we did a parameter sweep on the variance of the Gaussian distribution, considering . The range of parameters was chosen by empirical verification on 3 models from the airplanes, benches, and cars categories.
SVR [22]
As in the original paper, we use a Gaussian kernel to perform support vector regression. To generate occupancy samples, we augmented the input points with an ”inside“ and ”outside“ point by perturbing them by along the normal at that point. We used forall the models. For each model we did a joint parameter sweep over the regularization parameter and the variance parameter . The range of parameters was chosen by empirical verification on 3 models from the airplanes, benches, and cars categories.
Biharmonic RBF [9]
To generate occupancy samples, we augmented the input points with an ”inside“ and ”outside“ point by perturbing them by along the normal at that point. We used for all the models. The biharmonic function is very simple where and does require tuning parameters.
A.2 Detailed Description of Surface Reconstrucion Benchmark Experiment
A.3 Additional Figures
Figures 10 and 11 show at least one model reconstructed from each ShepeNet category using our method, Implicit Geometric Regularization (IGR) , SIREN , Fourier Feature Networks (FFN) , Screened Poisson Surface Reconstruction , Biharmonic RBF , and Support Vector Regression (SVR) . Figure 12 shows the reconstructions of all the models from the Surface Reconstruction Benchmark using our method, IGR , SIREN , and Fourier Feature Networks .
A.4 Quantitative Results Per ShapeNet Class
Tables 4 show the per ShapeNet category IoU and Chamfer distance statistics for the benchmark described in Section 4.1.
A.5 Per Model Performance Numbers
Table 5 shows the runtime and GPU usage required to reconstruct each model in the Surface Reconstruction Benchmark. For our model, we used 15k Nyström samples and a regularization of 1e-11. We do not report CPU memory usage since it is hard to profile exactly, however we observed that none of the methods used more than 4GiB of CPU memory. All timings were done on a machine with a single NVIDIA-V100 GPU with 16GiB of VRAM, 32GiB of CPU RAM, and and an 8 core Intel Xeon processor.
A.6 Empirical versus Analytical Kernel
Figure 14 compares results using the empirical kernel with neurons and using the analytical kernel. Figure 13 shows the convergence of the empirical Kernel to the analytic one as the number of neurons grows.
A.7 Quantitative Comparison between Gaussian and Uniform Kernels
Table 6 shows a quantitative comparison between Neural Splines using the Gaussian initialization (9) and the Uniform initialization (8) on the benchmark described in Section (4.1). The results in both cases are very close to each other in both Chamfer distance and in IoU.
Appendix B Derivation of the Infinite Width Kernels
We derive an explicit expression for the kernel in (5) in the case of uniform initialization (8). We first prove the following Lemma which we will use in our calculations.
where (for , etc.), is such that , and
We now consider the integral . For any two indices , we have
This is now a product of one-dimensional integrals. Since if is odd, we have that the integral (18) vanishes if . If instead , we use the fact that
This proves the diagonal part in our expression for . All remaining terms as well as the two integrals and follow from very similar (and slightly simpler) calculations. ∎
We now apply Lemma 5 to compute the kernel with the uniform initialization (8).
where , \tau=\arctan\bigg(\frac{||x||-||x^{\prime}||\cos(\alpha)}{||x^{\prime}||\sin(\alpha)}\bigg{missing}), is such that , and
The idea is to compute the integral with respect to the bias term and then split the result into homogeneous expressions where Lemma 5 can be applied. In particular, assuming that :
B.2 Gaussian Initialization
The Gaussian initialization (9) yields the following simpler formula for . The first term is well known and is derived in . The second and third terms are easily derived by taking derivatives of the first term.
If and , then
Appendix C Poisson Surface Reconstruction Kernel
In its simplest form, Poisson reconstruction of a surface , extracts the level set of a smoothed indicator function determined as the solution of
The vector field is obtained by interpolating the normals using a fixed-grid spline basis and barycentric coordinates of the sample points with respect to the grid cell containing it. This is equivalent to using a non-translation invariant non-symmetric locally-supported kernel :
Then . where only 8 terms corresponding to the vertices of the grid cube containing are nonzero. This yields the following expression for the kernel corresponding to Poisson reconstruction,
, the convolution of the Laplacian kernel and the gradient of . Using the identity , we can write this as the gradient of , defined as
To make it easier to understand the qualitative behavior of the kernel, replacing with a radial kernel , with qualitatively similar behavior (see Figure 3) yields a translation-invariant radial approximation of the kernel , as the convolution of two radial kernels is a radial function.
As both and the Laplace kernel are radial functions, their convolution is also radial. It can be expressed in a more explicit form using the relation between Fourier and Hankel transforms for radial functions. For , the Hankel transform is related to Fourier transform by
The Hankel transform is an involution, so the relationship for the inverse Fourier transform is similar. Writing , we obtain the expression for the radial convolution in terms of one-dimensional integrals,
where we use . and is just the gradient of this, i.e., a derivative times .
The RKHS norm for the space corresponding to this kernel is given by
with obtained using the Hankel transforms as above.
Appendix D RKHS Norm of the Neural Spline Kernel
We now discuss how in (12) is related to the Laplacian of the fucntion. If we make the mild assumption that our functions contain a linear and bias term (Lemma 8), then is the Radon Transform of the laplacian of the function. Thus, the least norm minimizers of the least squares problem (1) are related to the laplacian of the function and the RKHS norm corresponds to the integral of the laplacian over hyperplanes in the domain. In our experiments, we added an option to include the linear and bias terms to the solution. They appear to have no effect on the final reconstruction. The derivation below is borrowed from .
Then, can always be rewritten as
We can split the integral in into even and odd parts:
where and are the even and odd parts of respectively. Observing that and , we have that
Using Lemma 8, we will consider without loss of generality, neural networks of the form (25) with even measures . We now give a few useful definitions and lemmas.
where is a measure on the -hyperplane . Intuitively the Radon transform represents a function in terms of its integrals along all possible hyperplanes.
Since the hyperplane is the same as the hyperplane , the Radon transform is an even function. i.e. .
The Radon Transform satisfies the intertwining property. i.e. for any positive integer
where is the Radon Transform of . In particular, for ,
The Laplacian of in is (25)
which is precisely the Dual Radon Transform of . Since is even, and assuming it decays rapidly with , we can invert it using Lemma 12 yielding
The RKHS norm of the function is