Differentiable Spline Approximations
Authors: Minsu Cho
Aditya Balu
Ameya Joshi
Anjana Deva Prasad
Biswajit Khara
Soumik Sarkar
Baskar Ganapathysubramanian
Adarsh Krishnamurthy
Chinmay Hegde
Venue: NeurIPS 2021 Poster (2021)
Keywords: Differentiable programming
Spline approximation
NURBS
k-histogram

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 that the machine learning models 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 \url{https://github.com/idealab-isu/DSA}.

Differentiable Spline Approximations
Minsu Cholt
Aditya Balu
Ameya Joshi!
Anjana Deva Prasad²
Biswajit Khara²
Soumik Sarkar2
Baskar Ganapathyubrammaninn
Adarsh Krishnamurthy
Chinmay Hegde'
New York University!, Iowa State University² fmc8065, ameya.joshi chinmay. h) edu
{baditya, anjana, bkhara, soumiks, baskarg, adarsh)eiastate. edu
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 applicabílity. Our goal in this paper is use new, principled approach extend gradient-bass optimization to functions well modeled by splines, which encompass 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 diverse applications such as image segmentation, 3D point cloud reconstruction, and finite element analysis. We also open-source the code at https: //githbb..ttttt/// com idealab isu/DSA.
Introduction
Differentiable programming has been paradigm shift in algorithm design. The main idea to leverage gradient- based optimization to optimize the parameters 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., 2018 8a], physics engines [Degrave et al., 2017], computational simulations [Alnæs et al., 2015], and graphics [Li et al., 201 8b; 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 monotonicity, or piecewise constancy, are particularly prevalent in upplications 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.
Motivation:
For differentiable programming to work, all layers within the model must admit simple gradient calculations; however, this poses major limitation in many settings. For example, consider computer graphics applications for rendering 3D objects [Kindlmann al., 2003; Gross et al., 1995; Loop and Blinn, 2006]. A common primitive in such cases is spline (or piecewise polynomial) function which either exactly or approximately interpolates between discrete set of points to produce continuous shape surface. Similar spline (or other piecewise polynomial) approximations arise in
Challenges:
* Equal contribution.
35th Conference on Neural Information Processing Systems (NcurIPS 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 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 al., 2019; Deng et al., 2020].
Our contributions: We propose a principled approach for differentiable programming for spline functions without the use of continuous relaxation 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 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-lineal 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 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 al., 2018]. However, introducing more complex non-differentiable components requires careful treatment due to undefined badly behaved gradients. For example, in the case of sorting and ranking operators, 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. Berthe 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 [Pogančić 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.
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
Spline approximation:
While 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 neerrllnewwork-aae contour error prediction method for NURBS paths. Fey al. [2018] present 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 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., 2020J. Most often, the physics phenomena are represented using partial differential equations (PDEs) SSanchez-Gonzale et al., 2020; Holl 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 detailed survey of this (vast) area by Cai et al. [2021].
Differentiable Spline Approximation
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 casypto-comput 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 E T be a vector where the 2th element is denoted as fi. Let us use [m] = =1,,2.... n.) to denote the set of all coordinate indices. For a vector f E I" and an index set I C [n]. let fr be the restriction of f to I, i... for E I, we have fi) := fi. and fr(
0 for i & I. Now, consider any fixed partition of [ni] into set of disjoint intervals = {I...... J where the number of intervals |Z| = k. The f--norm of f is written as || == v E F? while the & distance between f,g g is written as ||ʃ - 9912-
We first define the notion of 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 special instance of this definition; we explain further below.)
Def. 2.1 (Discretized k-spline). A vector h R is called discretized k-spline with degree d if: (i) there exists a partition of [n] into is disjoint intervals I1,... I:: (ii) within each interval I,, the coefficients of hj. j E I, can be perfectly interpolated by some polynomial function of degree d.
Let us illustrate this by an example. Suppose that d = 1 and A = 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,hy) hy) for all j E [rz]). 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- - IR with A = 5 pieces (with suitably defined knot points, which are the location of the intervals 2), and evaluate / at any n equally spaced points in its domain. This gives us a vector h E R", 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 continuity or differentiability properties. Therefore, our definition encompasses all standard spline families (including intrpolaiig/appoximating splines such as smoothing- cubic-, and B B-splines).
specific
2.1 Spline Approximation
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 given vector space IR.,, consider Sh, 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 [r]] into k subsets, the family of discretized k-splines is a k(d + 1) ---dimensional subspace of IR.. Now suppose that the knot indices are allowed to vary. The number of possible partitions is finite (of the order of (*)), and therefore the set St is a finite union of subspaces, or a nonlinear submanifold, embedded in IRz..
Therefore, the problem of discretized k-spline approximation can be viewed as an orthogonal projection onto this nonlinear manifold. Consider any arbitrary vector z E IR" (we can think of (i, 2i)
3

---

as set of n data points to which we are trying to fit k-spline). Then, the best k-spline fit to æ (in the sense Eg distance) amounts solving the optimization problem:
F(x)= arg min
2 ||
1111
(x)
h2)² s.t. h
(I)
h
1=1
This operation resembles standard spline regression. But it strictly more general since this requires only optimizing piecewise spline parameters but also the knot indices. Crucially, note that F is both con-differentiable and non-convex map. Nevertheless, such an orthogonal projection can be computed in polynomial (in fact, nearly-linear) time [Jagadish al., 1998; Acharya et al., 2015] using many different techniques, including dynamic programming. This forms the forward pass of DSA "layer".
Our first main conceptual contribution 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-differentiabiiity built into the definition of the k- spline). Therefore, will instead he 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 corre- sponds to partition of [n]. that is, each element of h2 corresponds to a single interval, Lj. and (2) within each interval, the least-squares operation is continuous and differentiable. The first property ensures that every element 2; contributes to only a single piece in the output h. Given that the sub-functie fobfffnctions from the piecewise partitioning function smooth, also observe that the size of each block corresponds to the size of the partition, 4, Using this observation, we get:
Theorem 1.
The Jacobian of the operation F with respect to Rd can be expressed as block diagonal matrix, RRax,, whose (8, 1)th entry obeys:
(F(()))) @h(z), Sr,
Bhg (x),
if
otherwise
As concrete instantiation of this result, consider the case d 0. This is the case where wish to best approximate the entries of x with at most "horizontal" pieces, where the break-points are obtained during the forward pass Call this approximation h. Then, the Jacobian of h with respect to x forms the block-diagonal matrix E nxx.
JI
0
0
J2
0
J
JR.
where all entries of each block, J; TR. 4 are constant and equal to 1/[Ta], i.e., they are www/column-stochastiee 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 any degree 2 1, and generalization to 2D domains (surface approximation). In Section we demonstrate the utility of this approach for 2D segmentation (i.e., piecewise constant approximation) problem, similar to the setting studied in Djolonga and Krause [2017].
II the data summarization literature, this class of functions is sometimes called & k-histograms [Jagadish al., 1998]

---

2.2 Differentiable NURBS
We now switch to 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 continuous function C RR R defined as follows. Construct any knot vector u (i.e. non-decreasing sequence of real coordinate values) and fix degree d. Recursively define a sequence of basis functions, No R - IR computed using the Cox-de Boor formula:
(u)
24 u
No (w) + U1+d+1 (4), NV(()= {
uf S u S u4+1 0 otherwise
1
wi+d
14+d+1 U;+1
for d = 1, ..... In the uniform case (where the knots are equally spaced), each No can be viewed as being generated by recursively convolving box function with A The non-uniform case cannot be written as convolution, but the intuition is similar. With these basis functions in hand, the NURBS curve C is defined as the rational function:
Ei=o NT(u)),PP
C(u) = Σi= (1)
(5)
where P,, i == 0, 1..... t are called control points and uy 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 = + 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 called the partition of unity property). Therefore:
C(u)= = Σ (u)P
(6)
==0
In summary, the NURBS curve is parametrically defined via the control points and the knot positions. This discussion is for ID curves, but 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, Uy 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 u.). Once again, we see that the non- con-differentiabliity 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 N; (u) = (0,,0141) with respect to the knot points, which is not well defined. However, Nº (u) can be viewed as the difference between convolutions of the unit step function with on, and 6u.. where 8 is the Dirac delta defined over the real line. We smoothly approximate the delta function by Gaussian function with small enough bandwidth hyperparameter a: 8(4)) R g(u) = exp(-
u4)/20²). This function is now n² (()). 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 for pseudocode and the Appendix for details.
differentiable with respect to d; with g((u) 1-thh
Algorithm I Backward pass for NURBS Jacobian (for one curve point, C(u))
P', U': gradients of C' w.r.t. P.U Initialize: P', U' - 0
Retrieve Uapan Nd. C(u) calculated during forward pass /* Uapan is the index of knot position
/+ No is the basis function of degree d +/
/* C(u) is the evaluated curve point +/
for h 0: d do
P th No. // easy since C is linear function of P.
U th Ng. U Uhpanth // due to Gaussian approximation; discussion below.
5

---

Let us probe the structure of this Jacobian a bit further. Suppose we evaluate the curve Cat at n arbitrary domain points. There are slightly less than control points, and therefore the Jacobian is roughly of size x O(k). However, due to the recursive nature of the definition of basis functions, the spun (or support) of each basis function small and only touches d + 1 knots; for example, only 2 knots affect No. only 3 knots impact NJ and so on. This endows a natural sparse structure on the Jacobian. Moreover, for a fixed order parameter d + I, 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.
2.3 Differentiable Finite Element PDE Solvers
Next, we see how spline approximations can be used to improve finite element analysis for solv- ing 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 al., 2020]. However, these approaches come with challenges while used in conjunction with autodiff packages, and gradient pathologies pose major barrier (Wang et al., 2020].
Using our principles developed above, we propose an alternative PDE solution approach via differen- tiable 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:
[[) F(w), U E S2,
where N denotes the differential operator and U S2 - R is a continuous field variable; it is common specify additional boundary constraints on U. The Galerkin method converts solving for the best possible U (which is continuous variable) into a discrete problem by first looking at the weak form: R(U) fo V (UU 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. typical set of basis functions L is obtained by (piecewise) concatenation of polynomials, each defined over elements of a given partition of Q (also called mesh). Commonly used choices include Lagrange polynomials, defined by:
(
ΣV. U, II
2 My
Um s.t. Ty. E [-1,1]
Um
where {uo. M1..... Fae are finite set of nodes (akin to control points in our above discussion, except in this case the splines interpolate the control points) and U, is the corresponding coefficient. We use this collection of basis functions D to represent U:
#nodes
U
Σ AW
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 U = } into the definition of R. we get standard Finite Element form,
R(US, V9) = B(U V VVV-)(VV)
(10)
where B(U Vo) is the discrete form (bilinear for linear operators) that encodes the differential operator and L(U) is 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:
((U9) 2 B(U U9) - L(U)).
(11)
Optimization of this energy functional can now be performed using gradient- based iterations evaluated by automatic differentiation. This is powerful approach since formal techniques exist (e.g., Galerkin Least Squares Bochev and Gunzburger [2009]) that reformulate weak forms PDEs into equivalent energy functionals. The key aspect to note here is that differentiating "through" the differential operator A (embedded within B) requires derivative computations of the piecewise polynomial basis functions sss, and therefore our techniques developed above are applicable.

---

3 Experiments
We have implemented the DSA framework (and its different applications provided below) by ex- tending autograd functions in Pytorch. We also provide the capability run the code using CUDA for GPU support. All the experiments were performed using a local cluster with 6 com- pute nodes and each node having 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 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:// thub. com/ideale aideallbbbiie isuDDA..
Image segmentation: We begin with implementing 2D piecewise constant splines regression approach for the image segmentation problem using a UNet [Ronneberger et al., 2015]. For differ- entiation, 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 (Maaseline)
We train two models (MDSA Museliie) 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 syyee--arameters for both models (see Appendix for details.)
We observe that our DSA layer provides more consistent segmentation maps and higher Jaccard scores than the bascline model; Figure 1. For the Weizmann horse dataset, Mceun enforces the connectivity of the segmented objects while also limiting noise in the segmentation map. In the cell
Table 1: Summary of Experiments: We present three experiments in this paper with diverse applications, model architectures, basis functions, and degree of splines.
Application
Input
Output
Architecture
Spline Function
Degree
Image
2D piecewise constant
Box-Car
Image
U-Net
Segmentation
knot partitions
functions
Point cloud
3D control points
Dynamic Graph
BSpline
3D point cloud
reconstruction
and rational weights
CNN
polynomials
PDE-based
2D Mesh
2D physics field
Lagrange
U-Net
1,2,3
surrogate physics
(regular grid)
on the mesh grid
polynomials
Original
Ground
Mhaseline
MDSA
Original
Ground
M baseline
MDSA
Figure 1: Segmentation results. The two models, MDSA and Mbassiiee were trained with and without the DSA layer, respectively. Note that MDSA generates better segmentation masks with fewer holes and enforced connectivity. Note the sharper edges compared to the standard segmentation results. Additional figures are the Appendix.
7

---

Table 2: Results for the horse and cell segmentation dataset: Jaccard scores for the baseline and connected component models for the cell and horse segmentation task. From independent three runs with random seeds and the table reports mean and standard deviation. As the objects of interest (piecewise constant components) are smaller, the model with the DSA layer learns better representation. Predictions are thresholded at 0.5.
Dataset
Baseline (Masseine))
Baseline DSA (MDsa)
Weizmann Horse [Borenstein and Ullman] Broad Bioimage Benchmark [Ljosa et al.]
72.06 0.60
73.13 0.31
79.34 + 0.43
81.56 + 0.24
segmentation task, note that the number of segments is high while the objects small. Since the size of the components 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 complex benchmark surface represented by 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 x 256 points as shown on the left of Figure 2. For fitting 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 x 8. Using DSA, evaluate the NURBS surface for uniform grid of 256 x 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) do not update the knot vectors (i.e., no reparameterization) and (ii) we compute the gradients for the knot vectors and allow for reparameterizatioo (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 very minimal impact 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 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 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 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
0.005
0005
3.25
4,005
500
0J5
0.50
can
Predicted (no reparameterizstonn)
Pointwise Normalized
Predicted NURBS (with reparametertzation))
Normalized
MSE
MSE
Figure 2: NURBS surface fitting results: Surface fitting to point cloud generated using the Bukin's 's function N.6 given by = 100 0.01I2 + 0.01| + 10); -15 < x < -5, -3 < y 3. The center image shows the surface fit obtained without reparameterization of the knots. We obtain better fit by reparameterizing the knots.
8

---

Table 3: NURBS surface fitting results: Comparison of mean squared error between the target surface point cloud and the surface generated using DSA with and without reparaaettrizatton.
Number of Points
(wihhoit reparameterizzation) (wth reparamtterizzaionn
128 128
19.83 + 0.001
8.25 + 0.01
256 256
19.85 + 0.001
8.23 + 0.02
Table 4: Point-cloud reconstruction results: Comparison between the model proposed by Sharma et al. [2020] and its extension using DSA (with different number of control points), We compare the two-sided chamfer distance (scaled by 100) between the input point cloud and the fitted surface.
Mhaaeline
MDSA
MDSA
MDSA
Experiment
(20 20)
(20 x 20)
(5 x 5)
(4 x 4)
Chamfer Distance
1.18 0.10
0.03 + 0.02 0.14 + 0.07 0.02
0.01
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, significantly better performance with fewer control points. This because most of the surfaces in the dataset are simple curved surfaces that can be easily fit with fewer control points.
observe that get
PDE based surrogate physics priors: Finally, we leverage DSA in the context of solving PDEs as prior. In particular, we consider the Poisson equation solved for u:
(v(x) vu) fxx in D
(12)
ulaD = 0
(13)
where D = [0, 1j2, a 2D square domain, v is the diffusivity and I is the forcing function. We consider two experiments here: (1) validation of our approach with an analytically known solution, and (2) extending this learn the solutions for the parametric Poisson equation parameterized using v.
For the first experiment, we set v to 1 and the forcing f = f(z) = f(x y) = 2n²² sin(xx)) sin(ry). 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 uez (x, y) = sin(xx) sin(ry). We compare our results (uDSA) with the exact solution Mez- Also, 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.
perform this experiment with
Next, present results for training a deep learning network with prior for solving parametric Poisson's 's equation. The input to the network are different diffusivity maps v sampled from
v(x:w) = exp Σ (.AAE(())))())
(14)
where W; is an m-dimensional parameter, À is vector of real numbers with monotonically decreasing values arranged in order; and & and 1/ are functions of a and respectively. We take m 4, (1.72, 4.05, 6.85, 9.82). Also E,(x) 3/2 cos(a;2)
- 1-3, 3) and X; (I+0.2507) where sin (a,ar) and n73) = et cas(a,y) + sin(a.v). We generate several diffusivity maps by sampling this function with different values of w.
Table 5: Quantitative comparison of Solving PDEs: T.2 Norm between the analytical exact solution Vez and predicted u using PINNs [Raissi et al., 2019] and DSA with different degrees of the Lagrange polynomials.
Model
PINN
DSA (d = 1)
DSA (d = 2)
DSS (d 3)
128 x 128
3.72 + 0.20 E-4
3.32 + 0.05 E-5
2.16 + 0.04
2.37 + 0.10 E-5
Lg Norm
256 256
2.63 + 0.20 E-4
2.57 + 0.01 E-5
2.79 + 0.20 E-5
2.59 + 0.10 F-5

---

1.0
3.0
0.01
0.0
0.0
-0.01
0.0
0.0
V
1.0
UDSA
grround
UpSA -Uground
Figure 3: Learning a parametric faailly of PDE solutions: Poisson's 's equation with log permeabil- ity coefficients w (-0.26, 0.77, -0.37, -0.92) the Poisson's equation.
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 '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 (uDS and the actual solution field (Ugrownd)) 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.
4 Broader Impact and Discussion
We introduce principled approach to estimate gradients for spline approximations. Specifically, we derive the (weak) Jacobian in the form of 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 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 factor of two. Further, with the application of our methods finite element analysis, we show significantly better performance than state-of-the-an physies-ine physics-informed informed neural networks.
Our method is quite generic and may impact applications such as computer graphics, physics simulations, and engincering 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).
Acknowledgements
This work was supported in part by the National Science Foundation under grants CCF-2005804. HAAP-HI:0533600 CMMI: 1644441, CPS- POS-FRRNTIER: 1954556, USDAA-NIFA20221-6722--555 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.
and ARPA-E
10

---

References
Michael Innes. Algorithmic differentiation. In Machine Learning and Systems, pages 1-12, 2020.
Mike Innes, A. Edelman, K. Fischer, C. Rackauckas, E. Saba, V. B. Shah, and Will Tebbutt. A differentiable programming system to bridge machine learning and scientific computing. ArXiv, abs/ 907.07587, 2019.
Schafer, M. Kloc, C. Bruder, Lorch. A differentiable programming method for quantum control, ArXiv, 2020.
Tzu-Mao Li, Michaël Gharbi, Andrew Adams, Frédo Durand, and Jonathan Ragan-Kelley. Differentiable programming for image processing and deep learning in halide. ACM Transactions on Graphics, 37(4))1-13, 2018a.
J. Degrave, Michiel Hermans, J. Dambre, and Wyffels. A differentiable physics engine for deep learning in robotics. Frontiers Neurorobotics, 13, 2017.
Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Maric Rognes, and Garth Wells. The FEniCS project version 1.5. Archive of Numerical Software, 3(100), 2015.
Tzu-Mao Li, Miika Aittala, Frédo Durand, and Jaakko Lehtinen. Differentiable Monte-Cart ray tracing through edge sampling. ACM Transactions on Graphics, 37(6):1-11, 2010)) 2018b.
Wenzheng Chen, Huan Ling, Jun Gao, Edward Smith, Jaakko Lehtinen, Alec Jacobson, and Sanja Fidler. Learning to predict 3D objects with an interpolation- based differentiable renderer. In Advances Neural Information Processing Systems, volume 32, pages 1-11, 2019.
Gordon Kindlmann, Ross Whitaker, Tolga Tasdizen, and Torsten Moller. Curvature-hased transfer functions for direct volume rendering: Methods and applications. In IEEE Visualization, 2003. VIS 2003., pages 513-520. IEEE, 2003.
Markus H Gross, Lars Lippert, A Dreger, and R Koch. A new method approximate the volume-rendering equation using wavelet bases and piecewise polynomials. Computers Graphics, 19(1):77-2,, 1995.
Charles Loop and Jim Blinn. Real-time GPU rendering of piecewise algehraic surfaces. In SIGGRAPH, pages 6644670. ACM, 2006.
Thomas JR Hughes, John Cottrell, and Yuri Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comp. methods in Applied Mechanics and Engineering, 194(39-41): 4135-4195, 2005.
Anantharam Balakrishnan and Stephen C Graves. A composite algorithm for a concave-cost network flow problem. Networks, 19(2):175-203: 1989.
Mathieu Blondel, O. Teboul, Quentin Berthet, and Josip Djolonga. Fast differentiable sorting and ranking. ArXiv, abs/2002.08881,, 2020.
Marco Cuturi, Teboul, and Phillipe Vert. Differentiable ranking and sorting using optimal transport. In Neural Information Processing Systems, 2019.
A. Mensch and Mathieu Blondel. Differentiable dynamic programming for structured prediction and attention. ArXiv, 0011880006777 2018.
Josip Djolonga and Andreas Krause. Differentiable learning of submodular models. In Adv. Neural Inf. Proc. Sys. (NeurIPS), pages 1013-10233 1023-1022,1110 2017.
A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and Z. Kolter. Differentiable convex optimization layers. In Adv. Neural Inf. Proc. Sys. (NeuriPS), 2019.
Boyang Deng, Kyle Genova, Soroosh Yazdani, Sofien Bouaziz, Geoffrey Hinton, and Andrea Tagliasacchi. CvxNet: Learnable convex decomposition. In IEEE Conf. Comp. Vision and Pattern Recog. JEBE, 2020.
Yoshua Bengio. Estimating or propagating gradients through stochastic neurons. ArXiv, abs/ 30..000..22... 2013.
Penghang Yin, Jiancheng Lyu, Shuai Zhang. Stanley Osher, Yingyong Qi. and Jack Xin. Understanding straight- through estimator in training activation quantized neural nets. In Proc. Int, Conf. Learning Representations (ICLR), 2019.

---

Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: survey. J. Machine Learning Research, 18(153):1-43, 2018. URL .tmml
http:// /ttp://jmlr. wtgppaperssv18/17--11
Quentin Berthet, Mathieu Blondel, O. Teboul, Marco Cuturi, Phillipp Vert, and Francis Bach. Learning with differentiable perturbed optimizers. ArXiv, ahs/2002.08676, 2020.
Yujia Xie, Hanjun Dai, M. Chen, Bo Dai, Tuo Zhao, H. Zha, Wei Wei, and T. Pfister. Differentiable top-k operator with optimal transport. ArXiv, abs/2002.06504, 2020.
Hyunsung Lee, Yeongjae Jang. Jackwang Kim, and Honguk Woo. A differentiable ranking metric using relaxed sorting opeartion for top-k recommender systems. ArXiv, abs/2008. 13141, 2020.
Marin Vlastelica Pogančić, Anselm Paulus, Vit Musil, Georg Martius, and Michal Rolinck. Differentiation of blackbox combinatorial solvers. In Proc. Int. Conf. Learning Representations (ICLR), 2020. URL
Brandon Amos and J. Zico Kolter. Optnet: Differentiable optimization layer in neural networks. ArXiv, 1703.00443, 2019.
Sheroze Sheriffdeen, J. Ragusa, Morel, M. Adams, and T. Bui-Thanh. Accelerating PDE-constrained inverse solutions with deep learning and reduced order models. ArXiv, abs/1912.08864. 2019.
Ameya Joshi, Minsu Cho, Viraj Shah, B. Pokuri, Soumik Sarkar, Baskar Ganapathysubramanian, and Chinmay Hegde. InvNet: Encoding geometric and statistical invariances in deep generative models. In Association for the Advancement of Artificial Intelligence Conference, pages 1-8, 2020.
Les Piegl and Wayne Tiller. The NURBS Book (2nd Ed.). Springer- Verlag, Berlin, Heidelberg, 1997. ISBN 3540615458.
Ludovico Minto, Pictro Zanuttigh, and Giampaolo Pagnutti. Deep learning for 3D shape classification based on volumetric density and surface approximation clues. In VISIGRAPP (5: VISAPP), pages 317-324, 2018.
Krystian Erwinski. Marcin Paprocki, Andrzej Wawrzak, and Lech M Grzesiak. Neural network contour error predictor CNC control systems. In 2016 21 International Conference on Methods and Models in Automation and Robotics (MMAR), pages 537.542. IEEE, 2016.
Matthias Fcy, Jan Eric Lenssen, Frank Weichert, and Heinrich Müller. SplineCNN: Fast geometric deep learning with continuous B-spline kernels. In Conference on Computer Vision and Puttern Recognition, pages 869-877, 2018.
Gopal Sharma, Difan Liu, Subhransu Maji. Evangelos Kalogerakis, Siddhartha Chaudhuri, and Radomir Mčch. ParSeNet: A parametric surface fitting network for 3D point clouds, 2020.
Yuanming Hu, Luke Anderson, Tzu-Mao Li, Qi Sun, Nathan Carr, Jonathan Ragan- Kelley, and Frédo Durand. Difftaichi: Differentiable programming for physical simulation. arXiv preprint 1110v0991,.00935, 2019.
Yi-Ling Qiao, Junbang Liang, Vladlen Koltun, and Ming C Lin. Scalable differentiable physies for learning and control. arXiv preprint arXiv:2007.02/12 02.666....8........ 2020.
Alvaro Gonzalzz, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International Conference on Machine Learning, pages 8459-84440 8468. PMLR, 2020.
Philipp Holl, Vladlen Koltun, and Nils Thuerey. Learning to control pdes with diffcrentiable physics. arXiv preprint arXiv:2001.0744 2020.
M. Raissi, P. Perdikaris, and G.E. Kariiadakis. Physics-informen neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686- 707, 2019. ISSN 0021 doi:
http:///g///..1111.....00001/1 URL http: //tp://wwwww//ao wciencedirectt.....cm... com/ http///ccinneearrree /eciiine//stttcce///////////////
Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning of nonlinear partial
differential equations. Journal of Computational Physics, 357:125-1411 2018.
Ehsan Kharazmi, Zhongqiang Zhang, and George Em Karniadakis. hp- VPINNS: Variational physics physics-infommed informed neural networks with domain decomposition. Computer Methods in Applied Mechanics and Engineering 374: 13347, 2021.
12

---

Shengze Cai. Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis. Plysics-informed neural networks (pinns) for fluid mechanics: A review. arXiv preprint arXiv:2105.09506, 2021.
H. V. Jagadish, Nick Koudas, S. Muthukrishnan, Viswanath Poosala, Kenneth C. Sevcik, and Torsten Suel. Optimal histograms with quality guarantees. In Proc. Int. Conference Very Large Data Bases (VLDB), 1998.
Jayadev Acharya, Ilias Diakonikolas, Chinmay Hegde, Jerry Zheng Li, and Ludwig Schmidt. Fast and near- optimal algorithms for approximating distributions by histograms. In Proc. ACM SIGMOD- SIGACC-IGGTTAAIIA Symp. on Principles of Database Systems, 2015.
Sifan Wang. Yujun Teng, Paris Perdikaris. Understanding and mitigating gradient pathologies in physics- informed neural networks, arXiv preprint urXiv:2001.04554 2020.
Pavel Bochev and Max D Gunzburger. Least-squares finite element methods, volume 166. Springer Science Business Media, 2009.
Olaf Ronneberger. Philipp Fischer, Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computer-asiisted intervention, pages 1.4-241. Springer, 2015.
Eran Borenstein and Shimon Ullman. Learning to segment. In Euro. Conf. Comp. Vision, pages 315-328. Springer. 2004.
Vebjom Ljosa, Katherine L. Sokolnicki, and Anne E Carpenter. Annotated high-throughput microscopy image sets for validation. Nature methods, 9(7):637-66 637, 2012.
Sudipto Guha, Nick Koudas, and Kyuseok Shim. Approximation and streaming algorithms for histogram construction problems. ACM Truns. Dutabase Systems (TODS), 31(1))396-438) 2006.
Kesheng Wu, Ekow Oton, and Arie Shoshani. Optimizing connected component labeling algorithms. In Medical Imaging 2005: Image Processing, volume 5747, pages 1965-1976, International Society for Optics and Photonics, 2005.
Adarsh Krishnamurthy, Rahul Khardekar, Sara McMains, Kirk Haller, and Gershon Elber. Performing efficient NURBS modeling operations the GPU. JEEE Transactions on Visualization and Computer Graphics, 15 (4):530-544 543, 2009.
13


Reviews:
[{"review_id": "bcMArkZsP3w", "reviewer": "Reviewer_ZLi8", "summary": "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).", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is rather math dense and it took me some time to understand the main proposal set forth by this work. The main pros and cons are concluded below.\n\nPros:\n+ I think the differentiable approximation of the spline functions is a very interesting topic, potentially allowing the integration of spline function approximations into the differentiable programming paradigm. The current field of deep learning is heavily dominated by \"neural\" (i.e., Multi-Layer Perception) approximations for parameterizing high dimensional functions, which are flexible and extensible, but has its own set of limitations. I think generalizing the idea of differentiable splines as function approximations is attractive.\n+ The idea of smoothing the discontinuous delta functions using Gaussian with small bandwidth feels intuitive.\n+ I like the error intervals reported in the experiments.\n\nCons:\n- Though the paper is quite content heavy, most of the stuff is not new (e.g., when computing the Jacobian, only the gradients corresponding to the specific segments are non-zero). The main contribution seems to be that the authors proposed to smooth the step function by approximating it with a Gaussian kernel. However, this main aspect does not seem very thoroughly studied. \n    - For example, how does the bandwidth of the Gaussian affect the results (with empirical evidence)? I would assume a bandwidth too small will result in nasty gradients, and a bandwidth too large will produce inaccurate approximations. \n    - What if you replace the step function in the base case (Eqn. 4) with two sigmoids of different bandwidth instead of using a Gaussian? What about using two tanh? Different options can be explored in greater depth.\n- The experimental section is not very well executed.\n    -  Following the point above, ablation studies for the choice of smoothing function are lacking.\n    - The experiments lack details and context. For example, I find it difficult to understand how the differentiable spline layer is used in the image segmentation experiments. The only description seems to be \"We begin with implementing a 2D piecewise constant splines regression approach for the image segmentation problem using a UNet\" (l 242). What does that mean? Are the shape contours represented using NURBS? How is the output parameterized? What parameters are the neural network outputting after the UNet feature map? The baseline method is simply described as \"without the DSA layer\", with not enough details about the DSA layer for the reader to understand. The evaluation is a little hard to grasp. How is the Jaccard score defined? Is it the same as pixel-wise IoU score? For the PDE experiment, is my understanding correct that the first part looks at directly optimizing the Lagrange polynomial based representation (by minimizing PDE residuals), and the second part looks at training for a learned prior over a dataset (using a UNet)? Similar question to the first experiment, how are the U-Net feature maps associated with the spline polynomials?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "In general I think differentiable spline is an interesting topic. A few areas of improvement I would love to see:\n1. More comprehensive ablations on the smoothing part for differentiability.\n2. Clear details for the experiments. If space is an issue, consider moving more details to the suppmat.\n3. Some performance benchmarking would be nice-to-have. How does the algorithm performance scale with number of knots / degree? How fast is this layer compared with other NN layers?", "review_text": "", "extra_scores": {"ethical_concerns": "n/a", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3 hours"}, "extra_text": {}}, {"review_id": "Rq9P_znzUPA", "reviewer": "Reviewer_Q18T", "summary": "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. \nThe 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.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality:\n- To my knowledge, this is the first method that allows backpropagating gradients to the knot vectors of NURBS curves. The authors achieve this by framing the box function as convolutions with Dirac impulses and approximating those Dirac impulses using Gaussian functions. This seems to be an interesting idea, which possibly could be applied to other situations as well.\n- The observation that the same principles can be applied to backpropagate through functions in FEM solvers is interesting and has potential. There might be several useful applications.\n- Utilizing the sparsity of the Jacobian (induced by the local support property) to develop efficient backward algorithms through splines was already done by previous work (e.g. SplineCNN [1]). However, this is the first work that presents that concept in a principled way for general splines.\n\nQuality:\n- The paper is technically sound\n- The experiments show that the method is indeed applicable in several domains.\n- All three experiments are proofs of concept. I wonder what hinders the application to more sophisticated tasks. At least for image segmentation it should be possible to compare to SOTA methods on common segmentation benchmarks. Instead, the authors only compare to their own baseline.\n\nClarity:\n- The paper is well-written, well structured and nice to read.\n- The paper would profit from some explanatory figures for the general splines in Section 2/2.1 and the FEM method in Section 2.3.\n\nSignificance:\n- This is a paper which can spark a lot of applications in several areas, if the presented methods turn out to be practical. There are several interesting ideas and concepts.\n\nAll in all, I recommend to accept the paper for the mentioned reasons. The experiments certainly leave room for improvement, mostly because the authors are not comparing to other SOTA methods. However, the given experiments reach the bar in that they showcase the different potential applications as proofs of concept.\n\nQuestions:\n- Is the approximation of Basis functions (using Gaussians instead of Dirac impulses) only done in the backward step or also in the forward step? If it is done in both I wonder if the approach would also work if the normal basis is used in the forward computation, while the \"surrogate\" is used for gradient computation. \n- How efficient is the FEM solver? Would it be possible to apply it to larger scale tasks, for example volumetric reconstruction tasks as in SIREN [2]?\n- In the given applications, the knot vectors (and therefore the gradients computed for them) lie in 2- or 3-dimensional space. When performing gradient descent in such low dimensional spaces, one often runs into bad local minima. Did the authors face any of such challenges when optimizing their positions?\n\n\n[1] Fey et al.: SplineCNN: Fast Geometric Deep Learning with Continuous B-Spline Kernels, CVPR 2018\n\n[2] Sitzmann et al.: Implicit Neural Representations with Periodic Activation Functions, NeurIPS 2020\n\n----------------- After Rebuttal ---------------\n\nI thank the authors for answering my questions. I will keep my score.\n\nThe other reviewers correctly note that in this paper some contributions are oversold, as both ideas, utilizing sparsity of the Jacobian for efficient backprop and computing gradients through Dirac functions via smooth approximations are not novel but already present in previous work. I agree and encourage the authors to adjust the statements in future versions. \n\nHowever, all in all, I see a lot of value in this work, as it is the first that builds a full differentiable building block for general splines, which can be useful to many researchers.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors could further discuss the limitations of their method. This would help to assess the applicability of the presented approach for real-world tasks. I see no immediate societal impact.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "ApvOwLQQz3d", "reviewer": "Reviewer_qxvV", "summary": "The authors propose to derive a technique allowing gradient based learning of splines (or piecewise polynomials) and to apply such a method to standard tasks such as PDE solving/approximation and other function approximation problems.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Some notations are counter intuitive e.g. f_I(i) since F is a vector and F_I is a restriction (in term of dimension) of this version, it seems unnatural to call it as a function (especially since it is not well introduced/defined).\n\nThe problem of both optimizing the spline (per-region) function parameters and the knots of the partition is well studied yet no reference is provided (see all the literature on adaptive spline, and many DeVore references).\n\nTheorem1 is merely a direct calculus result that is straightforward and has been known in the spline community for a while (it is trivial to show the block structure based on the fact that each subset of parameters is only ‘’active’’ within its partition region)\n\nMany english formulation that are not appropriate ‘’the backward pass is a bit more tricky’’ and so on. Additionally the whole organization and choice of organization is odd (e.g. bold paragraph in the intro and so on).\n\nMany sentences seem incomplete e.g. “This is powerful since ….into equivalent energy functionals”.\n\nNo comparison is done against alternative spline fitting methods including the ones that allow adaptivity of the partition though mesh refinement (which is the whole direction of PDEs where mesh refinement is able to somewhat compete with adaptive knot splines).\n\n\n----- AFTER REBUTTAL ------\n\nI thank the authors for their detailed answers and acknowledging the drawbacks I have raised. I am raising my score from 4 to 5 for the following reason. In my opinion, some of the technical contributions that are the pillars of the paper have been obtained before (e.g. using Gaussian (or more generally narrow support function) approximation for gradient computation in https://arxiv.org/pdf/1706.04698.pdf https://arxiv.org/pdf/1901.09948.pdf, gradient sparsity which is (as also pointed out by ZLi8) known and already leveraged in various fields such as PDE approximations (the sparsity of the gradient is what allows most methods to actually be tractable on a fine mesh)). And when considering the remaining contributions, I am uncertain that this paper passes the acceptance threshold. However, some reviewers are pointing out that this paper offers an interesting application and find the submission useful, hence I am raising my score accordingly.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Congratulations, the paper is accepted to NeurIPS 2021!
Please reformulate your contributions in light of reviewer's qxvV comments. 
Please incorporate other corrections and additions from the rebuttal/reviews.
Decision: Accept (Poster)