1 Introduction↩︎

The multivariate nonparametric regression problem is a fundamental statistical primitive, with a vast history and many approaches. We adopt the following setup: given a dataset, \(\{(\boldsymbol{x}_i, y_i)\}_{i \in [n]}\), where \(\boldsymbol{x}_i \in \mathbb{R}^d\) are predictors, assumed drawn i.i.d. from a law \(P_X,\) and \(y_i\) are responses such that \(y_i = f(\boldsymbol{x}_i) + \varepsilon_i,\) for a bounded, centered, independent random noise \(\varepsilon_i\), and bounded \(f\), the goal is to recover an estimate \(\hat{f}\) of \(f\) such that on new data, the squared error \(\mathbb{E}[(y- \hat{f}(\boldsymbol{x}))^2]\) is small.

The statistical challenge of the problem lies in the fact that \(f\) is only weakly constrained - for instance, \(f\) may only be known to be differentiable, or Lipschitz. In addition, the problem is algorithmically challenging in high-dimensions, and many approaches to the univariate problem do not scale well with the dimension \(d\). For instance, piecewise linear regression methods typically involve a prespecified grid, and thus the number of grid points, or knots, grows exponentially with dimension. Similarly, methods like splines typically require both stronger smoothness guarantees and exponentially more parameters to fit with dimension in order to avoid singularities in the estimate.

This paper is concerned with regression over the class of functions that are differences of convex functions, i.e., DC functions. These are functions \(f\) that can be represented as \(f = \phi_1 - \phi_2\) for a choice of two convex functions. DC functions constitute a very rich class - for instance, they are known to contain all \(\mathcal{C}^2\) functions. Such functions have been applied in a variety of contexts including non-convex optimization [1], [2], sparse signal recovery [3] and reinforcement learning [4].

The principal contribution of this paper is a method for piecewise linear regression over the class of DC functions. At the heart of the method is a representation of piecewise linear DC functions via a set of linear constraints, in a manner that generalises the representations used for max-affine regression. The choice of the function is regularised for smoothness by a new seminorm that controls the \(\ell_\infty\)-Lipschitz constant of the resulting function. The resulting estimate is thus a piecewise linear function, represented as the difference of two piecewise linear convex functions, that are smooth in the sense of having bounded gradients.

The method enjoys two main advantages:

1. It is agnostic to any knowledge about the function, and requires minimal parameter tuning.

2. It can be implemented efficiently, via quadratic programming, even in high dimensions. For \(n\) data points in \(\mathbb{R}^d\), the problems has \(2n(2d+1) + 1\) decision variables, and \(n^2\) linear constraints, and can be solved in the worst case in \(O(d^2 n^5)\) time by interior-point methods. Furthermore the algorithm does not need to specify partitions for piece-wise linear parts and avoids ad-hoc generalizations of splines or piece-wise linear methods to multi-dimensions.

In addition, the method is shown to be statistically viable, in that it is shown to attain vanishing risk as the sample size grows at a non-trivial rate, under the condition that the ground truth has bounded DC seminorm. The risk further adapts to structure such as low dimensional supports. Lastly, the solution obtained is a piecewise-linear fit, and thus enjoys interpretability in that features contribution heavily to the value can be readily identified. Further, the fitting procedure naturally imposes \(\ell_1\) regularisation on the weights of the piecewise linear parts, thus enforcing a sparsity of local features, which further improves interpretability.

To establish practical viability, we implement the method on a number of regression and classification tasks. The settings explored are those of moderate data sizes - \(n \le 10^3\), and dimensions \(d \le 10^2\). We note that essesntially all non-parametric models are only viable in these settings - typical computational costs grow with \(n\) and become infeasible for large datasets, while for much higher dimensions, the sample complexity - which grows exponentially with \(d\) - cause small datasets to be non-informative. More pragmatically, all nonparametric methods we compare against have been evaluated on such data. Within these constraints, the method is shown to have better error performance and fluctuation with respect to popular methodologies such as multivariate adaptive regression splines, nearest neighbour methods, and two-layer perceptrons, evaluated on both synthetic and real world data-sets.

1.1 Connections to Existing Methodologies↩︎

Piecewise Linear Regression is popular since such regressors can can model the local features of the data without affecting the global fit. In higher than \(1\) dimensions, piecewise linear functions are usually fit via choosing a partition of the space and fitting linear functions on each part. The principle difficulty thus lies in choosing these partitions. The approach is usually a rectangular grid - for instance, a variable rectangular partition of the space is studied in [5] and solved optimally. However the rectangulization becomes prohibitive in high dimension as the number of parts grow exponentially with the dimension. Other approaches include Bayesian methods such as [6], which rely on computing posterior means for the parameters to be fit via MCMC methods. Max-Affine Regression is a nonparametric approximation to convex regression, originating in [7], [8] that recovers the optimal piece-wise linear approximant to a convex function with the form \(f = \max_{i\in [K]} \langle a_i, \boldsymbol{x}_i \rangle + b_i\) for some \(K\). Smoothness of the estimate can be controlled by constraining the convex function to be Lipschitz. The problem is generic in that it is easily argued that piecewise linear convex functions can uniformly approximate any Lipschitz convex function on a bounded domain. Parametric approaches, i.e., with a fixed \(K\), are popular, but can be computationally intensive due to the induced combinatorics of which of the \(K\) planes is maximised at which data point, and various heuristics and partitioning techniques have to be applied [9]–[11]. The nonparametric case, where \(K\) grows with \(n\), has been extensively analysed in the works [12], [13].

On the other hand, if \(K = n,\) i.e. if the number of affine functions used is the same as the number of datapoints, then the problem becomes amenable to convex programming techniques - when estimating the parameters \(a_i, b_i,\) one can remove the nonlinearity induced by the max, and instead enforce the same via \(n\) linear constraints. This simple fact allows efficient algorithmic approaches to max-affine functions. The heart of our method for DC function regression is an extension of this trick to DC functions. Smoothing splines are an extremely popular regression methodology in low dimensions. The most popular of these are the \(L_2\) smoothing splines, which, in one dimension, involve fixing a ‘knot’ at each data point, and estimating the gradients of the function at each point, with regularisation of the form \(\int |\hat{f}''|^2\). Unfortunately this \(L_2\) regularisation leads to singularities in \(d \ge 3\) dimensions, and methods such as thin plate splines generalising these to higher dimensions resort of regularising up to the \(d\)th derivative of the estimate, leading to an explosion in the number of parameters to be estimated [14].

Our method is closer in relation to \(L_1\) regularised splines, which in the univariate case regularise \(\int |\hat{f}''|\) - it is shown in Appx. 8 that in one dimension our method reduces to these. As a consequence, one may view this method as a new generalisation of the \(L_1\)-spline regressor.

A number of alternative methods for multivariate splines have been proposed, with several, such as general additive models modelling the data via low dimensional projections and assumptions. The most relevant multivariate spline methods are the adaptive regression splines, originating in [15], which is a greedy procedure for recursively refining a partition of the data, and fitting new polynomials over the parts. Previous DC Modeling Finally, let us mention that the final chapter of the doctoral thesis of [13] anticipates our study of DC function regression, but gives neither algorithms nor analyses. Subsequently, [16] introduces the same DC modeling as us in a broader context, where the loss function can also be a DC function. However their problem ends up being non-convex. They focus on developing a majorization-minimization algorithm to find an approximate solution with desirable guarantees such as convergence to a directional stationary solution.

2 A brief introduction to DC functions↩︎

Difference of convex functions are functions which can be represented as difference of two continuous convex function over a domain \(\Omega \subset \mathbb{R}^d\), i.e., the class \[\begin{align} \mathcal{DC}(\Omega) \triangleq \{f: \Omega \to \mathbb{R} \,| \exists &\phi_1,\phi_2 \textrm{ convex, continuous s.t.}\\ &f = \phi_1 - \phi_2 \}. \end{align}\] Throughout the text we will set \(\Omega = \{ \boldsymbol{x} : \| \boldsymbol{x} \|_\infty \le R\}\). We will assume that the noise and the ground truth function are bounded so that \(|\varepsilon|, \sup_{\boldsymbol{x} \in \Omega} |f(\boldsymbol{x})| \le M.\) As a consequence, \(|y| \le 2M\).

One of the first characterizations of DC functions is due to [17]: a univariate function is a DC function if and only if its directional derivatives each has bounded total variation on all compact subsets of \(\Omega\). For higher dimensions it is known that DC functions are a subclass of locally-Lipschitz functions and include \(C^{2}\) functions. Therefore, any continuous function can be uniformly approximated by DC functions. For a recent review see [18].

In the following section we show that D.C functions can fit any sample data. Thus, to allow a bias-variance tradeoff, we regularise the selection of DC functions via the DC seminorm \[\begin{align} \| f\| \triangleq \inf_{\phi_1,\phi_2}&\sup_{\boldsymbol{x}} \sup_{\boldsymbol{v}_i \in \partial_* \phi_i(\boldsymbol{x})} \|\boldsymbol{v}_1\|_1 + \|\boldsymbol{v}_2\|_1 \\ & \textrm{s.t. } \phi_1, \phi_2 \text{ are convex, } \phi_1 - \phi_2 = f, \end{align}\]

where \(\partial_*\phi_i\) denotes the set of subgradients of \(\phi_i\). The above function is not a norm because every constant function satisfies \(\|c\| = 0\). Indeed, if we equate DC functions up to a constant, then the above seminorm turns into a norm. We leave a proof of the fact that the above is a seminorm to Appx. 9.1. Note that the DC seminorm offers strong control on the \(\ell_\infty\)-Lipschitz constant of the convex parts of at least one DC representation of the function (and in turn on the Lipschitz constant of the function).

We will principally be interested in DC functions with bounded DC seminorm, and thus define \[\mathcal{DC}_L \triangleq \{ f \in \mathcal{DC}: \|f\| \le L\}.\]

The bulk of this paper concentrates on piecewise linear DC functions. This is justified because piecewise linear functions are known to uniformly approximate bounded variation functions, and structural results [19], [20] showing that every piecewise linear function may be represented as difference of two convex piecewise linear functions - i.e., max-affine functions. Since the term is used very often, we will abbreviate “piecewise linear DC” as PLDC, and symbolically define \[\begin{align} \textrm{pl-}\mathcal{DC}&\triangleq \{ f \in \mathcal{DC}: f\textrm{ is piecewise linear} \},\\ \textrm{pl-}\mathcal{DC}_{L} &\triangleq \{ f \in \mathcal{DC}_{L}: f\textrm{ is piecewise linear} \}. \end{align}\] The following bound on the seminorm of PLDC functions is useful. The proof is obvious, and thus omitted.

2.1 Expressive Power of Piecewise linear DC functions↩︎

We begin by arguing that PLDC functions can interpolate any finite data. The principle characterisation for DC functions is as follows:

The principal utility of the conditions stated above is that we can utilise these to enforce the shape restriction of getting a DC estimate in an efficient way when performing empirical risk minimisation. Indeed, suppose we wish to fit a DC function to some data \((\boldsymbol{x}_i, y_i)\). We may then introduce decision variables \(\hat{y}_i, z_i, \boldsymbol{a}_i, \boldsymbol{b}_i,\) where the \(\hat{y}_i\) represent the value of our fit at the various \(\boldsymbol{x}_i,\) and then enforce the linear constraints of the above proposition (with \(y_i\) replaced by \(\hat{y}_i\)) while minimising a loss of the form \(\sum (y_i - \hat{y}_i)^2\). Since these constraints are linear, the resulting program is thus convex, and can be efficiently implemented. This observation forms the core of our algorithmic proposal in §3

The above characterisation relies on existence of vectors that may serve as subgradients for the two convex functions. This condition can be removed, as in the following.

Next, in order to contextualise the expressiveness of DC functions, we argue that the popular parametric class of ReLU neural networks can be represented by PLDC functions, and vice versa. This is also argued in [16] for a \(2\) layer network.

where \(|\cdot|\) is the entry-wise absolute value. The above is proved in Appx. 9.2 via an induction over the layers using the relations \[\begin{align} &\max(\max(a,0)\!-\!\max(b,0),0) = \max(a,b,0)\! -\! \max(b,0)\\ &\max(a,b) + \max(c,d) = \max (a+b,a+d, b+c, b+d). \end{align}\]

The proof is constructed in Appx. 9.2 using the relations \[\begin{align} &\max(a,b,c,d) = \max \left(\max (a,b), \max(c,d)\right) \\ &\max(a,b) = \max(a\!-\!b,0) + \max(b,0) -\max(-b,0) . \end{align}\]

3 Algorithms↩︎

We begin by motivating the algorithmic approach we take. This is followed by separate section developing key portions of the algorithm.

Suppose we observe a data-set \(S_n = \{ (\boldsymbol{x}_i, y_i)\}_{i\in[n]}\) generated iid from some distribution \(\mathcal{X}\times\mathcal{Y}\) and a convex loss function \(\ell:\mathbb{R}\times\mathbb{R}\rightarrow \mathbb{R}_+\) bounded by \(c\). The goal is to minimize the expected risk \[\begin{align} \label{erm:dc2} &\min_{f\in\mathcal{DC}} \mathbb{E} [(f(\boldsymbol{x}) - y)^2] \end{align}\tag{3}\] by choosing an appropriate function \(f\) from the class of \(\mathcal{DC}\) functions. Note instead of the squared error, the above could be generalised to any bounded, Lipschitz losses \(\ell(f(\boldsymbol{x}), y)\). Note also that the squared loss is bounded in our setting because of our assumption that both the ground truth and noise are bounded.

There are two basic problems with the above - the distribution is unknown, so the objective above cannot be evaluated, and that the class of DC functions is far too rich, and so the problem is strongly underdetermined. In addition, directly optimising over all DC functions is an intractable problem.

To begin with, we reduce the problem by instead finding the values that a best fit DC function must take at the datapoint \(\boldsymbol{x}_i\), and then fitting a PLDC functions with convex parts that are max-affine over precisely \(n\) linear functionals on this. This has the significant advantage of reducing the optimisation problem to a set of convex constraints. Quantitatively, this step is justified in §4, which argues that the error induced by this approximation via PLDC functions is dominated by the intrinsic risk of the problem.

To handle the lack of knowledge of the distribution, we resort to uniform generalisations bounds in the literature. Our approach to relies on the following result, which mildly generalises the bounds of [21], and uniformly controls the generalisation error of an empirical estimate of the expectation (specialised to our context):

The above statement essentially arises from a doubling trick over a Rademacher complexity based bound for a fixed \(\|f\|\). The broad idea is that since \(\mathcal{DC}_L \subset \mathcal{DC}_{L'}\) for \(L \le L',\) we can separately develop Rademacher complexity based bounds over \(L\) of the form \(2^j\), each having the more stringent high-probability requirement \(\delta_j = \delta 2^{-j}\). A union bound over these then extends these bounds to all of \(\mathcal{DC}= \bigcup_{j \ge 1} \mathcal{DC}_{2^j},\) and for a particular \(f\), the bound for \(j = \lceil \log_2 \|f\| \rceil\) can be used. See §10 for details.

Optimising the empirical upper bound on \(\mathbb{E}[(f(X) - Y)^2]\) implied by the above directly leads to a structural risk minimization over the choice of \(L\). The crucial ingredient in the practicality of this is that for DC functions, \(\hat{D}_n(\mathcal{DC}_L) = L \hat{D}_n(\mathcal{DC}_1)\), and further, \(\hat{D}_n(\mathcal{DC}_1)\) can be computed via linear programming. Thus, the term \(\hat{D}_n\) above serves as a natural, efficiently computable penalty function, and acts exactly as a regularisation on the DC seminorm.

3.1 Computing empirical maximum discrepancy.↩︎

Throughout this and the following sections, we use \(\hat{y}_i\) to denote \(\hat{f}(\boldsymbol{x}_i)\), where \(\hat{f}\) is the estimated function.

The principle construction relies on the characterisation of Proposition 2.

3.2 Structural Risk Minimisation↩︎

To perform the SRM, we again utilize the structural result of Proposition \(1\) to determine the values that the optimal estimate takes at each of the \(\boldsymbol{x}_i\). The choice of the values is penalised by the seminorm as \(\lambda L,\) where \(\lambda = 24M\hat{D}_n(\mathcal{DC}_1),\) which may be computed using the program (?? ). Note that the logarithmic term in the generalisation bound (4 ) is typically small, and is thus omitted in the following. This also has the benefit of rendering the objective function convex, as it avoids the \(\log L\) term that would instead enter. If desired, an convex upper bound may be obtained for the same, for instance by noting that \(\sqrt{ \max(1, \log_2 \|f\|)} \le 1 + \|f\|\). This has the effect of bumping up the value of \(\lambda\) required by \(O(M^2/\sqrt{n})\). However, theoretically this term is strongly dominated by the \(\hat{D}_n\) (see §4), while practically even the value \(\lambda = \hat{D}_n(\mathcal{DC}_1)\) tends to produce overly smoothened solutions (see §5).

With the appropriate choice of \(\lambda,\) this yields the following convex optimisation problem, \[\begin{align} \label{program:SRM} &\qquad \min_{\hat{y}_i, z_i, \boldsymbol{a}_i, \boldsymbol{b}_i, L} \sum_{i = 1}^n ( \hat{y_i} - y_i)^2 + \lambda L\\ &\textrm{s.t.} \begin{cases} \hat{y}_i - \hat{y}_j + z_i - z_j\geq \langle \boldsymbol{a}_j, \boldsymbol{x}_i-\boldsymbol{x}_j\rangle & i,j \in [n] \\ z_i - z_j \geq \langle \boldsymbol{b}_j, \boldsymbol{x}_i-\boldsymbol{x}_j \rangle & i,j \in [n] \\ \|\boldsymbol{a}_i\|_1 + \|\boldsymbol{b}_i\|_1 \leq L & \quad i\in[n] \end{cases}\notag \end{align}\tag{5}\]

Once again, in the above constraints, the first two are necessary and sufficient to ensure that a DC function taking the values \(\hat{y}_i\) at \(\boldsymbol{x}_i\) exists, with the vectors \(\boldsymbol{a}_i, \boldsymbol{b}_i\) serving as witnesses for the subgradients of the convex parts of such a function at the \(\boldsymbol{x}_i\), and the third constraint enforces that the function has seminorm at most \(L\). Notice that the third condition effectively imposes \(\ell_1\)-regularisation on the weight vectors \(\boldsymbol{a}_i, \boldsymbol{b}_i\). This causes these weights to be sparse.

Finally, we may use the witnessing values \(\hat{y}_i, z_i\) and \(\boldsymbol{a}_i, \boldsymbol{b}_i,\) to construct, in the manner of Proposition 2, the following PLDC function, which interpolates the values \(\hat{y}_i\) to the entirety of the domain. Notice that since this function has the same loss as any DC function that satisfies \(f(\boldsymbol{x}_i) = y_i\), this construction enjoys risk bounds constructed above. \[\label{eqn:pldc95estimate} \begin{align} \hat{f}(\boldsymbol{x}) \triangleq & \max_{i\in [n]} \langle \boldsymbol{a}_i, \boldsymbol{x}-\boldsymbol{x}_i \rangle + \hat{y}_i + z_i \\ &\quad -\max_{i\in[n]} \langle \boldsymbol{b}_i, \boldsymbol{x}-\boldsymbol{x}_i \rangle + z_i \end{align}\tag{6}\]

3.3 Computational Complexity↩︎

Training First, we note that we may replace the constraints on the \(1\)-norms of the vectors \(\boldsymbol{a}_i, \boldsymbol{b}_i\) in the above by linear constraints via the standard trick of writing the positive and negative parts of each of their components separately. Overall, this renders the program (?? ) as a linear programming problem over \(2n(2d + 1)\) variables, and with \(n^2\) non-trivial constraints. Note that in our setting, one typically requires that \(n \ge d\) - that one has more samples than dimensions. Via interior point methods, this problem may be solved in \(O(n^5)\) time. For the least squares loss \(\ell(y,\hat{y}) = (y - \hat{y})^2,\) the second program (5 ) is a convex quadratic program when the \(1\)-norm constraints are rewritten as above. The decision variables are the same as the first problem, with the addition of the single \(L\) variable, and the constraints remain identical. Again, via interior point methods, these programs can be solved in \(O(d^2n^5)\) time (see Ch. 11 of [22]). The latter term dominates this complexity analysis. We note that in practice, these problems can be solved significantly faster than the above bounds suggest.

Speeding up training via a GPU implementation an iterative solver for a modified version of the SRM in program (5 ) is given in Algorithm ([program:SRM95ADMM]) in the Appx 12 via the ADMM method [23]. Each iteration of this algorithm can be distributed to \(O(n^2+nd)\) parallel processors. A python implementation is given in our GitHub repository 2 which is compatible with GPU’s. We note that a similar algorithm for that of Lipschitz convex regression is provided in [13], [24] however not all the ADMM blocks are solved in closed form and require additional optimization in each iteration.

Prediction By appending a \(1\) to the input, and the constants \(y_i + z_i - \langle \boldsymbol{a}_i , \boldsymbol{x}_i \rangle\) and \(z_i - \langle \boldsymbol{b}_i, \boldsymbol{x}_i\rangle\) to \(\boldsymbol{a}_i\) and \(\boldsymbol{b}_i\), we can reduce the inference time problem to solving two maximum inner product search problems over \(n\) vectors in \(\mathbb{R}^{d+1}\). This is a well developed and fast primitive, e.g. [25] provide a locality sensitive hashing based scheme that solves the problem in time that is sublinear in \(n\).

4 Analysis↩︎

We note that this analysis makes extensive use of the work of [12], [13] on convex and max-affine regression, with emphasis on the latter thesis, which contains certain refinements of the former paper.

In this section, we assume that the true model \(y = f(\boldsymbol{x}) + \varepsilon_i\) holds for a \(f\) that is a DC function, and that we have explicit knowledge that \(\| f\| \le L.\) Also recall our assumption that the distribution is supported on a set that is contained in the \(\ell_\infty\) ball of radius \(R\). We begin with a few preliminaries

4.0.0.1 A lower bound on risk

The minimax risk of estimation under the squared loss is \(\Omega(n^{-2/d+2}).\) This follows by setting \(p = 1\) (i.e., Lipschitzness) in Theorem 3.2 of [26], which can be done since the standard constructions of obstructions to estimators used in showing such lower bounds all have regularly varying derivatives, and thus finite DC seminorms.

4.0.0.2 PLDC solutions are not lossy

Lemma 5.2 of [13] argues that for every convex \(L\)-Lipschitz functions \(\phi\) with Lipschitz constant¹, \(\sup_x \|\nabla_* \phi\|_1\) bounded by \(L\), there exists a Lipschitz max-affine function \(\phi_{\mathrm{PL}}\) with maximum over \(n\) pieces such that \(\| \phi - \phi_{\mathrm{PL}}\| \le 36LR n^{-2/d}\). Recall that PLDC functions can be represented as differences of max-affine functions, and DC functions as differences of convex functions. Since the DC seminorm controls the \(1\)-norm of the subgradients, which dominates the \(2\)-norm, it follows that for every \(\mathcal{DC}_L\) function \(f\), there exists a \(\textrm{pl-}\mathcal{DC}_{2L}\) function \(f_{\mathrm{PL}}\) with \(\| f - f_{\mathrm{PL}}\|_\infty\! =\! O(n^{-2/d}).\) Note that the resulting excess risk in the squared error due to using PLDC functions can thus be bounded as \(O(n^{-4/d}),\) which is \(o(n^{-2/d+2}),\) i.e., it is dominated by minimax risk.

4.1 Statistical risk↩︎

The bound (4 ) provides a instance specific control on the generalisation error of an estimate via the empirical maximum discrepancy \(\hat{D}_n\). This section is devoted to providing generic bounds on the same under the assumption of i.i.d. sampled data. We adapt the analysis of [13] for convex regression in order to do so. The principal result is as follows.

4.1.0.2 On suboptimality

Comparing to the lower bound, we find that the above rate is close to being minimax, although the (multiplicative) gap is about \(n^{4/d^2}\) and diverges polynomially with \(n\). In analogy with the observations of [13] for max-affine regression, we suspect that this statistical suboptimality can be ameliorated by restricting the PLDC estimate to have some (precisely chosen) \(K_n < n\) hyperplanes instead of the full \(n\). However, as discussed in §1.1, such restrictions lead to increase in the computational costs of training, and thus, we do not pursue them here.

5 Experiments↩︎

In this section we apply our method to both synthetic and real datasets for regression and multi-class classification. The datasets were chosen to fit in the regime of \(n\le 10^3, d \le 10^2\) as described in the introduction. All results are averaged over \(100\) runs and are reported with the \(95\%\) confidence interval.

For the DC function fitting procedure, we note that that the theoretical value for the regularization weight tends to over-smooth the estimators. This behaviour is expected since the bound (4 ) is designed for the worst case. For all the subsequent experiments we make two modifications - since none of the values in the datasets observed are very large, we simply set \(12M = 1,\) and further, we choose the weight, i.e. \(\lambda\) in (5 ), by cross validation over the set \(2^{-j} \hat{D}_n(\mathcal{DC}_1)\) for \(j \in [-8:1]\). Fig. 1 presents both these estimates in a simple setting where one can visually observe the improved fit. Note that this tuning is still minimal - the empirical discrepancy of \(\mathcal{DC}_1\) fixes a rough upper bound on the \(\lambda\) necessary, and we explore only a few different scales.

For the regression task we use the \(L_1\) empirical loss in our algorithm, instead of \(L_2\). That is, the objective in (5 ) is replaced by \(\sum |y_i - \hat{y}_i|\). This change allows us to implement the fitting program as a linear programming problem and significantly speeds up computation. However, in the following we will only report the \(L_2\) error of the solutions obtained thi way. We compare our method to a multi-layer perceptron (neural network) with variable number of hidden layers chosen from \(1:10\) by \(5\)-fold cross validation, Multivariate Adaptive Regression Splines (MARS) and \(K\)-nearest neighbour(\(K\)-NN) where Best value of \(K\) was chosen by \(5\)-fold cross validation from \(1:10\).

For the multi-class classification task we adopt the multi-class hinge loss to train our model, i.e. the loss \[\begin{align} \sum_i \sum_{j\neq y_i}^m\max( f_j(\boldsymbol{x}_i)-f_{y_i}(\boldsymbol{x}_i) +1, 0), \end{align}\] where \(m\) is the number of classes and \(f_j\)’s are DC functions. We compare with the same MLP as above but trained with the cross entropy Loss, KNN and a one versus all SVM.

In both cases, we have used MATLAB Statistics and Machine learning Library for their implementation of MLP, KNN and SVM. For MARS we used the open source implementation in ARESLab toolbox implemented in MATLAB. Our code along with the other algorithms is available in our GitHub repository²
. In each of the following tasks, we observe that our method performs competitively to all considered alternatives in almost every dataset, and often outperforms them, across the variation in dimensionality and dataset sizes.

5.0.0.1 Regression on Synthetic Data

We generated data from the function, \[\begin{align} y &= f(\boldsymbol{x}) + \varepsilon\\ f(\boldsymbol{x}) &= \sin\big(\frac{\pi}{\sqrt{d}}\sum_{j=1}^d x_{j}\big)+\big(\frac{1}{\sqrt{d}}\sum_{j=1}^d x_{j}\big)^2, \end{align}\] where the \(\boldsymbol{x}\) is sampled from a standard Gaussian, while \(\varepsilon\) is a centred Gaussian noise with standard deviation of \(0.25\).

We generate \(50\) points for training. For testing we estimate the Mean Squared Error based on a test set of \(5000\) points without the added noise. We normalize the MSE by variance of the values of test data and multiply by \(100\). Results are presented in Fig. 2. Observe that our algorithm consistently outperforms the competing methods, especially as we increase the dimension of the data. Furthermore our algorithm has lower error variance across the runs.

5.0.0.2 Regression on Real Data

We apply the stated methods to various moderately sized regression datasets that are available in the MATLAB statistical machine learning library. The results are presented in Fig. 3.

In the plot, the datasets are arranged so that the dimension increases from left to right. We observe that we do comparably to the other methods for some datasets and outperform in others. See Appx. D for a description of each of the datasets studied.

5.0.0.3 Multi-class classification

We used popular UCI classification datasets for testing our classification algorithm. We repeated the experiments \(100\) times We present the mean and \(95\%\) C.I.s on a 2-fold random cross validation set, in Fig. 4. We observe to perform comparably to other algorithms on some datasets and outperform in others.

6 Discussion↩︎

The paper proposes an algorithm to learn piecewise linear functions using difference of max-affine functions. Our model results in linear or convex programs which can be solved efficiently even in high dimensions. We have shown several theoretical results on expressivity of our model, as well as its statistical properties, including good risk decay and adaptivity to low dimensional structure, and have demonstrated practicality of our resulting model under regression and classification tasks.

6.0.0.1 Wider context

Non-parametric methods are most often utilised in settings with limited data in moderate dimensions. Within this context, along with strong accuracy, it is often desired that the method is fast, and is interpretable, especially in relatively large dimensions.

In settings with large datasets, accuracy is relatively easy to achieve, and speed at inference is more important. In low dimensional settings, this makes methods such as MARS attractive, due to their fast inference time, while in high-dimensions MLPs are practically preferred. In settings with low dimensionality and small datasets, interpretability and speed take a backseat due to the small number of features, while accuracy become the critical requirement, promoting use of kernelised or nearest neighbour methods.

On the other hand, for small datasets in moderate to high dimensions, interpretability gains an increased emphasis. Within this context, our method results in a piecewise linear fit for which it is easy to characterise the locally important features, and further does so with relatively sparse weights via the \(\ell_1\) regularisation on \(\boldsymbol{a}_i, \boldsymbol{b}_i\). Further, since the suboptimality of our statistical risk is controlled as \(n^{O(1/d^2)},\) as the dimension increases, our method gets closer to the optimal accuracy. This thus represents the natural niche where application of DC function regression is appropriate.

7 Acknowledgment↩︎

This research was supported by NSF CAREER Award 1559558, CCF-2007350 (VS), CCF-2022446 (VS), CCF-1955981 (VS) and the Data Science Faculty Fellowship from the Rafik B. Hariri Institute.

References↩︎