Differentiable Spline Approximations
Year: 2022

Abstract:
The paradigm of differentiable programming has significantly enhanced the scope of machine learning via the judicious use of gradient-based optimization. However, standard differentiable programming methods (such as autodiff) typically require the machine learning models to be differentiable, limiting their applicability. Our goal in this paper is to use a new, principled approach to extend gradient-based optimization to functions well modeled by splines, which encompass a large family of piecewise polynomial models. We derive the form of the (weak) Jacobian of such functions and show that it exhibits a block-sparse structure that can be computed implicitly and efficiently. Overall, we show that leveraging this redesigned Jacobian in the form of a differentiable “layer” in predictive models leads to improved performance in diverse applications such as image segmentation, 3D point cloud reconstruction, and finite element analysis. We also open-source the code at https://github.com/idealab-isu/DSA.

heading: 1 Introduction
text: Motivation: Differentiable programming has been a paradigm shift in algorithm design. The main idea is to leverage gradient-based optimization to optimize the parameters of the algorithm, allowing for end-to-end trainable systems (such as deep neural networks) to exploit structure in data and achieve better performance. This approach has found use in a large variety of applications such as scientific computing [Innes, 2020; Innes et al., 2019; Schafer et al., 2020], image processing [Li et al., 2018a], physics engines [Degrave et al., 2017], computational simulations [Alnæs et al., 2015], and graphics [Li et al., 2018b; Chen et al., 2019]. One way to leverage differentiable programming modules is to encode additional structural priors as “layers” in a larger machine learning model. Inherent structural constraints such as monotonicity, or piecewise constancy, are particularly prevalent in applications such as physics simulations, graphics rendering, and network engineering. In such applications, it may be beneficial to build models that obey such priors by design.
Challenges: For differentiable programming to work, all layers within the model must admit simple gradient calculations; however, this poses a major limitation in many settings. For example, consider computer graphics applications for rendering 3D objects [Kindlmann et al., 2003; Gross et al., 1995; Loop and Blinn, 2006]. A common primitive in such cases is a spline (or a piecewise polynomial) function which either exactly or approximately interpolates between a discrete set of points to produce a continuous shape or surface. Similar spline (or other piecewise polynomial) approximations arise in
⇤†Equal contribution. 35th Conference on Neural Information Processing Systems (NeurIPS 2021).
partial differential equation (PDE) solvers [Hughes et al., 2005], network flow problems [Balakrishnan and Graves, 1989], and other applications.
For such problems, we would like to compute gradients “through” operations involving spline approximation. However, algorithms for spline approximation often involve discontinuous (or even discrete) co-domains and may introduce undefined (or even zero) gradients. Generally, embedding such functions as layers in a differentiable program, and running automatic differentiation on this program, requires special care. A popular solution is to relax these non-differentiable, discrete components into continuous approximations for which gradients exist. This has led to recent advances in differentiable sorting [Blondel et al., 2020; Cuturi et al., 2019], dynamic programming [Mensch and Blondel, 2018], and optimization [Djolonga and Krause, 2017; Agrawal et al., 2019; Deng et al., 2020].
Our contributions: We propose a principled approach for differentiable programming for spline functions without the use of continuous relaxation2. For the forward pass, we leverage fast algorithms for computing the optimal projection of any given input onto the space of piecewise polynomial functions. For the backward pass, we leverage a fundamental locality property in splines that every piece (or basis function) in the output approximation only interacts with a few other elements. Using this, we derive a weak form of the Jacobian for the spline operation and show that it exhibits a particular block-structured form. While we focus on spline approximation in this paper, our approach can be generalized to any algorithmic module with piecewise smooth outputs. Our specific contributions are as follows:
1. We propose the use of spline function approximations as “layers” in differentiable programs. 2. We derive efficient (nearly-linear time) methods for computing forward and backward passes for
various spline approximation problems, showing that the (weak) Jacobian in each case can be represented using a block sparse matrix that can be efficiently used for backpropagation. 3. We show applications of our approach in three stylized applications: image segmentation, 3D point cloud reconstruction, and finite element analysis for the solution of partial differential equations.
Related Work Before proceeding, we briefly review related work. Extensions of autodiff: Automatic differentiation (autodiff) algorithms enable gradient computations over basic algorithmic primitives such as loops, recursion, and branch conditions [Baydin et al., 2018]. However, introducing more complex non-differentiable components requires careful treatment due to undefined or badly behaved gradients. For example, in the case of sorting and ranking operators, it can be shown that the corresponding gradients are either uninformative or downright pathological, and it is imperative the operators obey a ‘soft’ differentiable form. Cuturi et al. [2019] propose a differentiable proxy for sorting based on optimal transport. Blondel et al. [2020] improve this by proposing a more efficient differentiable sorting/ranking operator by appealing to isotonic regression. Berthet et al. [2020] introduce the use of stochastic perturbations to construct smooth approximations to discrete functions, and other researchers have used similar approaches to implement end-to-end trainable top-k ranking systems [Xie et al., 2020; Lee et al., 2020]. Several approaches for enabling autodiff in optimization have also been researched [Pogančić et al., 2020; Amos and Kolter, 2019; Agrawal et al., 2019; Mensch and Blondel, 2018].
Structured priors as neural “layers”: As mentioned above, one motivation for our approach arises from the need for enforcing structural priors for scientific computing applications. Encoding non-differentiable priors such as the solutions to specific partial differential equations [Sheriffdeen et al., 2019], geometrical constraints [Joshi et al., 2020; Chen et al., 2019], and spatial consistency measures [Djolonga and Krause, 2017] perform well but typically require massive amounts of structured training examples.
Spline approximation: Non-Uniform Rational B-splines (NURBS) are commonly used for defining spline surfaces for geometric modeling [Piegl and Tiller, 1997]. NURBS surfaces offer a high level of control and versatility; they can also compactly represent the surface geometry. The versatility of NURBS surfaces enables them to represent more complex shapes than Bèzier or B-splines. Several frameworks that leverage deep learning are beginning to use NURBS representations. Minto et al. [2018] use NURBS surfaces fitted over the 3D geometry as an input representation for the object
2While tricks such as straight-through gradient estimation [Bengio, 2013] also avoid continuous relaxation, they are heuristic in nature and may be inaccurate for specific problem instances [Yin et al., 2019].
classification task of ModelNet10 and ModelNet40 datasets. Erwinski et al. [2016] presented a neural-network-based contour error prediction method for NURBS paths. Fey et al. [2018] present a new convolution operator based on B-splines for irregular structured and geometric input, e.g., graphs or meshes. Very recently, Sharma et al. [2020] perform point cloud reconstruction to predict a B-spline surface, which is later processed to obtain a complete CAD model with other primitives “stitched” together.
Differentiable PDE solvers: With the advent of deep learning, there has been a recent rise in the development of differentiable programming libraries for physics simulations [Hu et al., 2019; Qiao et al., 2020]. Most often, the physics phenomena are represented using partial differential equations (PDEs) [Sanchez-Gonzalez et al., 2020; Holl et al., 2020]. Considerable effort has gone into designing physics-informed loss functions [Raissi et al., 2019; Raissi and Karniadakis, 2018; Kharazmi et al., 2021] whose optimization leads to desired solutions for PDEs. Due to space limitations, we defer to a detailed survey of this (vast) area by Cai et al. [2021].

heading: 2 Differentiable Spline Approximation
text: We now introduce our framework, Differentiable Spline Approximation (DSA), as an approach to estimate gradients over piecewise polynomial operations. Our main goal will be to estimate easy-to-compute forms of the (weak) Jacobian for several spline approximation problems, enabling their use within backward passes in general differentiable programs.
Setup. We begin with some basic definitions and notation. Let f 2 Rn be a vector where the ith element is denoted as fi. Let us use [n] = {1, 2, . . . , n} to denote the set of all coordinate indices. For a vector f 2 Rn and an index set I ✓ [n], let fI be the restriction of f to I , i.e., for i 2 I , we have fI(i) := fi, and fI(i) := 0 for i /2 I . Now, consider any fixed partition of [n] into a set of disjoint intervals I = {I1, . . . , Ik} where the number of intervals |I| = k. The `2-norm of f is written as kfk2 := pPn i=1 f 2 i while the `2 distance between f, g is written as kf gk2.
We first define the notion of a discretized k-spline. Note that the use of “spline" here is non-standard and somewhat more general than what is typically encountered in the literature. (Indeed, the spline concept used in computer graphics is a special instance of this definition; we explain further below.) Def. 2.1 (Discretized k-spline). A vector h 2 Rn is called a discretized k-spline with degree d if: (i) there exists a partition of [n] into k disjoint intervals I1, . . . , Ik; (ii) within each interval Ii, the coefficients of hj , j 2 Ii, can be perfectly interpolated by some polynomial function of degree d.
Let us illustrate this by an example. Suppose that d = 1 and k = 5. Then, h is a discretized k-spline with degree d if, in a “line plot” of the vector h (i.e., we interpolate the 2D points (j, hj) for all j 2 [n]), we see up to k = 5 distinct linear pieces. A different way to interpret this definition is that we start with a piecewise degree-d polynomial function H : R ! R with k = 5 pieces (with suitably defined knot points, which are the location of the intervals I), and evaluate H at any n equally spaced points in its domain. This gives us a vector h 2 Rn, which we call a discretized k-spline. In contrast with traditional splines, we allow H to be arbitrarily defined at the knot points and require no specific continuity or differentiability properties. Therefore, our definition encompasses all standard spline families (including interpolating/approximating splines such as smoothing-, cubic-, and B-splines).

heading: 2.1 Spline Approximation
text: Our focus in this paper is the problem of computing the best possible spline fit to a given set of data points (where both the parameters of the spline as well as the knot vectors are allowed to be variable).
We provide an algebraic interpretation of this problem. For a given vector space Rn, consider Skd , the set of all discretized k-splines with degree d. Since (standard) splines are vector spaces for a fixed set of knots, one can easily see that for any fixed partition of [n] into k subsets, the family of discretized k-splines is a k(d+ 1)-dimensional subspace of Rn. Now suppose that the knot indices are allowed to vary. The number of possible partitions is finite (of the order of
n k ), and therefore the set Skd is a
finite union of subspaces, or a nonlinear submanifold, embedded in Rn. Therefore, the problem of discretized k-spline approximation can be viewed as an orthogonal projection onto this nonlinear manifold. Consider any arbitrary vector x 2 Rn (we can think of (i, xi)
as a set of n data points to which we are trying to fit a k-spline). Then, the best k-spline fit to x (in the sense of `2 distance) amounts to solving the optimization problem:
F (x) = argmin h
1 2 kx hk22 = 1 2
nX
i=1
(xi hi)2 s.t. h 2 Skd (1)
This operation resembles standard spline regression. But it is strictly more general since this requires not only optimizing piecewise spline parameters but also the knot indices. Crucially, we note that F is both a non-differentiable and a non-convex map. Nevertheless, such an orthogonal projection can be computed in polynomial (in fact, nearly-linear) time [Jagadish et al., 1998; Acharya et al., 2015] using many different techniques, including dynamic programming. This forms the forward pass of our DSA “layer”.
Our first main conceptual contribution is a formal derivation of the backward pass of the orthogonal projection operation. Strictly speaking, the Jacobian is not well-defined due to the non-differentiable nature of the forward pass (owing to the non-differentiability built into the definition of the kspline). Therefore, we will instead be deriving the so-called “weak” form of the Jacobian (borrowing terminology from Blondel et al. [2020]).
We leverage two properties of the projection operation: (1) the output of the forward pass h corresponds to a partition of [n], that is, each element of hj corresponds to a single interval, Ij , and (2) within each interval, the least-squares operation is continuous and differentiable. The first property ensures that every element xi contributes to only a single piece in the output h. Given that the sub-functions from the piecewise partitioning function are smooth, we also observe that the size of each block corresponds to the size of the partition, Ii. Using this observation, we get:
Theorem 1. The Jacobian of the operation F with respect to x 2 Rn can be expressed as a block diagonal matrix, J 2 Rn⇥n, whose (s, t)th entry obeys:
Jx(F (x))(s, t) = @h(x)s @xt =
( @hIi (x)s
@xt if s, t 2 Ii
0 otherwise (2)
As a concrete instantiation of this result, consider the case d = 0. This is the case where we wish to best approximate the entries of x with at most k “horizontal” pieces, where the break-points are obtained during the forward pass3. Call this approximation h. Then, the Jacobian of h with respect to x forms the block-diagonal matrix J 2 Rn⇥n:
J =
2
664 J1 0 . . . 0 0 J2 . . . 0 ... ... . . .
... 0 0 . . . Jk
3
775 (3)
where all entries of each block, Ji 2 R|Ii|⇥|Ii| are constant and equal to 1/|Ii|, i.e., they are row/column-stochastic. Note that the sparse structure of the Jacobian allows for fast computation and that computing the Jacobian vector product JT ⌫ for any input ⌫ requires O(n) running time. As an additional benefit, the decoupling induced by the partition enables further speed up in computation via parallelization. See the Appendix for proofs, as well as derivations of similar Jacobians for k-spline approximation of any degree d 1, and generalization to 2D domains (surface approximation). In Section 3 we demonstrate the utility of this approach for a 2D segmentation (i.e., piecewise constant approximation) problem, similar to the setting studied in Djolonga and Krause [2017].
3In the data summarization literature, this class of functions is sometimes called k-histograms [Jagadish et al., 1998]

heading: 2.2 Differentiable NURBS
text: We now switch to a slightly different setting involving a special spline family known as non-uniform rational B-splines (NURBS), which are common in geometric modeling. Mathematically, a NURBS curve is a continuous function C : R ! R defined as follows. Construct any knot vector u (i.e. a non-decreasing sequence of real coordinate values) and fix degree d. Recursively define a sequence of basis functions, Ndi : R ! R computed using the Cox-de Boor formula:
Ndi (u) = u ui
ui+d ui Nd 1i (u) + ui+d+1 u ui+d+1 ui+1 Nd 1i+1 (u), N 0 i (u) = ⇢ 1 if ui  u  ui+1 0 otherwise
(4) for d = 1, 2, . . .. In the uniform case (where the knots are equally spaced), each Ndi can be viewed as being generated by recursively convolving a box function with Nd 1i . The non-uniform case cannot be written as a convolution, but the intuition is similar. With these basis functions in hand, the NURBS curve C is defined as the rational function:
C(u) =
Pn i=0 N
d i (u)wiPiPn
i=0 N d i (u)wi
, (5)
where Pi, i = 0, 1, . . . , t are called control points and wi are corresponding non-negative weights. The number of control points is related to the number of knots k and curve degree d as follows: k = t + d + 1. For simplicity, assume that all weights are equal to one. The basis functions in NURBS add up to one uniformly for each u (this is called the partition of unity property). Therefore:
C(u) = tX
i=0
Ndi (u)Pi, (6)
In summary, the NURBS curve is parametrically defined via the control points and the knot positions. This discussion is for 1D curves, but an extension to higher-order surfaces is conceptually similar.
Consider implementing NURBS as a differentiable “layer” where the inputs are the knot positions and control points. The forward pass through this layer simply consists of evaluating Equation 6 via the recursive Equation 4, and storing the various basis functions (and their spans) for further use.
However, the backward pass is a bit more complicated, once again due to the non-differentiable nature of C. The gradient with respect to the control point coordinates, P is straightforward since the mapping from P to C is linear. However, the gradient with respect to the knot positions, ui, is not well-defined due to the non-differentiable nature of the base cases of the recursion (which are box functions specified in terms of ui). Once again, we see that the non-differentiability of NURBS is built into its very definition, and this affects the numerics.
To resolve this, we propose the following approach to compute an (approximate) Jacobian of C. The main source of the issue is the derivative of the box-car function N0i (u) = 1[ui,ui+1) with respect to the knot points, which is not well defined. However, N0i (u) can be viewed as the difference between convolutions of the unit step function with ui and ui+1 , where is the Dirac delta defined over the real line. We smoothly approximate the delta function by a Gaussian function with small enough bandwidth hyperparameter : (ui) ⇡ g(u) = exp( (u ui)/2 2). This function is now differentiable with respect to ui, with g0(u) = u ui 2 g(u). Convolutions and differences are linear, and hence the derivative is the basis function times a multiplicative factor. Finally, a similar approach as the Cox-de Boor recursion (Equation 4) can be used to reconstruct the derivatives for all basis functions of higher order. See Algorithm 1 for pseudocode and the Appendix for details. Algorithm 1 Backward pass for NURBS Jacobian (for one curve point , C(u)) P0, U0: gradients of C w.r.t. P, U Initialize: P0,U0 ! 0 Retrieve uspan, Ndi , C(u) calculated during forward pass /* uspan is the index of knot position */ /* Ndi is the basis function of degree d */ /* C(u) is the evaluated curve point */ for h = 0 : d+ 1 do
P0uspan+h = N d h // easy since C is a linear function of P. U0uspan+h = N d h Uuspan+h // due to Gaussian approximation; see discussion below.
Let us probe the structure of this Jacobian a bit further. Suppose we evaluate the curve C at n arbitrary domain points. There are slightly less than k control points, and therefore the Jacobian is roughly of size n ⇥ O(k). However, due to the recursive nature of the definition of basis functions, the span (or support) of each basis function is small and only touches d+ 1 knots; for example, only 2 knots affect N0i , only 3 knots impact N1i , and so on. This endows a natural sparse structure on the Jacobian. Moreover, for a fixed order parameter d+ 1, the span is constant [Piegl and Tiller, 1997]; therefore, assuming evenly spaced evaluation points, we have the same number of nonzeros. Therefore, the Jacobian exhibits an interesting Toeplitz structure (unlike the block diagonal matrix in the case of Equation 3), thereby enabling efficient evaluation during any gradient calculations. We show below in Section 3 that automatic differentiation using this approach surpasses existing NURBS baselines.

heading: 2.3 Differentiable Finite Element PDE Solvers
text: Next, we see how spline approximations can be used to improve finite element analysis for solving PDEs. Popular recent efforts for solving PDEs using autodiff construct “physics-informed” solvers [Raissi et al., 2019; Raissi and Karniadakis, 2018], while other efforts have been made to utilize variational [Kharazmi et al., 2021] or adjoint-based derivative methods [Holl et al., 2020]. However, these approaches come with challenges while used in conjunction with autodiff packages, and gradient pathologies pose a major barrier [Wang et al., 2020].
Using our principles developed above, we propose an alternative PDE solution approach via differentiable finite elements. PDE solvers based on Finite Element Methods (FEM) are ubiquitous, and we provide a very brief primer here. Consider a domain ⌦ and a differential system of equations:
N [U(u)] = F (u), u 2 ⌦, (7) where N denotes the differential operator and U : ⌦ ! R is a continuous field variable; it is common to specify additional boundary constraints on U. The Galerkin method converts solving for the best possible U (which is a continuous variable) into a discrete problem by first looking at the weak form: R(U) = R ⌦ V [N (U) F ] du, where V is called a test function (and the weak form may involve some integration by parts), and rewriting this weak form in terms of a finite set of basis coefficients. A typical set of basis functions j is obtained by (piecewise) concatenation of polynomials, each defined over elements of a given partition of ⌦ (also called a mesh). Commonly used choices include Lagrange polynomials, defined by:
pri,d(u) = dX
r=1
Ur Y
0md m 6=r
u um ur um s.t. xr 2 [ 1, 1] (8)
where {u0, u1, . . . , ud} are a finite set of nodes (akin to control points in our above discussion, except in this case the splines interpolate the control points) and Ur is the corresponding coefficient. We use this collection of basis functions j to represent U:
U(u) = #nodesX
j=1
j(u)U d j (9)
and likewise for V . (The resemblance with Equation 5 above should be clear, and indeed NURBS basis functions could be an alternative choice.) Plugging the discrete coefficient representation Uc := {Ucj} into the definition of R, we get a standard Finite Element form,
R(Uc, V c) = B(Uc, V c) L(V c) (10) where B(Uc, V c) is the discrete form (bilinear for linear operators) that encodes the differential operator and L(v) is a linear functional involving the forcing function. For most PDE operators (including linear elliptic operators), one can form the energy functional by using U as the test function:
J(Uc) = 1
2 B(Uc,Uc) L(Uc). (11)
Optimization of this energy functional can now be performed using gradient-based iterations evaluated by automatic differentiation. This is a powerful approach since formal techniques exist (e.g., Galerkin Least Squares Bochev and Gunzburger [2009]) that reformulate weak forms of PDEs into equivalent energy functionals. The key aspect to note here is that differentiating “through" the differential operator N (embedded within B) requires derivative computations of the piecewise polynomial basis functions js, and therefore our techniques developed above are applicable.

heading: 3 Experiments
text: We have implemented the DSA framework (and its different applications provided below) by extending autograd functions in Pytorch. We also provide the capability to run the code using CUDA for GPU support. All the experiments were performed using a local cluster with 6 compute nodes and each node having 2 GPUs (Tesla V100s with 32GB GPU memory). All training was done using a single GPU. We summarize all our experiments in Table 1. Each experiment shown below is performed multiple times with different random seeds, and the average value with error bars is provided. Due to limited space, we provide three interesting applications of spline approximations here (see Appendix for additional examples). We have also open-sourced the code at https://github.com/idealab-isu/DSA.
Image segmentation: We begin with implementing a 2D piecewise constant splines regression approach for the image segmentation problem using a UNet [Ronneberger et al., 2015]. For differentiation, we use the formulation of splines discussed in Section 2. We analyze the efficacy of our approach by adding a piecewise constant DSA layer as the final layer of our network (MDSA). We compare this approach with the baseline model without the piecewise constant layer (Mbaseline).
We train two models (MDSA, Mbaseline) on two different segmentation tasks: the Weizmann horse dataset [Borenstein and Ullman, 2004] and the Broad Bioimage Benchmark Collection dataset [Ljosa et al., 2012] (publicly available under Creative Commons License). We split both the Weizmann horse and Broad Bioimage Benchmark Collection datasets into train and test with 85% and 15% of the dataset. We use binary cross-entropy error between the ground truth and the predicted segmentation map. We use the same architecture and hyper-parameters for both models (see Appendix for details.)
We observe that our DSA layer provides more consistent segmentation maps and higher Jaccard scores than the baseline model; see Figure 1. For the Weizmann horse dataset, Mconn enforces the connectivity of the segmented objects while also limiting noise in the segmentation map. In the cell
segmentation task, we note that the number of segments is high while the objects are small. Since the size of the components is small, our DSA layer Jacobian exhibits substantial differences from the commensurate identity gradient for the baseline models. Table 2 also shows the further improvement in Jaccard score on cell segmentation tasks over the Weizmann horse dataset.
3D point cloud reconstruction using NURBS: Next, we provide results for two experiments using DSA with NURBS discussed in Section 2.2. The first application is surface fitting for a complex benchmark surface represented by a mesh of surface points obtained by evaluating the benchmark test function at these points. We use Bukin function N.6 (publicly available here) for generating a grid of 256 ⇥ 256 points as shown on the left of Figure 2. For fitting a NURBS surface from the defined target point cloud, we initialize a uniform clamped knot vector for a cubic basis function and random control points of size 8⇥ 8. Using DSA, we evaluate the NURBS surface for a uniform grid of 256⇥ 256 parametric points. We now evaluate the surface and use mean squared error for fitting the surface point cloud using NURBS. We consider two scenarios: (i) we do not update the knot vectors (i.e., no reparameterization), and (ii) we compute the gradients for the knot vectors and allow for reparameterization (i.e., change of knot locations). We provide the comparison of these scenarios in Table 3. We see that the reparameterization helps in reducing the error in fit by half. Also, we notice that the density of points evaluated has a very minimal impact on the performance (see more details in Appendix). Visually, in Figure 2, we see that two knots in the “v” direction come close to each other around 0.06, enabling a sharp edge in the evaluated surface.
The next experiment we present involves surface reconstruction from point clouds using a graph convolutional neural network and DSA for unsupervised training. We use the SplineNet method proposed by Sharma et al. [2020] to be the baseline for point cloud reconstruction using splines. SplineNet uses a dynamic graph convolutional neural network (DGCNN) to predict the control points for a spline surface. The authors use a supervised control point loss to perform the training and include regularizations such as the Laplacian loss and a patch distance (using Chamfer distance) loss. Instead, we perform this training in an unsupervised manner by not using the control points prediction loss and only using DSA to evaluate the surface and then apply regularization of minimizing the Laplacian of the surface. Since we can train this unsupervised, we can even use an arbitrary number of control points and are not restricted to the target control points.
For a fair comparison, we use the same network, dataset, and hyperparameters as Sharma et al. [2020] and change the loss functions by removing the control point regression loss. For comparison, we compute the chamfer distance between the input point cloud and the NURBS surface fit by the DGCNN model (MDSA) (see Appendix for details of training). We use the Spline Dataset, which is
a subset of surfaces extracted from the ABC dataset (available for public use under this license). In Table 4, we provide a comparison of chamfer distance obtained between the predicted surface points from splines and the input point cloud for the test dataset. In our experiments, we observe that we get significantly better performance with fewer control points. This is because most of the surfaces in the dataset are simple curved surfaces that can be easily fit with fewer control points.
PDE based surrogate physics priors: Finally, we leverage DSA in the context of solving PDEs as a prior. In particular, we consider the Poisson equation solved for u:
r · (⌫(x)ru) = f(x) in D (12) u|@D = 0 (13)
where D = [0, 1]2, a 2D square domain, ⌫ is the diffusivity and f is the forcing function. We consider two experiments here: (1) validation of our approach with an analytically known solution, and (2) extending this to learn the solutions for the parametric Poisson equation parameterized using ⌫.
For the first experiment, we set ⌫ to 1 and the forcing f = f(x) = f(x, y) = 2⇡2 sin(⇡x) sin(⇡y), and minimize the residual using the approach described in Section 2.3. We know that for this PDE and the conditions provided, the exact solution is given by uex(x, y) = sin(⇡x) sin(⇡y). We compare our results (uDSA) with the exact solution uex. Also, we perform this experiment with Lagrange polynomials of different degrees. Further, we compare our results with results obtained using PINNs [Raissi et al., 2019]. We obtain significantly better performance (lesser `2-error by an order of magnitude) compared to PINNs, owing to more accurate gradients computed using our DSA approach. The performance improvement with increase in degree of polynomial in lower resolutions is more pronounced than at higher resolutions.
Next, we present results for training a deep learning network with a prior for solving a parametric Poisson’s equation. The input to the network are different diffusivity maps ⌫ sampled from
⌫(x;!) = exp
mX
i=1
!i i⇠i(x)⌘i(y)
! (14)
where !i is an m-dimensional parameter, is a vector of real numbers with monotonically decreasing values arranged in order; and ⇠ and ⌘ are functions of x and y respectively. We take m = 4, ! = [ 3, 3]4 and i = 1(1+0.25a2i ) , where a = (1.72, 4.05, 6.85, 9.82). Also ⇠i(x) = ai 2 cos(aix) + sin(aix) and ⌘(y) = ai2 cos(aiy) + sin(aiy). We generate several diffusivity maps by sampling this function with different values of !.
We use a UNet [Ronneberger et al., 2015] that takes these diffusivity maps and predicts the solution u, which is further optimized with the residual minimizing prior to the Poisson’s equation. Thus, we obtain a trained neural network that predicts the solution field u for any unknown diffusivity maps from the data distribution. We provide the predicted result along with its comparison with traditional numerical FEM results in Figure 3. Visually, we see both the predicted solution field map (uDSA) and the actual solution field (uground) obtained using traditional numerical methods match each other. The right most image shows the difference between both with the maximum deviation to be 0.01, showing the accuracy of our (easy-to-implement) DSA-based FEM solver.

heading: 4 Broader Impact and Discussion
text: We introduce a principled approach to estimate gradients for spline approximations. Specifically, we derive the (weak) Jacobian in the form of a block-sparse matrix based on the partitions generated by any spline approximation algorithm (which serves as the forward pass). The block structure allows for fast computation of the backward pass, thereby extending the application of differentiable programs (such as deep neural networks) to tasks involving splines. Our methods show superior performance than the state-of-the-art curve fitting methods by reducing the chamfer distance by an order of magnitude and the mean squared error in the case of surface fitting by a factor of two. Further, with the application of our methods in finite element analysis, we show significantly better performance than state-of-the-art physics-informed neural networks.
Our method is quite generic and may impact applications such as computer graphics, physics simulations, and engineering design. Care should be taken to ensure that these applications are deployed responsibly. Future works include further algorithmic understanding of the inductive bias encoded by DSA layers and dealing with splines having a dynamically chosen number of parameters (control points and knots).

heading: Acknowledgements
text: This work was supported in part by the National Science Foundation under grants CCF-2005804, LEAP-HI:2053760, CMMI:1644441, CPS-FRONTIER:1954556, USDA-NIFA:2021-67021-35329 and ARPA-E DIFFERENTIATE:DE-AR0001215. Any information provided and opinions expressed in this material are those of the author(s) and do not necessarily reflect the views of, nor any endorsements by, the funding agencies.

Review summaries:
The paper presents sparse Jacobians/backward algorithms for several spline formulations and motivates their use in differentiable algorithms. The formulations include general splines, NURBS, and basis functions for FEM solvers.  The resulting differentiable modules are applied to the three proof of concept tasks of image segmentation, point cloud reconstruction, and solving Poisson equations, where they show promising results.

The paper presented a differentiable approximation to splines for use as a differentiable layer in differential programming. The main idea is to exploit the sparsity of the resulting Jacobian matrix for the spline approximation operation, and to use a Guassian approximation for the non-differentiable Direc delta spike functions from the Cox-de Boor base function. The authors used this differentiable spline approximation in three application settings: 2d image segmentation; 3d surface fitting for point clouds using NURBS, and solving a PDE (Poisson Equation).