We study a natural lasso variant of the multivariate adaptive regression splines (MARS) method (see Friedman (1991)o rHastie, Tibshirani and Friedman (2009), Section 9.4) for nonparametric regression. To understand the relationship between a response variable y and d explanatory variables x 1 ,...,x d based on observed data (x (1) ,y 1 ),...,(x (n) ,y n ) with x (i) ∈ R d and y i ∈ R, MARS fits a function y = fmars (x 1 ,...,x d ) where fmars is a sparse linear combination of functions of the form (1)
b j (x j ) α j = j :α j =1 b j (x j ) with α = (α 1 ,...,α d ) ∈{0, 1} d and b j (x j ) = (x jt j ) + or (t jx j ) + for some real number t j .
Here • + := max{•, 0} indicates the ReLU function. As a concrete example, to understand the relationship between the logarithm of Weekly Earnings (y) and the two variables, Years of Education (x 1 ) and Years of Experience (x 2 ), from a standard dataset (ex1029 in the R library Sleuth3) collected from 25,437 full-time male workers in 1987, the default implementation of MARS from the R package earth (with the maximum degree of interaction set to two) fits the function 5.83 + 0.0695(x 1 -7) + -0.0370(11x 1 ) + + 0.0155(x 2 -13) + -0.0600(13x 2 ) + -0.0164(x 2 -30) + -0.0114(x 1 -11) + (x 2 -40) + + 0.00148(x 1 -11) + (40
which is clearly a linear combination of the eight functions each of the form (1). MARS fits a nonlinear function to the observed data that is simple enough to be interpretable because it is built from the basic ReLU functions (x jt j ) + and (t jx j ) + .Furthermore, MARS fits non-additive functions because of the presence of products in (1), which enables interactions between the explanatory variables x 1 ,...,x d . Indeed, the term (1) can be interpreted as an interaction term of order |α| between the variables in the set S(α).Here and in the rest of the paper, we use the notation The exact methodology that MARS uses involves a greedy algorithm similar to stepwise regression methods. Specifically, one adds in basis functions of the form (1) starting with a constant function using a goodness of fit criterion. Typically, one only considers terms (1)for which the interaction order |α| is smaller than a pre-chosen integer s ≤ d (most commonly s = 1ors = 2). Once a reasonably large number of basis functions with |α|≤s are added, a backward deletion procedure is applied to remove superfluous basis functions. We refer the reader to Hastie, Tibshirani and Friedman (2009), Section 9.4 for more details on MARS.
Our goal in this paper is to propose and study a lasso variant of the MARS method where we consider all the basis functions of the form (1) with |α|≤s and apply the lasso method of Tibshirani (1996). As is well known, lasso is an attractive alternative to stepwise regression methods in usual linear models. In order to apply lasso in the MARS setting, we first assume that the explanatory variables x 1 ,...,x d all take values in the interval [0, 1].Inotherwords, the domain of the regression function is assumed to be [0, 1] d . In practical settings, this can be achieved by subtracting the minimum possible value and dividing by the range for each explanatory variable. This scaling puts all the variables on a comparable footing enabling the application of lasso. Without such a scaling, x 1 ,...,x d might be on very different scales in which case penalizing or constraining the sum of the absolute values of the coefficients corresponding to the terms (x jt j ) + would be unnatural (e.g., think of the setting where x 1 is years of education and x 2 is days of experience). After fitting a function in the transformed domain [0, 1] d , we can simply invert the transformation to find the equation of the fitted function in the original domain (see Section 6 for some real data applications).
In the rest of the paper, we assume that the explanatory variables x 1 ,...,x d all belong to [0, 1]. The observed data is (x (1) ,y 1 ),...,(x (n) ,y n ) where x (i) ∈[0, 1] d and y i ∈ R.T othis data, we fit functions of the form y = f(x 1 ,...,x d ) where f :[0, 1] d → R via the application of lasso with the MARS basis functions. The restriction x j ∈[0, 1] allows us to make two simplifications to the usual MARS setup:
(i) Instead of considering both kinds of functions (x jt j ) + and (t jx j ) + , we only take into account (x jt j ) + , because as each x j isassumedtobein[0, 1], we can write (t jx j ) + = (x jt j ) +x j + t j = (x jt j ) + -(x j -0) + + t j , which implies that every linear combination of functions of the form (1) is also a linear combination of functions of the same form (1)whereb j (x j ) = (x jt j ) + for some t j .
(ii) We assume that t j ∈[0, 1) for each j . This is because when t j ≥ 1, the function (x jt j ) + becomes 0 as x j ∈[0, 1],a n df o rt j < 0, the function (x jt j ) + = x jt j = (x j -0) +t j is a linear combination of (x j -0) + and the constant function 1.
Because there are an uncountable number of functions of the form (1)( a st j can be any number in [0, 1)), we work with an infinite-dimensional version of lasso. Infinite-dimensional lasso formulations have been used in many papers including Rosset et al. (2007), de Castro and Gamboa (2012), Bredies and Pikkarainen (2013), Candès and Fernandez-Granda (2014), Duval and Peyré (2015), De Castro et al. (2017), Denoyelle et al. (2020), and Condat (2020), which studied various inverse problems in spaces of measures. The main idea is to consider infinite linear combinations of basis functions that are parametrized by signed measures and to measure complexity in terms of the variations of the involved signed measures. In the MARS context, infinite linear combinations of the basis functions (1) with |α|≤s are
where 0 := (0,...,0), a 0 ∈ R, ν α is a finite (Borel) signed measure on [0, 1) |α| ,a n dt (α) indicates the vector (t j ,j ∈ S(α)) for each binary vector α ∈{0, 1} d \{0} with |α|≤s.W e will denote the collection of all such functions f a 0 ,{ν α } by F d,s ∞-mars (the subscript ∞ highlights the fact that F d,s ∞-mars contains infinite linear combinations of the functions (1)). The usual MARS functions are special cases of (3) corresponding to discrete signed measures ν α . Indeed, when each ν α is supported on a finite set {(t
x jt (α) lj + .
Our infinite-dimensional lasso estimator minimizes the least squares criterion over
∞-mars with a constraint on the complexity of f a 0 ,{ν α } . The complexity measure involves the sum of the variations of the underlying signed measures ν α a n di sa ni n fi n i t edimensional analogue of the usual L 1 norm of the coefficients used in finite-dimensional lasso. Recall that, for a signed measure ν on and a measurable subset E ⊆ ,t h ev a r i ation of ν on E is denoted by |ν|(E) and is defined as the supremum of A∈π |ν(A)| over all partitions π of E into a countable number of disjoint measurable subsets. Using the variation of the involved signed measures, we define our complexity measure for functions
We are excluding 0 = (0,...,0) in the variation of ν α because we want to only penalize those basis functions that include at least one nonlinear term and leave unpenalized basis functions that are products of linear functions (note that (x jt j ) + = x j is linear when t j = 0b ecause x j ∈[0, 1]). In Section 11.1 of the Supplementary Material (Ki, Fang and Guntuboyina (2024)), we show that this complexity measure is well defined by proving the uniqueness of the representation f = f a 0 ,{ν α } .
To see why ( 6) is a generalization (to infinite linear combinations) of the notion of L 1 norm over the coefficients of the finite linear combination (5), just note that when ν α is the discrete signed measure given by (4), we have
which is simply the sum of the absolute values of the coefficients in (5) corresponding to the basis functions that have at least one nonlinear term in their product.
Our estimator is thus given by ( 7)
for a tuning parameter V>0. We prove that f d,s n,V exists and can be computed by applying finite-dimensional lasso algorithms to the finite basis of functions obtained by placing the following restrictions on the knots t j in (1):
Here we use the notation
d ) for the ith design point x (i) . As the finitedimensional lasso estimation procedure usually zeros out many regression coefficients, it enables us to obtain f d,s n,V that is a sparse linear combination of (1). Therefore, our estimation procedure can be seen as an alternative to the usual MARS procedure. It is interesting to note that the usual MARS algorithm also works with the restriction (8) on the knots, although typically no theoretical justification is provided for this reduction. We also introduce a computationally more efficient approximate version f d,s n,V of f d,s n,V which seems to work nearly as well in practice. The approximate version f d,s n,V is obtained by restricting the knots t j as
for some pre-selected positive integers N 1 ,...,N d .Forlargen, f d,s n,V can be computed much more efficiently than f d,s n,V . We study the theoretical accuracy of these estimators for an unknown regression function f * under the standard regression model:
where ξ i are mean zero errors whose distributions satisfy certain restrictions. We work with both the fixed design setting where x (1) ,...,x (n) form a lattice in [0, 1] d , as well as the random design setting where x (1) ,...,x (n) are assumed to be realizations of i.i.d. random variables. In the former lattice design setting, which is restrictive but standard in nonparametric function estimation (see, e.g., Nemirovski (2000)), we analyze the non-asymptotic accuracy of f d,s n,V and f d,s n,V . In the latter random design setting, we study their accuracy asymptotically. Our theoretical results show that these estimators achieve rates of convergence of the form n -4/5 (log n) as+b for some constants a and b. It is already known that in the univariate case (d = s = 1), the estimator f 1,1 n,V achieves the rate n -4/5 (see, e.g., Mammen and van de Geer (1997), Theorem 10 or Guntuboyina et al. (2020), Theorem 2.1). Thus, our results imply that in going from the univariate to the multivariate setting, the rate of convergence only deteriorates by a logarithmic multiplicative factor. This suggests that our lasso method for MARS fitting avoids the usual curse of dimensionality to some extent and can thus be an effective function estimation technique in higher dimensions.
We can see why our estimators achieve the dimension-free rates (up to the logarithmic multiplicative factors) in part from an alternative characterization of f d,s n,V . We can characterize f d,s n,V alternatively as a least squares estimator over a class of functions whose smoothness order, in a certain sense, grows with the dimension d. A key role in this characterization is played by mixed partial derivatives of order 2. For an integer k ≥ 1 and a real-valued function f defined on [0, 1] m , by the mixed partial derivatives of f of order k, we mean ( 9)
, where β is an m-dimensional nonnegative integer vector with max j β j = k. Whenever we use the notation f (β) , we inherently assume that f is sufficiently smooth, so that the right-hand side of ( 9) is irrespective of the order of differentiation and f (β) is well defined. Using mixed partial derivatives, we prove the following alternative characterization of f d,s n,V :
where
The main condition here is that the sum of the L 1 norms of mixed partial derivatives of order 2 is at most V .T h es e t T (β) appearing in the integral signifies that the integral of the mixed partial derivative f (β) is only over those coordinates x l for which β l = max j β j (the remaining coordinates are set to zero). Also, the condition f (α) = 0for|α| >s rules out interactions of order greater than s. This characterization shows that the maximum total order β 1 +•••+β d of the mixed partial derivatives appearing in the constraint equals 2d.Int h i s sense, the smoothness order of the constraint can be taken to be 2d, which increases with the dimension d and explains the dimension-free (up to the logarithmic multiplicative factors) rates of convergence. It should be noted however that not all (in fact, only one) mixed partial derivatives of total order 2d are considered in the constraint, and this keeps the function class being too small or restrictive. Also, it should be mentioned that it is well known from approximation theory that L p norm constraints on mixed partial derivatives are advantageous and allow one to overcome the curse of dimensionality to some extent from the perspective of metric entropy, approximation, and interpolation (see, e.g., Bungartz and Griebel (2004), Dũng, Temlyakov andUllrich (2018), Temlyakov (2018)).
In fact, the smoothness characterization (10) is not fully rigorous. Functions of the form (2) clearly belong to the constraint set in (7), but they do not belong to the constraint set in (10) because mixed partial derivatives of order 2 do not exist for these functions. We fix this problem by interpreting the L 1 norms of mixed partial derivatives of order 2 in terms of the Hardy-Krause variations of particular derivatives that we will define in Section 4. Hardy-Krause variation (see, e.g., Aistleitner and Dick (2015), Owen (2005)) is a multivariate generalization of total variation of univariate functions (we review the definition of Hardy-Krause variation and its properties in Section 7 of the Supplementary Material (Ki, Fang and Guntuboyina (2024))). Thus, even though (10) is not fully rigorous because mixed partial derivatives of order 2 do not exist for many important MARS functions, it is still helpful for understanding how the curse of dimensionality can be avoided by our estimators.
The characterization (10) also connects our estimators to other related methods from the literature. In the univariate case (d = s = 1), we have
which is a constrained analogue of the locally adaptive regression splines estimator of M a m m e na n dv a nd eG e e r(1997) when the order k (in their notation) equals 2. Hence, our estimator f d,s n,V can be seen as a multivariate generalization of this univariate estimator of Mammen and van de Geer (1997). Furthermore, if s = d and the condition max j β j = 2 is replaced with max j β j = 1i n( 10), then one obtains the Hardy-Krause variation denoising estimator of Fang, Guntuboyina and Sen (2021). Therefore, we can also view f d,d n,V as a second-order Hardy-Krause variation denoising estimator. Further connections to related work are detailed in Section 5.
We would like to point out here that the theoretical rates of convergence as well as the smoothness characterization have been made possible due to our infinite-dimensional lasso formulation of MARS. In contrast, to the best of our knowledge, no rates of convergence are known for the usual MARS method. Also, there exist no prior connections between the usual MARS method and nonparametric regression methods based on smoothness assumptions.
In addition to the above theoretical contributions, we also implement our method with a cross-validation scheme for the selection of the tuning parameter V and compare our estimators to the usual MARS estimator using simulated and real data.
The rest of the paper is organized as follows. In Section 2, we present results on the existence and computation of f d,s n,V and also introduce the approximate version f d,s n,V . Theoretical accuracy results for f d,s n,V and f d,s n,V are in Section 3. Section 4 is devoted to the alternative characterization (10) based on smoothness. In Section 5, we discuss connections between our method and other related methods. In Section 6, we illustrate the performance of our method in simulated and real data settings and compare its performance to that of the usual MARS algorithm.
In this section, we prove the existence of our infinite-dimensional lasso estimator f d,s n,V (defined in ( 7)) and show that it can be computed via finite-dimensional lasso algorithms. We also introduce a computationally more efficient approximate version of our estimator.
We start with the observation that the objective function of the optimization problem defined in (7) only depends on the function f through its values at the design points
As proved in the next lemma, this observation allows us to restrict our attention to the finitedimensional subclass of F d,s ∞-mars consisting of the functions (3) where each ν α is a discrete signed measure supported on the lattice generated by the design points. For each k ∈[d],let U k denote the finite subset of [0, 1] consisting of the points 0,x (1) k ,...,x (n) k , 1 (recall here that x (i) k denotes the kth coordinate of the ith design point
. As there could be ties among 0,x (1) k ,...,x (n) k , 1, we will write, for some
0 is always 0, and that u (k) n k is always 1. The next lemma (proved in Section 11.2.1 of the Supplementary Material (Ki, Fang and Guntuboyina (2024))) implies that, for the optimization problem (7), we can restrict to the functions of the form (3) where each ν α is a discrete signed measure supported on the finite set ( k∈S(α) U k ) ∩[0, 1) |α| . LEMMA 2.1. Suppose we are given a real number a 0 and a collection of finite signed measures {ν α } where ν α is defined on [0, 1) |α| for each α ∈{0, 1} d \{0} with |α|≤s. Then, there exists a collection of discrete signed measures {μ α } where μ α is concentrated on
When ν α is concentrated on ( k∈S(α) U k ) ∩[0, 1) |α| for each α, the function f a 0 ,{ν α } can be written as ( 12)
where we use the notation [p : q]:={p, p + 1,...,q} for two integers p ≤ q. Also, its complexity measure becomes
The above function ( 12) is a linear combination of the basis functions (1) whose knots t k are chosen from U k \{1}={0,x (1) k ,...,x (n) k }, and its complexity measure equals the absolute sum of the coefficients of the involved basis functions with at least one nonlinear term. Thus, if we additionally assume that f in the problem ( 7) is constructed from discrete signed measures as above, then (7) reduces to a finite-dimensional lasso problem. Lemma 2.1 then implies that every solution to this finite-dimensional lasso problem is also a solution to (7). A precise statement is given in the following result, which we prove in Section 11.2.2 of the Supplementary Material (Ki, Fang and Guntuboyina (2024)
and let M be the n ×|J | matrix with columns indexed by (α, l) ∈ J such that
where 1 := (1,...,1) and y = (y i ,i ∈[n]) is the vector of observations. Then, the function
is a solution to the problem (7). The problem (7) can have multiple solutions, but every solution
Because the set
is closed and convex, there exists a solution to the finite-dimensional lasso problem (13).
Hence, the existence of solutions to our estimation problem ( 7) is guaranteed by Proposition 2.2. Also, once we find a solution to the problem (13) via any optimization algorithms, we can construct a solution f d,s n,V to the problem (7) through the equation ( 14). However, solving the finite-dimensional lasso problem ( 13) can be computationally intensive if n is large because the number of columns of M equals
which is of order O(n s ) (ignoring a multiplicative factor in d) in the worst case when each
In the current implementation of our method, we utilize the optimization software MOSEK as a black-box tool for solving the problem (13) (see Section 6 for more details). Using this black-box tool involves creating the whole matrix M, and thus, when n is large, our current implementation not only requires a large amount of space for this matrix but also often consumes most of running time constructing it. This limitation motivates us to come up with the following approximate method. As we have seen above, Lemma 2.1 ensures that we only need to consider discrete signed measures ν α supported on the lattices ( k∈S(α) U k ) ∩[0, 1) |α| for our estimation problem (7). In the approximate method, we instead restrict our attention to discrete signed measures ν α supported on the lattices generated by
for some pre-selected positive integers N 1 ,...,N d , and we only take into consideration the basis functions corresponding to those signed measures. Note that in contrast to U k whose cardinality are of order O(n) in the worst case, the cardinality of each set Ũk is always N k + 1 regardless of the design points x (1) ,...,x (n) .
We then consider the finite-dimensional optimization problem to which the problem (7) reduces when we additionally impose such restrictions on signed measures ν α . We call this problem the approximate (finite-dimensional optimization) problem. The approximate problem has the same form as (13) but with different M and J .Here
and M is the n ×|J | matrix with columns indexed by (α, l) ∈ J such that
As opposed to the original finite-dimensional problem ( 13), the number of columns of M in this problem is always fixed and not affected by the design points x (1) ,...,x (n) . Hence, this approximate problem can be solved much efficiently than (13), especially when n is large.
Once we find a solution to the approximate problem, we can construct an estimator of the true underlying function f * through the equation ( 14) as before. We denote this estimator by f d,s n,V and call it an approximate version of f d,s n,V . In the next section, we study the theoretical accuracy of f d,s n,V along with f d,s n,V . We will see that if we choose N k appropriately, the approximate method is as accurate as the original method, while it significantly improves computational efficiency.
3. Risk analysis. This section is dedicated to the study of the theoretical accuracy of f d,s n,V and f d,s n,V as an estimator for unknown regression functions. We first consider the nonasymptotic accuracy of f d,s n,V and f d,s n,V in the fixed design setting and then study their asymptotic accuracy in the random design setting. The proofs of all the results in this section are provided in Section 11.3 of the Supplementary Material (Ki, Fang and Guntuboyina (2024)).
3.1. Fixed design. Hereweassumethatx (1) ,...,x (n) form a lattice (15)
x (1) ,...,
where for every k ∈
for some constant ρ>0. We also assume that y 1 ,...,y n are generated according to the regression model ( 16)
where f * :[ 0, 1] d → R is an unknown regression function and ξ i are independent sub-Gaussian errors with mean zero and with a sub-Gaussian parameter σ ,thatis,
for all λ ∈ R. We measure the accuracy of an estimator fn of f * via the squared empirical
and define its risk as
where the expectation is taken over y 1 ,...,y n .
Our first result states an upper bound of the risk of f d,s n,V under the assumption
∞-mars and V mars (f * ) ≤ V and assume the lattice design (15). The estimator f d,s n,V then satisfies that
for some positive constant C ρ,d depending on ρ and d.
Note that for fixed ρ, d, σ ,andV and sufficiently large n, the first term is the dominant term on the right-hand side of (18), so that
, where the multiplicative constant underlying O(•) depends on ρ, d, σ ,andV . In the univariate case, we can deduce from Guntuboyina et al. (2020), Theorem 2.1 that
where C ρ is a positive constant depending on ρ.Inotherwords,
where the multiplicative constant underlying O(•) depends on ρ, σ ,andV . Thus, what Theorem 3.1 tells us is that for general d and s, f d,s n,V can achieve the same rate n -4/5 , although it slightly deteriorates by a logarithmic multiplicative factor depending on s. This suggests that our lasso method for MARS fitting can avoid the curse of dimensionality to some extent and be a useful estimation technique in higher dimensions.
The key step of our proof of Theorem 3.1 is to find an upper bound of the metric entropy of D m (under the L 2 norm), which is defined as the collection of all the functions of the form for every 0 < < m . The logarithmic multiplicative factor can be omitted when m = 1.
REMARK 3.3. If the class D m is altered by replacing (xt) + with 1{x ≥ t} and restricting ν to probability measures, one obtains the collection of all the functions of the form
This class of functions is indeed the collection of all probability distributions on [0, 1] m , whose upper bounds on the metric entropy were derived in Blei, Gao and Li (2007). Thus, we are basically extending the argument in Blei, Gao and Li (2007)
Theorem 3.2 is novel to the best of our knowledge even though we use standard tools and techniques for proving it. We first connect upper bounds of the metric entropy of D m to lower bounds of the small ball probability of integrated Brownian sheet based on ideas from Blei, Gao and Li (2007), Section 3 and Gao (2008), Section 3 and results from Li and Linde (1999), Theorem 1.2 and Artstein et al. (2004), Theorem 5. The small ball probability of integrated Brownian sheet here refers to the quantity
where >0a n dX m is an m-dimensional integrated Brownian sheet (a description of integrated Brownian sheet is given in Section 11.3.2 of the Supplementary Material (Ki, Fang and Guntuboyina (2024))). Required bounds on this small ball probability are then obtained using results from Dunker et al. (1999), Theorem 6 and Chen and Li (2003), Theorem 1.2. Specifically, we show that there exist positive constants c m and m depending on m such that log P sup
for every 0 < < m . This result, along with the connection between the metric entropy and the small ball probability, leads to Theorem 3.2, which is the main ingredient in our proof of Theorem 3.1. Now, we turn to the result for the approximate version f d,s n,V . The next theorem presents an upper bound of the risk of f d,s n,V under the same assumption as in Theorem 3.1. Recall that N k are the pre-selected integers for the approximate method.
∞-mars and V mars (f * ) ≤ V and assume the lattice design (15). The estimator f d,s n,V then satisfies that
for some positive constant C ρ,d depending on ρ and d, where N = min k N k .
Theorem 3.4 shows that f d,s n,V has almost the same risk upper bound as f d,s n,V . The only difference is the existence of the approximation error term 8V 2 /N 2 , which converges to 0 as N goes to infinity. Hence, for sufficiently large N , f d,s n,V achieves the same rate as f d,s n,V . Indeed, if we set each N k to be of order at least n 2/5 ,then
, where the multiplicative constant underlying O(•) depends on ρ, d, σ ,andV .
3.2. Random design. Hereweassumethatx (1) ,...,x (n) are realizations of i.i.d. random variables X (1) ,...,X (n) with a probability density function p 0 on [0, 1] d that is bounded by some constant B ≥ 1, that is, p 0 ∞ ≤ B. Also, we assume that (X (1) ,y 1 ),...,(X (n) ,y n ) are generated according to the regression model ( 20)
where ξ i are i.i.d. errors independent of X (1) ,...,X (n) with mean zero and with finite L 5,1 norm; that is, (21)
Note that the condition ( 21) is stronger than the finite fifth-moment condition ξ i 5 < ∞,but weaker than the finite (5 + ) th -moment condition ξ i 5+ < ∞ for every >0 (see, e.g., Ledoux and Talagrand (1991), Chapter 10). In this setting, we measure the accuracy of an estimator fn of f * by
The next theorem presents the rate of convergence of f d,s n,V under the assumption f * ∈ F d,s
∞-mars and V mars (f * ) ≤ V . Note that f d,s n,V still achieves the rate n -4/5 as in the fixed lattice design setting, although the exponent of the logarithmic multiplicative factor is slightly bigger when s>2.
∞-mars and V mars (f * ) ≤ V , then we have
.
AsthemetricentropyofD m played a central role in our proof of Theorem 3.1, the bracketing entropy of D m is the key ingredient of our proof of Theorem 3.5. The following theorem states an upper bound of the bracketing entropy of D m . THEOREM 3.6. There exists a positive constant C m depending on m such that
for every >0, where
REMARK 3.7. Theorem 3.6 also provides an upper bound of the metric entropy of D m (under the L 2 norm). Since
for every >0. However, this upper bound is weaker than the one we achieved in Theorem 3.2. Although it has the same order for , the exponent of the logarithmic multiplicative factor is bigger. We can obtain from this result an upper bound of the risk of f d,s n,V under the fixed lattice design setting, but it will lead to a bound looser than the one in Theorem 3.1.
We can prove a similar result as in Theorem 3.5 for the approximate version f d,s n,V .A s we state in the following theorem, f d,s n,V achieves the same rate of convergence as f d,s n,V if N 1 ,...,N d are sufficiently large. Together with (19), this result suggests that the approximate method with appropriately chosen N 1 ,...,N d can be as accurate as the original method. THEOREM 3.8. Suppose f * ∈ F d,s ∞-mars and V mars (f * ) ≤ V . Also, assume that N = min k N k = (n 4/15 ), that is, there exists a positive constant c B,d,V possibly depending on B, d, and V such that
Then, the estimator f d,s n,V satisfies that
8(s-1) 5
.
Our next result shows that the logarithmic multiplicative factor in (23) can not be completely removed in the minimax sense. Specifically, we bound the minimax risk defined as
where the expectation is taken over (X (1) ,y 1 ),...,(X (n) ,y n ) of ( 20)a n di n f fn denotes the infimum over all estimators fn of f * based on (X (1) ,y 1 ),...,(X (n) ,y n ).H e r ew ef u r t h e r restrict that ξ i in the model ( 20) are independent Gaussian errors with mean zero and variance σ 2 and that the probability density function p 0 of X (i) is bounded below by some positive constant b,thatis, p 0 ∞ ≥ b. Our result shows that the supremum risk of every estimator indeed requires a logarithmic multiplicative factor depending on s in addition to the n -4/5 term. Note though that there is still a gap between the exponent 8(s -1)/5oflogn in the rate of convergence of f d,s n,V and the exponent 4(s -1)/5oflogn in the minimax lower bound.
THEOREM 3.9. There exist positive constants C b,B,s depending on b, B, and s and c B,s depending on B and s such that
Our proof of Theorem 3.9 is based on Assouad's lemma with a finite set of functions in {f * ∈ F d,s ∞-mars : V mars (f * ) ≤ V } that is constructed by an extension of the ideas in Blei, Gao and Li (2007), Section 4. Results similar to Theorem 3.9 can be proved under the fixed design setting, but we do not go into detail in this paper.
In this section, we provide alternative characterizations of F d,s ∞-mars , V mars (•),a n d f d,s n,V in terms of smoothness. To motivate the results for general d and s, let us first consider the univariate case d = s = 1. We include the proofs of all the results in this section in Section 11.4 of the Supplementary Material (Ki, Fang and Guntuboyina (2024)).
∞-mars consists of all the functions f :[0, 1]→R of the form
where a 0 is a real number and ν is a finite signed measure on [0, 1), and the complexity measure of f is given by the variation of ν on (0, 1);thatis,V mars (f ) =|ν|((0, 1)).
The following simple arguments show that F 1,1 ∞-mars can be characterized in terms of smoothness. First, by replacing (xt) + with 1 0 1{t ≤ s ≤ x} ds in the integral in (24)an d changing the order of integration, we obtain
where the function g :[0, 1]→R is given by ( 26)
It can be readily verified that the function g in ( 26) is right-continuous on [0, 1] and leftcontinuous at 1, and the total variation V(g)of g is finite and can be represented as ( 27) V(g)=|ν| (0, 1) .
Here the total variation of a function h :[0, 1]→R is defined by
where the supremum is over all integers k ≥ 1 and partitions 0
Conversely, every function g :[0, 1]→R that is right-continuous on [0, 1],l e f tcontinuous at 1, and has finite total variation can be written as (26) for a unique signed measure ν on [0, 1) (see, e.g., Aistleitner and Dick (2015), Theorem 3). Putting these observations together, we can argue that F 1,1 ∞-mars has the following alternative characterization:
g is right-continuous on [0, 1], left-continuous at 1, V(g)<∞, and f(x)= a 0 + x 0 g(t)dt for all x ∈[0, 1] .
(28)
Moreover, we can see that the complexity measure V mars (f ) for f ∈ F 1,1 ∞-mars is equal to the total variation V(g)of the function g appearing in (28).
For every function f ∈ F 1,1 ∞-mars , we can show that g satisfying the conditions in (28)i s unique, and thus, we can consider such g as a particular derivative of f satisfying (25). If we denote it by D (1) f , the estimator f 1,1 n,V then can be alternatively written as
The representation (25) implies, by the Lebesgue differentiation theorem (see, e.g., Rudin (1987), Theorem 7.10), that f exists and is equal to D (1) f almost everywhere (with respect to the Lebesgue measure) on [0, 1]. Hence, we can also describe f 1,1 n,V somewhat loosely as
The corresponding penalized version (29) argmin
was proposed by M a m m e na n dv a nd eG e e r(1997) as part of the class of estimators collectively called locally adaptive regression splines. In the univariate case, f 1,1 n,V can thus be seen as a constrained analogue of the locally adaptive regression spline estimator of Mammen and van de Geer (1997) when the order k (in their notation) equals 2. Steidl, Didas and Neumann (2006) used the terminology second-order total variation regularization,andKim et al. ( 2009)andTibshirani (2014) used the terminology first-order trend filtering for (29). Therefore, our estimator f d,s n,V can be considered as a multivariate generalization of piecewise linear (second-order) locally adaptive regression splines, second-order total variation regularization, or first-order trend filtering.
From the alternative characterization of F 1,1 ∞-mars given above, it follows that every sufficiently smooth function f :[0, 1]→R belongs to F 1,1 ∞-mars . Indeed, if f and f exist everywhere and are continuous on [0, 1],thenwehave
for all x ∈[0, 1],and
Thus, in this case, f belongs to F 1,1 ∞-mars and the complexity measure of f can be written as
This formula highlights the role of the second derivative f in the determination of V mars (f ) for each sufficiently smooth function f . 4.2. Smoothness characterization for general d and s. As we have seen in the previous subsection, in the univariate case, F 1,1 ∞-mars consists of all the functions f satisfying (25) with some function g having finite total variation and some one-sided continuity. An analogous characterization holds for general d and s. For general d and s, the role of total variation in the univariate case is played by Hardy-Krause variation, which is an extension of total variation of univariate functions to higher dimensions. In Section 7 of the Supplementary Material (Ki, Fang and Guntuboyina (2024)), we review the definition of Hardy-Krause variation and its properties that we will use for proving the results in this subsection. Standard references for Hardy-Krause variation are Aistleitner and Dick (2015)andOwen ( 2005). Here we use Hardy-Krause variation anchored at 0, which we denote by V HK0 (•).
The following result provides an alternative characterization of F d,s ∞-mars and V mars (•) in terms of smoothness. Recall that we use the notation t (α) to indicate the vector (t j ,j ∈ S(α)) for each α ∈{0, 1} d \{0} with |α|≤s. g α t (α) dt (α) for some a 0 ∈ R and some collection of functions {g α : α ∈{0, 1} d \{0} and |α|≤s}, where for each α ∈{0, 1} d \{0} with |α|≤s, ) |α| with respect to all the j th coordinates where x j = 1.
Furthermore, the complexity of f in (31) can be written in terms of the Hardy-Krause variations of g α as
Proposition 4.1 is completely analogous to (28) for the case d = s = 1. Specifically, the condition ( 31) is analogous to the univariate condition (25). The condition V HK0 (g α )<∞ for each α ∈{0, 1} d \{0} with |α|≤s corresponds to the univariate condition V(g)<∞.
The coordinatewise right-continuity of each g α on [0, 1] |α| is matched with the univariate right-continuity on [0, 1]. Lastly, the coordinatewise left-continuity of each g α at each x (α) ∈ [0, 1] |α| \[0, 1) |α| (with respect to all the j th coordinates where x j = 1) is a counterpart of the univariate left-continuity at 1. It is also interesting to note that V mars (f ) equals the sum of the Hardy-Krause variations of g α over α ∈{0, 1} d \{0} with |α|≤s.
For every function f ∈ F d,s
∞-mars , it can be easily checked that g α appearing in Proposition 4.1 are uniquely determined by f . As in the case d = s = 1, we can thus consider such g α as particular derivatives of f satisfying (31). Let us denote them by D (α) f for α ∈{0, 1} d \{0}. We can then write our estimator f d,s n,V alternatively as
Hence, our estimator can be viewed as a least squares estimator under a specific smoothness constraint involving the sum of the Hardy-Krause variations of the particular derivatives defined via D (α) .
Recall that the univariate condition (25) implies that f exists and equals to D (1) f almost everywhere. Similarly, the condition (31) imposes a certain kind of smoothness on f and characterizes the corresponding derivatives in terms of D (α) f . For each α ∈{0, 1} d \{0} and for
if the limit exists, where
if the limit exists. Note that in contrast to mixed partial derivatives f (α) ( d e fi n e di n( 9)), in which partial derivatives ∂/∂x j are taken sequentially, here all the j th coordinates where α j = 1 are considered simultaneously. Also, note that the remaining coordinates where α j = 0are set to zero for D α f . As in the case d = s = 1, we can show that D α f exist and equal to D (α) f almost everywhere (with respect to the Lebesgue measure) on [0, 1] |α| . The precise statement is given in the following result. Proposition 4.1 also implies that every sufficiently smooth function belongs to F d,d ∞-mars . This is proved in the next result, which also gives an expression for V mars (f ) in terms of the L 1 norms of the mixed partial derivatives of f , for sufficiently smooth functions f .
The following notation will be used below. For each α ∈{0, 1} d \{0},weletJ α be the set of all β ∈{0, 1, 2} d such that (34) max
Also, recall the notation T (β) from ( 11).
LEMMA 4.3. Suppose f : α) exists and is continuous on [0, 1] d for every α ∈{0, 1} d , and (ii) f (β) exists and is continuous on T (α) for every β ∈ J α , for every α ∈{0, 1} d \{0}.
Then, f ∈ F d,d ∞-mars and
f (β) .
Furthermore, if f (α) = 0 for all α ∈{0, 1} d with |α| >s in addition, then f ∈ F d,s ∞-mars and
f (β) .
Note that the integrals on the right-hand side of ( 35)and( 36) are only over those coordinates x l for which β l = max j β j (the remaining coordinates are set to zero).
The formula ( 35) is a multivariate generalization of the univariate formula (30), stating that for sufficiently smooth functions f , V mars (f ) is the sum of the L 1 norms of the mixed partial derivatives of f of order 2, where we take at most two partial derivatives along each coordinate and exactly two partial derivatives along at least one coordinate. Note that mixed partial derivatives of total order (β 1 +•••+β d )upto2d appear in (35). From this perspective, we can think of the smoothness order of our complexity measure V mars (•) as 2d,w h i c hi s proportional to the dimension d. This gives an intuitive explanation on why our estimators achieve the dimension-free rate n -4/5 (up to the logarithmic multiplicative factors), as we observedinSection3. However, we should also note that the maximum total order 2d is solely achieved by the mixed partial derivative f (2,...,2) , which prevents our function class from being too small and restrictive. For these reasons, we can say that our complexity measure is an effective constraint that leads to estimators avoiding the curse of dimensionality (to some extent) while keeping the corresponding function class reasonably large.
There is also an interesting connection between V mars (•) and V HK0 (•) via (35). Specifically, if the condition max j β j = 2i n( 35) is replaced with max j β j = 1, then one obtains the formula for V HK0 (f ) for sufficiently smooth functions f (see Lemma 7.4 of the Supplementary Material (Ki, Fang and Guntuboyina (2024))). Hence, we can also view V mars (•) as second-order Hardy-Krause variation (anchored at 0).
In Section 8 of the Supplementary Material (Ki, Fang and Guntuboyina ( 2024)), we describe the results of this section by specializing to the case d = s = 2. We encourage readers who find the results in this section too abstract to refer to Section 8 of Ki, Fang and Guntuboyina (2024) for more explicit formulae.
Here are some connections between our paper and existing works on nonparametric regression.
As mentioned earlier, our method can be viewed as a multivariate generalization of the piecewise linear locally adaptive regression spline estimator of Mammen and van de Geer (1997)( s e ea l s oSteidl, Didas and Neumann (2006)). There are other ways of generalizing the piecewise linear locally adaptive regression splines estimator to the multivariate setting as well (see, e.g., Parhi andNowak (2021, 2023) for some recent work). Also, we have seen that our method can be considered as a multivariate extension of first-order trend filtering (see, e.g., Kim et al. (2009), Tibshirani (2014)). Ortelli and van de Geer (2021)a n dSadhanala et al. ( 2021) recently studied different multivariate extensions of trend filtering. Although they covered all orders of trend filtering in contrast to our method, their methods were however restricted to lattice designs. Moreover, they imposed weaker penalties on their models, which resulted in their estimators converging to the true underlying function at dimensiondependent rates. In Section 9 of the Supplementary Material (Ki, Fang and Guntuboyina (2024)), we describe the method of Ortelli and van de Geer (2021) and compare their estimator to the discrete formation of our estimator in the equally spaced lattice design setting (see Remark 9.4 of Ki, Fang and Guntuboyina (2024) for more details).
We have also mentioned that V mars (•) can be viewed as second-order Hardy-Krause variation (anchored at 0). Our estimation strategy can thus be seen as second-order Hardy-Krause variation denoising. In Fang, Guntuboyina and Sen (2021), first-order Hardy-Krause variation denoising (i.e., squares estimation over functions with bounded Hardy-Krause variation) was studied. First-order Hardy-Krause variation denoising leads to piecewise constant fits while our method leads to MARS fits (linear combinations of products of ReLU functions of individual variables). Fang, Guntuboyina and Sen (2021) also proved that their estimator achieves a dimension-free (up to a logarithmic multiplicative factor) rate of convergence. However, it should be noted that their result is only proved in the fixed lattice design setting. Moreover, unlike our method, interaction order restriction is not considered in Fang, Guntuboyina and Sen (2021).
As is clear from the form of our functions (3) and our complexity measure (6), our method can also be considered as a multivariate ANOVA modeling method based on total variation constraints. There are a few works that utilize total variation penalties in multivariate ANOVA modeling. Petersen, Witten and Simon (2016), Yang and Tan (2018), and Sadhanala and Tibshirani (2019) utilized total variation of univatiate functions in additive modeling, which can be seen as a special case of ANOVA modeling where the interaction between covariates is not allowed. Also, Yang and Tan (2021) used for multivariate ANOVA modeling a class of penalties characterized in terms of certain hierarchical notions of total variation. Their hierarchical total variations are defined using a pre-fixed grid of points. Interestingly, for functions f that are sufficiently smooth, one of their hierarchical total variations (corresponding to m = 2in their notation) converges to V mars (f ) as the grid resolution becomes arbitrarily small.
It should be mentioned that Lin (1998Lin ( , 2000) ) also studied a multivariate ANOVA modeling method, but instead of L 1 norms as in our paper, they worked with penalties that are related to the squared L 2 -Sobolev norms. Relevance of their works to our paper is therefore not from the type of penalties but from tensor product structures on their basis functions. As our basis functions (1) are the tensor products of univariate ReLU functions, their basis functions are also the tensor products of univariate functions whose smoothness is constrained by the L 2 -Sobolev norms. It is notable that the multivariate function spaces considered in Lin (1998Lin ( , 2000) ) are defined as an appropriate completion of the pre-Hilbert space given by the tensor product of the univariate L 2 -Sobolev spaces. We are curious whether we can also view our function classes (e.g., F d,d ∞-mars ) as an appropriate completion of a tensor product space. However, it is unclear to us at this point what norm should be chosen for completion as a counterpart of the L 2 -Sobolev norms of Lin (1998Lin ( , 2000)). We believe this is an interesting direction to extend our work, which can provide a new perspective on our function spaces.
In addition, van der Laan, Benkeser and Cai (2023) discussed function classes similar to F d,d ∞-mars and norms similar to V mars (•) and used them for estimation in settings that are different from our classical nonparametric regression framework.
In this section, we provide the results of some numerical experiments illustrating the performance of our estimators from either the original or the approximate method. The performance of our estimators is compared to the performance of the usual MARS estimator in our experiments. The results for simulated data are presented first, and those for real data follow next. Our methods are implemented in the R package regmdc, which is available at https://github.com/DohyeongKi/regmdc. Our R package regmdc employs the R package Rmosek (based on interior point convex optimization) to solve the finite-dimensional lasso problem (13). For the usual MARS estimator, we use the R package earth based on Friedman (1991)andFriedman (1993). 6.1. Simulation studies. Mammen andvan de Geer (1997)andTibshirani (2014) demonstrated that their locally adaptive regression splines and trend filtering excel at estimating functions with locally varying smoothness. Considering that our estimator is a multivariate generalization of the piecewise (second-order) locally adaptive regression splines estimator and the first-order trend filtering estimator, it is natural to examine whether our methods also excel at adapting to the local variation of functions. To test the local adaptivity of our methods and compare it to the local adaptivity of the usual MARS method, we exploit in our simulation studies the following four functions (Function 1, Function 2, Function 3, and Function 4) whose smoothness varies significantly over the domain.
The definitions of the four functions (Function 1, Function 2, Function 3, and Function 4) are presented below. For each function, we consider the uniform design where X (i) are uniformly distributed on [0, 1] d and y i are generated according to the model (20). For Function 1 and Function 3, we also consider the equally spaced lattice design where
for some integers n k ≥ 2. In all cases, independent Gaussian errors with standard deviation 1 are added to function evaluations. Simulations are performed at various sample sizes for each function as presented in Table 1.
• Function 1 (L1). The first function f * :[0, 1] 2 → R is defined by
This function represents a two-dimensional damped sinusoidal wave. • Function 2 (L2). The second function f * :[0, 1] 5 → R is defined as in Function 1 but additionally includes three dummy variables x 3 , x 4 ,andx 5 .
This function is a (scaled) two-dimensional version of the Doppler function used for simulation studies in MammenandvandeGeer(1997)andTibshirani (2014). We add 0.001 to avoid division by zero. • Function 4 (L4). The fourth function f * :[0, 1] 5 → R is defined as in Function 3 but additionally includes three dummy variables x 3 , x 4 ,andx 5 .
We also use for comparison the following four Friedman's functions (see Friedman (1991), Section 4.3 and 4.4), which have been frequently utilized to measure the performance of nonparametric regression methods (see, e.g., Meyer, Leisch and Hornik (2003), Potts andSchmischke (2021, 2022)). The four functions (Function 5, Function 6, Function 7, and Function 8) are defined as below, and for each function, we generate data (X (1) ,y 1 ),...,(X (n) ,y n ) according to the model ( 20)whereX (i) are uniformly distributed on [0, 1] d . For Function 5 and Function 6, we consider Gaussian errors with standard deviation 1 as in Friedman (1991), Section 4.3. For Function 7 and Function 8, we consider Gaussian errors with standard deviation 125 and 0.1, which yield a 3:1 signal to noise ratio as designed in Friedman (1991), Section 4.4 (see also Meyer, Leisch and Hornik (2003)). Also, for each function, we conduct experiments on three different sample sizes as in Friedman (1991): 50, 100, and 200 for Function 5 and Function 6; 100, 200, and 400 for Function 7 and Function 8.
• Function 5 (F1). The first function f * :[0, 1] 10 → R is defined by f * (x 1 ,...,x 10 ) = 10 sin(πx 1 x 2 ) + 20(x 3 -0.5) 2 + 10x 4 + 5x 5 for x 1 ,...,x 10 ∈[0, 1].
• Function 6 (F2). The second function f * :[0, 1] 5 → R is defined as in Function 5, but here the five variables x 6 ,...,x 10 are removed.
w h e r eT 1 , T 2 ,a n dT 4 are the linear maps defined by T
whereT 1 , T 2 ,andT 4 are defined as in Function 7.
In all simulations, we compare the estimator f d,2 n,V and its approximate version f d,2 n,V to the usual MARS estimator whose order of interaction is restricted to 2. The approximate method is considered only when explanatory variables are uniformly distributed. The positive integers N k for the approximate method is set to 25 for each coordinate. Also, the tuning parameter V is selected by 10-fold cross-validation in all cases. For each function, we repeat the above data generation and estimator construction processes 25 times and average out computed losses. We use the squared empirical L 2 norm (17) for the fixed lattice design, and for the uniform design, we generate 1000 new samples and approximate the prediction error (22).
Table 1 presents the average loss of our estimators and the usual MARS estimator over 25 repetitions for each function. In Table 1, we can first observe that our estimators outperform the usual MARS estimator in capturing the local variation of the first four functions, L1, L2, L3, and L4, under both the lattice and the uniform design. The only exceptions are when we estimate the function L2 with 100 samples and the function L4 with 100, 200, and 400 samples. Even for these functions, our methods start performing better when they are given more samples. It turns out that our methods usually benefit more from increases in sample sizes compared to the usual MARS method. By comparing L1 with L2 and L3 with L4, we can also see that our methods however suffer more from the addition of dummy variables.
Table 1 also shows that our methods estimate the Friedman's functions F2, F3, and F4 much better, while the usual MARS method do better in the estimation of the function F1. However, the performance gap in estimating F1 narrows as the sample size increases, which reaffirms that increases in sample sizes are often more favorable for our methods than the usual MARS method. Also, comparing the results of F2 with those of F1, again we can see that our methods can benefit more from removing unnecessary covariates in advance. This hints that having appropriate variable selection as a pre-processing step can greatly improve the performance of our methods. We think it can be a promising avenue for future work.
Moreover, we can observe from Table 1 that the original and approximate method yield almost the same outputs in most cases. Only for the Doppler function L3, which has extremely wild fluctuation around the origin, the two methods show a notable difference. This tells us there may not be much degradation in practice from using the approximate method in place of the original one, and there even can be some gains, as we observed in the case of L3. From these simulation results, we can expect that although the usual MARS estimator might perform better on small sample sizes, our estimators usually outdo the usual MARS estimator when given enough samples. At what sample sizes such a transition occurs will definitely vary across functions, but in the above examples, they were reasonably small. 6.2. Real datasets. Here we use a few standard real datasets (Earnings, Airfoil Self-Noise, Abalone, Concrete, Ozone, Red Wine, and White Wine dataset) to compare our estimators with the usual MARS estimator. Brief descriptions of each dataset is presented in Section 10 of the Supplementary Material (Ki, Fang and Guntuboyina (2024)). For each dataset, we first linearly transform each explanatory variable into [0, 1].Ingeneral, there can be multiple options for linear transformations. If the domain of an explanatory variable is known as [m, M], then it is natural to subtract m from the variable and then divide it by Mm. In case the two extreme values of the domain, m and M, are unknown, we can consider the same linear transformation after simply setting m as the minimum and M as the maximum among observed values. Here we choose the latter option for all datasets.
We also split each dataset into a training set and a test set and use the mean squared error on the test set for comparison. We use 80% observations as training data and the remaining 20% observations as test data. Also, for every dataset, we focus on interactions between explanatory variables of order up to 2. In other words, we use f d,2 n,V or its approximate version f d,2 n,V for estimating regression functions and compare it to the usual MARS estimator whose order of interaction is restricted to 2. For the Earnings and the Airfoil Self-Noise datasets, we employ the original method, and for the other datasets, we employ the approximate method, with setting each N k to 25. Furthermore, the tuning parameter V is chosen by 10-fold crossvalidation in all cases.
In Table 2, we present for each dataset the number of explanatory variables, the number of data, and the type of our method we use (whether the original or the approximate method). Table 3 shows the average mean squared error of our estimator and the usual MARS estimator over 25 random training and test set splits for each dataset. In Table 3, we can see that our method outperforms the usual MARS method in almost all examples, while showing great improvement on the Airfoil Self-Noise, Concrete, Ozone, and White Wine dataset. The Red Wine dataset is the only exception for which the usual MARS method produces better or equally good fits. However, for its sibling (the White Wine dataset), which has about three times as many observations, our method is clearly a better option than the usual MARS method. This again proves that the usual MARS method may beat our method on small sample sizes, but our method reclaims the throne once enough data are provided. We can conclude from all the results in this section that not only can our estimators be an alternative to the usual MARS estimator, but also they can supplement each other.