# MARS via LASSO

Dohyeong<EMAIL_ADDRESS>[ Billy
<EMAIL_ADDRESS>[ Adityanand
<EMAIL_ADDRESS>[ University of California at
Berkeley and Google  University of California Berkeley
Evans Hall, Berkeley, CA 94720

E-mail:  Google LLC
Mountain View, CA 94043

E-mail:  University of California Berkeley
Evans Hall, Berkeley, CA 94720

E-mail:

###### Abstract

MARS is a popular method for nonparametric regression introduced by Friedman
in 1991. MARS fits simple nonlinear and non-additive functions to regression
data. We propose and study a natural LASSO variant of the MARS method. Our
method is based on least squares estimation over a convex class of functions
obtained by considering infinite-dimensional linear combinations of functions
in the MARS basis and imposing a variation based complexity constraint. We
show that our estimator can be computed via finite-dimensional convex
optimization and that it is naturally connected to nonparametric function
estimation techniques based on smoothness constraints. Under a simple design
assumption, we prove that our estimator achieves a rate of convergence that
depends only logarithmically on dimension and thus avoids the usual curse of
dimensionality to some extent. We implement our method with a cross-validation
scheme for the selection of the involved tuning parameter and show that it has
favorable performance compared to the usual MARS method in simulation and real
data settings.

62G08,

Constrained least squares estimation,

curse of dimensionality,

Hardy-Krause variation,

infinite-dimensional optimization,

integrated Brownian sheet,

locally adaptive regression spline,

L1 penalty,

metric entropy bounds,

mixed derivatives,

nonparametric regression,

piecewise linear function estimation,

small ball probability,

tensor products,

total variation regularization,

trend filtering,

###### keywords:

[class=MSC]

###### keywords:

, , and

t1Supported by NSF Grant DMS-2023505 t3Supported by NSF CAREER Grant
DMS-1654589

##  1 Introduction

We study a natural LASSO variant of the Multivariate Adaptive Regression
Splines (MARS) method (see Friedman [16] or Hastie et al. [20, Section 9.4])
for nonparametric regression. To understand the relation between a response
variable $y$ and $d$ explanatory variables $x_{1},\dots,x_{d}$ based on
observed data $(x^{(1)},y_{1}),\dots,(x^{(n)},y_{n})$ with
$x^{(i)}\in\mathbb{R}^{d}$ and $y_{i}\in\mathbb{R}$, MARS fits a function
$y=\hat{f}_{\text{mars}}(x_{1},\dots,x_{d})$ where $\hat{f}_{\text{mars}}$ is
a sparse linear combination of the functions of the form:

$\prod_{j=1}^{d}(b_{j}(x_{j}))^{\alpha_{j}}=\prod_{j:\alpha_{j}=1}b_{j}(x_{j})$
(1)

with $\alpha=(\alpha_{1},\dots,\alpha_{d})\in\\{0,1\\}^{d}$ and

$b_{j}(x_{j})=(x_{j}-t_{j})_{+}\text{ or }(t_{j}-x_{j})_{+}\text{ for some
real number }t_{j}.$

As a concrete example, to understand the relation between Weekly Earnings
($y$) and the two variables, Years of Experience ($x_{1}$) and Years of
Education ($x_{2}$), from a standard dataset ($\mathsf{ex1029}$ from the R
library $\mathsf{Sleuth3}$) collected from 25,437 full-time male workers in
1987, the default implementation of MARS fits the function

$\displaystyle
389.12-15.257(22-x_{1})_{+}+57.435(x_{2}-6)_{+}+47.947(x_{2}-13)_{+}$
$\displaystyle\quad-6.2336(8-x_{1})_{+}(x_{2}-13)_{+}-1.0916(x_{1}-8)_{+}(x_{2}-13)_{+}$
$\displaystyle\quad-0.47938(40-x_{1})_{+}(x_{2}-6)_{+}-2.7663(x_{1}-40)_{+}(x_{2}-6)_{+},$

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_{j}-t_{j})_{+}$ and $(t_{j}-x_{j})_{+}$. Further, MARS fits non-additive
functions because of the presence of products in (1), which enables
interactions between the explanatory variables $x_{1},\dots,x_{d}$. 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. Once a
reasonably large number of basis functions are added, a backward deletion
procedure is applied in order to remove superfluous basis functions. We refer
the reader to Hastie et al. [20, Section 9.4] for details.

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) and apply the
LASSO method of Tibshirani [36]. 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},\dots,x_{d}$ all take values in the interval $[0,1]$. In
other words, the domain of the regression function is assumed to be
$[0,1]^{d}$. In a practical setting, this can be achieved by subtracting the
minimum possible value and dividing by the range for each explanatory
variable. This scaling puts all the explanatory variables on a comparable
footing enabling the application of LASSO. Without such a scaling,
$x_{1},\dots,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_{j}-t_{j})_{+}$ would be unnatural (e.g., think
of the setting where $x_{1}$ is years of experience and $x_{2}$ is days of
education). 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 7.2 for two real data
applications).

In the rest of the paper, we assume that the explanatory variables
$x_{1},\dots,x_{d}$ all belong to $[0,1]$. The observed data is
$(x^{(1)},y_{1}),\dots,(x^{(n)},y_{n})$ where $x^{(i)}\in[0,1]^{d}$ and
$y_{i}\in\mathbb{R}$. We shall fit, to this data, functions of the form
$y=\hat{f}(x_{1},\dots,x_{d})$ where $\hat{f}:[0,1]^{d}\rightarrow\mathbb{R}$.
Our methodology is based on the application of LASSO with the MARS basis
functions. The restriction $x_{j}\in[0,1]$ allows us to make two
simplifications to the MARS setup:

  1. (i)

Instead of considering both kinds of functions $(x_{j}-t_{j})_{+}$ and
$(t_{j}-x_{j})_{+}$, we only take into account $(x_{j}-t_{j})_{+}$. This is
because, as each $x_{j}$ is assumed to be in $[0,1]$, we can write

$(t_{j}-x_{j})_{+}=(x_{j}-t_{j})_{+}-x_{j}+t_{j}=(x_{j}-t_{j})_{+}-(x_{j}-0)_{+}+t_{j},$

which implies that every linear combination of the functions of the form (1)
can also be written as a linear combination of the functions of the same form
(1) where $b_{j}(x_{j})=(x_{j}-t_{j})_{+}$ for some $t_{j}$.

  2. (ii)

We assume that $t_{j}\in[0,1)$ for each $j$. This is because when $t_{j}\geq
1$, the function $(x_{j}-t_{j})_{+}=0$ as $x_{j}\in[0,1]$, and for $t_{j}<0$,
the function $(x_{j}-t_{j})_{+}=x_{j}-t_{j}=(x_{j}-0)_{+}-t_{j}$ is a linear
combination of $(x_{j}-0)_{+}$ and the constant function 1.

We shall use LASSO to fit a sparse linear combination of basis functions from
the set:

$\mathcal{B}_{\text{mars}}^{d}:=\bigg{\\{}(x_{1},\dots,x_{d})\mapsto\prod_{j\in
S(\alpha)}(x_{j}-t_{j})_{+}:\alpha\in\\{0,1\\}^{d}\text{ and
}t_{j}\in[0,1)\text{ for }j\in S(\alpha)\bigg{\\}}.$ (2)

Here and in the rest of the paper, we use the notation

$S(\alpha):=\\{j\in[d]:\alpha_{j}=1\\},\qquad\text{where
$[d]:=\\{1,\cdots,d\\}$.}$

Because $\mathcal{B}_{\text{mars}}^{d}$ contains an uncountable number of
basis functions, we work with an infinite-dimensional version of LASSO.
Infinite-dimensional LASSO formulations have been used in many papers
including Rosset et al. [31], Bredies and Pikkarainen [7], Condat [11], Unser
et al. [39], Flinth and Weiss [15], Boyer et al. [5], Bredies and Carioni [6],
and Unser [38]. The main idea behind infinite-dimensional LASSO estimators is
to consider infinite linear combinations of basis functions that are
parametrized by finite signed measures and to measure complexity in terms of
the variations of the involved signed measures. In the MARS context, by
infinite linear combinations of basis functions in
$\mathcal{B}_{\text{mars}}^{d}$, we mean the functions of the form

$f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}(x_{1},\dots,x_{d}):=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{[0,1)^{|\alpha|}}\prod_{j\in
S(\alpha)}(x_{j}-t_{j})_{+}\,d\nu_{\alpha}(t^{(\alpha)}),$ (3)

where $\mathbf{0}:=(0,\dots,0)$, $a_{\mathbf{0}}\in{\mathbb{R}}$,
$\nu_{\alpha}$ is a finite signed measure on $[0,1)^{|\alpha|}$ with

$|\alpha|:=|S(\alpha)|=\sum_{j=1}^{d}\mathbf{1}\\{\alpha_{j}=1\\},$

and $t^{(\alpha)}$ indicates the vector $(t_{j},j\in S(\alpha))$ for each
binary vector $\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$. We shall
denote the collection of all such functions
$f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}$ by
$\mathcal{F}_{\infty-\text{mars}}^{d}$. The subscript $\infty$ highlights the
fact that $\mathcal{F}_{\infty-\text{mars}}^{d}$ contains infinite linear
combinations of functions from $\mathcal{B}_{\text{mars}}^{d}$. The usual MARS
functions (which are finite linear combinations of functions in
$\mathcal{B}_{\text{mars}}^{d}$) are special cases of (3) corresponding to
discrete measures $\nu_{\alpha}$. Indeed, when each $\nu_{\alpha}$ is
supported on a finite set $\\{(t_{lj}^{(\alpha)},j\in
S(\alpha)):l=1,\dots,k_{\alpha}\\}$ with

$\nu_{\alpha}\big{(}\big{\\{}\big{(}t_{lj}^{(\alpha)},j\in
S(\alpha)\big{)}\big{\\}}\big{)}=b_{l}^{(\alpha)}\qquad\text{for
$l=1,\dots,k_{\alpha}$},$ (4)

the function $f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}$ becomes

$(x_{1},\dots,x_{d})\mapsto
a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\sum_{l=1}^{k_{\alpha}}b_{l}^{(\alpha)}\prod_{j\in
S(\alpha)}\big{(}x_{j}-t_{lj}^{(\alpha)}\big{)}_{+}.$ (5)

We shall denote the collection of all the finite linear combinations of the
form (5) by $\mathcal{F}_{\text{mars}}^{d}$. As a function class,
$\mathcal{F}_{\infty-\text{mars}}^{d}$ is strictly larger than
$\mathcal{F}_{\text{mars}}^{d}$.

Our infinite-dimensional LASSO estimator minimizes the least squares criterion
over
$f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}\in\mathcal{F}_{\infty-\text{mars}}^{d}$
with a constraint on the complexity of
$f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}$. Our complexity measure involves the
sum of the variations of the underlying signed measures $\\{\nu_{\alpha}\\}$
and is an infinite-dimensional analogue of the usual $L_{1}$ norm of the
coefficients used in finite-dimensional LASSO. Recall that, for a signed
measure $\nu$ on $\Omega$ and a measurable subset $E\subseteq\Omega$, the
variation of $\nu$ on $E$ is denoted by $|\nu|(E)$ and is defined as the
supremum of $\sum_{A\in\pi}|\nu(A)|$ over all partitions $\pi$ of $E$ into a
countable number of disjoint measurable subsets. Using the variation of the
involved signed measures, we define the complexity measure for a function
$f=f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}\in\mathcal{F}_{\infty-\text{mars}}^{d}$
by

$V_{\mathrm{mars}}(f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}):=\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}|\nu_{\alpha}|\big{(}[0,1)^{|\alpha|}\setminus\\{\mathbf{0}\\}\big{)}.$
(6)

We are excluding the point $\mathbf{0}=(0,\dots,0)$ in the variation of
$\nu_{\alpha}$ because we want to only penalize those basis functions in
$\mathcal{B}_{\text{mars}}^{d}$ that include at least one nonlinear term and
leave unpenalized basis functions that are products of linear functions (note
that $(x_{j}-t_{j})_{+}=x_{j}$ is linear when $t_{j}=0$ because
$x_{j}\in[0,1]$). In Section 8.1.1, we prove that this complexity measure is
well-defined by showing the uniqueness of the representation
$f=f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}$.

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 $\nu_{\alpha}$ is the discrete signed measure given
by (4), we have

$V_{\mathrm{mars}}(f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}})=\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\sum_{l=1}^{k_{\alpha}}|b_{l}^{(\alpha)}|\cdot\mathbf{1}\big{\\{}\big{(}t_{lj}^{(\alpha)},j\in
S(\alpha)\big{)}\neq(0,\dots,0)\big{\\}},$

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.

We are now ready to introduce our infinite-dimensional LASSO estimator for
MARS fitting. It is defined as

$\hat{f}^{d}_{n,V}\in\underset{f}{\mathop{\rm
argmin}}\bigg{\\{}\sum_{i=1}^{n}\big{(}y_{i}-f(x^{(i)})\big{)}^{2}:f\in\mathcal{F}_{\infty-\text{mars}}^{d}\text{
and }V_{\mathrm{mars}}(f)\leq V\bigg{\\}}$ (7)

for a tuning parameter $V>0$. The study of $\hat{f}^{d}_{n,V}$ is the goal of
this paper, and we make the following contributions towards this goal.

Contribution 1: We provide an alternative characterization of
$\mathcal{F}_{\infty-\text{mars}}^{d}$ and $V_{\mathrm{mars}}(f)$ in terms of
the Hardy-Krause variations of mixed partial derivatives of $f$ of orders up
to 1. Hardy-Krause variation [1, 28] is a multivariate generalization of total
variation of univariate functions (we review basic properties of Hardy-Krause
variation in Section 2.2). The precise form of the connection between
$V_{\mathrm{mars}}(\cdot)$ and Hardy-Krause variation is given in Proposition
2.5. Here let us just additionally note that when $f$ is sufficiently smooth
(in the sense of Lemma 2.6), we have

$V_{\mathrm{mars}}(f)=\sum_{\begin{subarray}{c}\beta\in\\{0,1,2\\}^{d}\\\
\max_{j}\beta_{j}=2\end{subarray}}\int_{\bar{T}^{(\beta)}}|f^{(\beta)}|$ (8)

where $f^{(\beta)}$ is the partial derivative

$f^{(\beta)}:=\frac{\partial^{\beta_{1}+\dots+\beta_{d}}f}{\partial
x_{1}^{\beta_{1}}\dots\partial x_{d}^{\beta_{d}}}$ (9)

and $\bar{T}^{(\beta)}:=\bar{T}_{1}\times\dots\times\bar{T}_{d}$ with
$\bar{T}_{l}=[0,1]$ if $\beta_{l}=\max_{j}\beta_{j}$ and $\bar{T}_{l}=\\{0\\}$
if $\beta_{l}<\max_{j}\beta_{j}$ for each $l\in[d]$. The integral in (8) is
over only those coordinates $x_{l}$ for which $\bar{T}_{l}=[0,1]$, i.e.,
$\beta_{l}=\max_{j}\beta_{j}$. Also, only those vectors $\beta$ with
$\max_{j}\beta_{j}=2$ appear in the sum in (8). This means that
$V_{\mathrm{mars}}(f)$ is the sum of the $L_{1}$ norms of the partial
derivatives of $f$, where in each partial derivative, we take at most two
derivatives with respect to every variable and exactly two derivatives with
respect to at least one variable. As a consequence, $\hat{f}_{n,V}^{d}$ can be
seen as a second-order smoothness constrained least squares estimator that
minimizes least squares over the smoothness constraint that bounds the sum of
the $L_{1}$ norms of the mixed partial derivatives of orders up to 2.

In the univariate ($d=1$) case, (8) reduces to
$V_{\mathrm{mars}}(f)=\int_{0}^{1}|f^{\prime\prime}|$, and thus
$\hat{f}_{n,V}^{1}$ can be thought of as a constrained analogue of the locally
adaptive regression splines estimator of Mammen and van de Geer [25] when the
order $k$ (in the notation of [25]) equals 2. Thus, for $d\geq 2$,
$\hat{f}_{n,V}^{d}$ can be seen as a multivariate generalization of this
univariate estimator of [25]. Further, if the condition $\max_{j}\beta_{j}=2$
in (8) is replaced with $\max_{j}\beta_{j}=1$, then one obtains the Hardy-
Krause variation of $f$. As a result, we can refer to $V_{\mathrm{mars}}(f)$
as a second-order Hardy-Krause variation anchored at $(0,\dots,0)$ and the
estimation procedure leading to $\hat{f}_{n,V}^{d}$ as second-order Hardy-
Krause variation denoising. For comparison, a first-order Hardy-Krause
variation denoising estimator was studied in Fang et al. [14].

Contribution 2: We prove that $\hat{f}_{n,V}^{d}$ exists and can be computed
by applying finite-dimensional LASSO algorithms to the finite basis of
functions

${\mathcal{B}}^{d}_{\text{design}:\text{mars}}=\bigg{\\{}\prod_{j\in
S(\alpha)}(x_{j}-t_{j})_{+}:\alpha\in\\{0,1\\}^{d},\
t_{j}\in\\{0\\}\cup\big{\\{}x^{(i)}_{j}:i\in[n]\big{\\}}\text{ for }j\in
S(\alpha)\bigg{\\}},$ (10)

where we use the notation $x^{(i)}=(x^{(i)}_{1},\dots,x^{(i)}_{d})$ for the
$i^{th}$ design point $x^{(i)}$.
${\mathcal{B}}^{d}_{\text{design}:\text{mars}}$ is clearly a finite set of
functions, and the subscript in
${\mathcal{B}}^{d}_{\text{design}:\text{mars}}$ indicates the dependence on
the design points $x^{(1)},\dots,x^{(n)}$. Recall that finite-dimensional
LASSO estimation procedure usually zeros out many regression coefficients.
Thus, even though our estimator is defined as a solution to an optimization
problem over $\mathcal{F}_{\infty-\text{mars}}^{d}$, we can find
$\hat{f}_{n,V}^{d}$ that is a sparse linear combination of functions from
${\mathcal{B}}^{d}_{\text{design}:\text{mars}}$ so that
$\hat{f}_{n,V}^{d}\in\mathcal{F}_{\text{mars}}^{d}$. It is interesting to note
that the usual MARS algorithm also works with the finite basis
${\mathcal{B}}^{d}_{\text{design}:\text{mars}}$ instead of the full basis
$\mathcal{B}_{\text{mars}}^{d}$ although typically no theoretical
justification is provided for this reduction. Our estimation procedure can
therefore be seen as an alternative to the usual MARS procedure.

Contribution 3: We study the theoretical accuracy of $\hat{f}_{n,V}^{d}$ as an
estimator for an unknown regression function under a standard set of
assumptions. We assume that the data $(x^{(1)},y_{1}),\dots,(x^{(n)},y_{n})$
are generated according to the model

$y_{i}=f^{*}(x^{(i)})+\xi_{i}\qquad\text{where
$\xi_{i}\overset{\text{i.i.d}}{\sim}N(0,\sigma^{2})$}$ (11)

for an unknown regression function $f^{*}:[0,1]^{d}\rightarrow{\mathbb{R}}$
belonging to the class $\mathcal{F}_{\infty-\text{mars}}^{d}$. We also assume
that the design points $x^{(1)},\dots,x^{(n)}$ form an equally spaced lattice
in $[0,1]^{d}$. In this fixed-design setting, we measure the accuracy of
$\hat{f}_{n,V}^{d}$ with the loss function:

$\mathcal{L}\big{(}\hat{f}_{n,V}^{d},f^{*}\big{)}:=\frac{1}{n}\sum_{i=1}^{n}\big{(}\hat{f}_{n,V}^{d}(x^{(i)})-f^{*}(x^{(i)})\big{)}^{2}.$

Under these assumptions, we prove in Theorem 5.1 that

${\mathbb{E}}\mathcal{L}\big{(}\hat{f}_{n,V}^{d},f^{*}\big{)}=O\big{(}n^{-\frac{4}{5}}(\log
n)^{\frac{3(2d-1)}{5}}\big{)}$ (12)

provided the tuning parameter $V$ in our estimator satisfies $V\geq
V_{\mathrm{mars}}(f^{*})$. The multiplicative constant underlying $O(\cdot)$
in (12) depends on $d,\sigma$, and $V$ (see Theorem 5.1 for the explicit
expression). This result is interesting from the following perspective.
Already in the univariate case ($d=1$), the estimator $\hat{f}_{n,V}^{1}$
achieves the rate $n^{-4/5}$ (see, e.g., Mammen and van de Geer [25, Theorem
10] or Guntuboyina et al. [19, Theorem 2.1]). Thus, (12) implies that in going
from the univariate to the multivariate setting, the rate only deteriorates by
a logarithmic multiplicative factor, which suggests that our LASSO method for
MARS fitting can be an effective function estimation technique in higher
dimensions. The assumption that the design points $x^{(1)},\dots,x^{(n)}$ form
an equally spaced lattice in $[0,1]^{d}$ is admittedly restrictive, although
this assumption is standard in theoretical analysis in nonparametric function
estimation (see, e.g., Nemirovski [26]), and it is also directly relevant to
applications involving image denoising. However, we are currently unable to
handle more general design assumptions.

The key ingredient behind our proof of (12) is the following metric entropy
upper bound. For $m\geq 1$, let ${\mathcal{D}}_{m}$ denote the class of all
the functions of the form

$(x_{1},\dots,x_{m})\mapsto\int(x_{1}-t_{1})_{+}\dots(x_{m}-t_{m})_{+}\,d\nu(t)$

where $\nu$ is a signed measure on $[0,1]^{m}$ with the total variation
$|\nu|([0,1]^{m})\leq 1$. We prove that there are two positive constants
$C_{m}$ and $\epsilon_{m}$ depending only on $m$ such that

$\log N(\epsilon,{\mathcal{D}}_{m},\|\cdot\|_{2})\leq
C_{m}\epsilon^{-\frac{1}{2}}\Big{(}\log\frac{1}{\epsilon}\Big{)}^{\frac{3(2m-1)}{4}}\qquad\text{for
every $0<\epsilon<\epsilon_{m}$}$ (13)

where $N(\epsilon,{\mathcal{D}}_{m},\|\cdot\|_{2})$ is the $\epsilon$-covering
number of ${\mathcal{D}}_{m}$ under the $L_{2}$ metric. Our proof of (13) is
based on establishing a connection between the metric entropy of
${\mathcal{D}}_{m}$ and the small ball probability of integrated Brownian
sheet (using ideas and results from Blei et al. [3], Gao [18], Li and Linde
[22], and Artstein et al. [2]) and then deriving the small ball probability of
integrated Brownian sheet (using ideas and results from Dunker et al. [13] and
Chen and Li [10]).

Contribution 4: Under the model assumption (11) and the equally spaced lattice
design assumption, we prove the minimax rate under the loss function
$\mathcal{L}$ over the class

$\big{\\{}f\in\mathcal{F}_{\infty-\text{mars}}^{d}:V_{\mathrm{mars}}(f)\leq
V\big{\\}}$ (14)

is bounded from below by

$n^{-\frac{4}{5}}(\log n)^{\frac{4(d-1)}{5}}.$

This result implies that the logarithmic term in (12) cannot be completely
eliminated. In other words, every estimator has to pay a logarithmic price
that grows with the dimension $d$ in the worst case over the class (14).

Contribution 5: In practical implementations of the usual MARS method, the
considered set of functions is often further constrained by allowing in (5)
only those $\alpha$ with $|\alpha|\leq s$ for some fixed $s\in[d]$. Motivated
by this, we consider the class

$\mathcal{F}_{\infty-\text{mars}}^{d,s}:=\big{\\{}f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}\in\mathcal{F}_{\infty-\text{mars}}^{d}:\nu_{\alpha}=0\text{
whenever }|\alpha|>s\big{\\}}$ (15)

and the corresponding constrained least squares estimator:

$\hat{f}^{d,s}_{n,V}\in\underset{f}{\mathop{\rm
argmin}}\bigg{\\{}\sum_{i=1}^{n}\big{(}y_{i}-f(x^{(i)})\big{)}^{2}:f\in\mathcal{F}_{\infty-\text{mars}}^{d,s}\text{
and }V_{\mathrm{mars}}(f)\leq V\bigg{\\}}.$ (16)

Intuitively, $\mathcal{F}_{\infty-\text{mars}}^{d,s}$ is the subcollection of
$\mathcal{F}_{\infty-\text{mars}}^{d}$ consisting of all the functions whose
order of interaction in variables is restricted to $s$. Under the same
assumptions leading to (12) as well as the additional assumption
$f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{d,s}$, we prove that

${\mathbb{E}}\mathcal{L}\big{(}\hat{f}_{n,V}^{d,s},f^{*}\big{)}=O\big{(}n^{-\frac{4}{5}}(\log
n)^{\frac{3(2s-1)}{5}}\big{)}$ (17)

where the multiplicative constant underlying $O(\cdot)$ depends on $d,\sigma$,
and $V$. The logarithmic multiplicative factor in (17) now only depends on $s$
and not on $d$. This suggests that $\hat{f}_{n,V}^{d,s}$ with $s\ll d$ can be
an even more effective estimator than $\hat{f}_{n,V}^{d}$ when $d$ is large
and $f^{*}$ does not have higher-order interactions in variables. Note however
that the constant multiplicative factor underlying the rate in (17) still
depends on the dimension $d$.

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 show that it has favorable performance compared to the usual
MARS method in simulation and real data settings.

The rest of the paper is organized as follows. In Section 2, we provide
alternative characterizations of
$\mathcal{F}_{\infty-\text{mars}}^{d},V_{\mathrm{mars}}(\cdot)$, and
$\hat{f}_{n,V}^{d}$ in terms of smoothness. Here we also discuss connections
between our estimator and some existing estimators from the literature. In
Section 3, we present our results on the existence of $\hat{f}_{n,V}^{d}$ and
its computation via reduction to a finite-dimensional LASSO optimization
problem. In Section 4, we consider the special case where the design points
$x^{(1)},\dots,x^{(n)}$ form an equally spaced lattice in $[0,1]^{d}$, and for
this design assumption, we provide another finite-dimensional optimization
problem for the computation of $\hat{f}_{n,V}^{d}$. In Section 5, we present
our risk results for the equally spaced lattice design. Section 6 is dedicated
to the study of the estimator $\hat{f}_{n,V}^{d,s}$. In Section 7, we
illustrate the performance of our method in simulated and real data settings
and compare its performance to that of the usual MARS algorithm. Proofs of all
our results can be found in Section 8.

##  2 Characterization in terms of Smoothness

In this section, we provide alternative characterizations of
$\mathcal{F}_{\infty-\text{mars}}^{d},V_{\mathrm{mars}}(\cdot)$, and
$\hat{f}^{d}_{n,V}$ in terms of smoothness. To motivate the results for
general $d\geq 1$, let us first consider the univariate case $d=1$.

###  2.1 Smoothness Characterization for $d=1$

For $d=1$, $\mathcal{F}_{\infty-\text{mars}}^{1}$ consists of all the
functions $f:[0,1]\rightarrow\mathbb{R}$ of the form

$f(x)=a_{0}+\int_{[0,1)}(x-t)_{+}\,d\nu(t)$ (18)

for a real number $a_{0}$ and a finite signed measure $\nu$ on $[0,1)$, and
the complexity measure of such function $f$ is given by the variation of $\nu$
on $(0,1)$: $V_{\mathrm{mars}}(f)=|\nu|((0,1))$.

The following well-known argument shows that one can characterize
$\mathcal{F}_{\infty-\text{mars}}^{1}$ in terms of smoothness. Replacing
$(x-t)_{+}$ with $\int_{0}^{1}\mathbf{1}\\{t\leq s\leq x\\}ds$ in the integral
in (18) and changing the order of integration, we obtain

$f(x)=a_{0}+\int_{0}^{x}g(t)\,dt$ (19)

where the function $g:[0,1]\rightarrow\mathbb{R}$ is given by

$g(t)=\nu([0,t]\cap[0,1)).$ (20)

The equation (19) indicates that $f$ is absolutely continuous. Also, we can
note that the function $g$ from (20) is right-continuous on $[0,1]$ and left-
continuous at 1, and the total variation $V(g)$ of $g$ is finite and satisfies

$V(g)=|\nu|((0,1)).$ (21)

Here the total variation of the function $g:[0,1]\rightarrow\mathbb{R}$ is
defined by

$V(g):=\sup_{0=u_{0}<u_{1}<\dots<u_{k}=1}\sum_{i=0}^{k-1}\big{|}g(u_{i+1})-g(u_{i})\big{|}$

where the supremum is over all $k\geq 1$ and partitions
$0=u_{0}<u_{1}<\dots<u_{k}=1$ of $[0,1]$.

Conversely, every function $g:[0,1]\rightarrow\mathbb{R}$ that is right
continuous on $[0,1]$, left continuous at 1, and has finite total variation
can be written as (20) for a unique signed measure $\nu$ on $[0,1)$ (see,
e.g., Aistleitner and Dick [1, Theorem 3]). From these arguments, we can see
that $\mathcal{F}_{\infty-\text{mars}}^{1}$ has the following alternative
characterization:

$\displaystyle\begin{split}&\mathcal{F}_{\infty-\text{mars}}^{1}=\Big{\\{}f:[0,1]\rightarrow\mathbb{R}:\exists
a_{0}\in\mathbb{R}\text{ and }g:[0,1]\rightarrow\mathbb{R}\text{ s.t. }\\\
&\qquad\qquad\qquad\qquad\qquad g\text{ is right-continuous on }[0,1],\text{
left-continuous at 1, }\\\ &\qquad\qquad\qquad\qquad\qquad V(g)<\infty,\text{
and }f(x)=a_{0}+\int_{0}^{x}g(t)\,dt\text{ for all
}x\in[0,1]\Big{\\}}.\end{split}$ (22)

Also, we can see that the complexity measure $V_{\mathrm{mars}}(f)$ for
$f\in\mathcal{F}_{\infty-\text{mars}}^{1}$ is equal to the total variation
$V(g)$ of the function $g$ appearing in (19).

The function $g$ appearing in the representation (19) is identical to the
derivative $f^{\prime}$ if $f$ is sufficiently smooth. The representation (19)
implies, by the fundamental theorem of calculus (see, e.g., Rudin [32, Theorem
7.10]), that $f^{\prime}$ exists and is equal to $g$ almost everywhere (w.r.t
the Lebesgue measure) on $[0,1]$. Hence, if
$f\in\mathcal{F}_{\infty-\text{mars}}^{1}$ is smooth in the sense that
$f^{\prime}$ exists everywhere and continuous on $[0,1]$, then $g$ is
identical to $f^{\prime}$ and $V_{\mathrm{mars}}(f)=V(f^{\prime})$. For these
reasons, somewhat loosely, we can describe the estimator $\hat{f}_{n,V}^{1}$
as

$\hat{f}_{n,V}^{1}\in\mathop{\rm
argmin}_{f}\bigg{\\{}\sum_{i=1}^{n}\big{(}y_{i}-f(x^{(i)})\big{)}^{2}:V(f^{\prime})\leq
V\bigg{\\}}.$

The corresponding penalized version

$\mathop{\rm
argmin}_{f}\bigg{\\{}\sum_{i=1}^{n}\big{(}y_{i}-f(x^{(i)})\big{)}^{2}+\lambda
V(f^{\prime})\bigg{\\}}$ (23)

was first proposed by Mammen and van de Geer [25] as a part of a class of
estimators collectively called locally adaptive regression splines. Thus, in
the univariate case, $\hat{f}_{n,V}^{1}$ is a constrained analogue of the
locally adaptive regression spline estimator of [25] when the order $k$ (in
the notation of [25]) equals 2\. Steidl et al. [33] used the terminology
second-order total variation regularization, and Kim et al. [21] and
Tibshirani [37] used the terminology trend filtering for (23). For $d\geq 1$,
our estimator $\hat{f}_{n,V}^{d}$ is thus a multivariate generalization of
second-order total variation regularization or trend filtering.

From the alternative characterization of
$\mathcal{F}_{\infty-\text{mars}}^{1}$ given above, it follows that every
sufficiently smooth function $f:[0,1]\rightarrow\mathbb{R}$ belongs to
$\mathcal{F}_{\infty-\text{mars}}^{1}$. Indeed, if $f^{\prime}$ and
$f^{\prime\prime}$ exist everywhere and are continuous on $[0,1]$, then we
have

$f(x)=f(0)+\int_{0}^{x}f^{\prime}(t)\,dt$

for all $x\in[0,1]$, and

$V(f^{\prime})=\int_{0}^{1}|f^{\prime\prime}|<\infty.$

Thus, in this case, $f$ belongs to $\mathcal{F}_{\infty-\text{mars}}^{1}$ and
the complexity measure of $f$ can be written as

$V_{\mathrm{mars}}(f)=\int_{0}^{1}|f^{\prime\prime}|.$ (24)

The formula above highlights the role of the second derivative
$f^{\prime\prime}$ in the definition of $V_{\mathrm{mars}}(f)$ in the
univariate case.

###  2.2 Hardy-Krause Variation

As we have seen in the previous subsection, in the univariate case,
$\mathcal{F}_{\infty-\text{mars}}^{1}$ consists of all the functions $f$
satisfying (19) with some function $g$ having finite total variation. An
analogous characterization holds for all $d\geq 1$, and it will be presented
in the next subsection. For general $d$, the role of total variation in $d=1$
will be played by the Hardy-Krause variation, which is an extension of total
variation of univariate functions to higher dimensions. Thus, here we briefly
introduce Hardy-Krause variation and recall some properties of the variation
that we shall use in the following subsection. Standard references for Hardy-
Krause variation are Aistleitner and Dick [1] and Owen [28].

Hardy-Krause variation is based on another variation called Vitali variation.
We first recall Vitali variation and some related concepts and then define
Hardy-Krause variation using them. We shall use the following notation below.
For $m\geq 1$ and two $m$-dimensional vectors $u=(u_{1},\cdots,u_{m})$ and
$v=(v_{1},\cdots,v_{m})$ with $u_{i}\leq v_{i}$ for all $i\in[m]$, we let
$[u,v]$ be the axis-aligned closed rectangle:

$[u,v]:=\\{x\in\mathbb{R}^{m}:u_{i}\leq x_{i}\leq v_{i}\text{ for all
}i\in[m]\\}.$

The quasi-volume of a function $g$ on $[u,v]$ is defined by

$\Delta(g;[u,v])=\sum_{\delta\in\\{0,1\\}^{m}}(-1)^{\delta_{1}+\cdots+\delta_{m}}g(v_{1}+\delta_{1}(u_{1}-v_{1}),\cdots,v_{m}+\delta_{m}(u_{m}-v_{m})).$

We say that a collection $\mathcal{P}$ of subsets of $[0,1]^{m}$ is an axis-
aligned split if

$\mathcal{P}=\bigg{\\{}\prod_{k=1}^{m}\big{[}u^{(k)}_{l_{k}-1},u^{(k)}_{l_{k}}\big{]}:l_{k}\in[n_{k}]\text{
for }k\in[m]\bigg{\\}}$

for some $0=u^{(k)}_{0}<\cdots<u^{(k)}_{n_{k}}=1$ for $k\in[m]$. Note that
every element of $\mathcal{P}$ is an axis-aligned closed rectangle of the form

$\Big{[}\big{(}u^{(1)}_{l_{1}-1},\cdots,u^{(m)}_{l_{m}-1}\big{)},\big{(}u^{(1)}_{l_{1}},\cdots,u^{(m)}_{l_{m}}\big{)}\Big{]}$

for some $l_{k}\in[n_{k}]$ for $k\in[m]$.

We next define Vitali variation which coincides with the usual total variation
when $m=1$.

######  Definition 2.1 (Vitali variation).

For a real-valued function $g$ on $[0,1]^{m}$, we define the Vitali variation
of $g$ by

$V^{(m)}(g)=\sup_{\mathcal{P}}\sum_{A\in\mathcal{P}}|\Delta(g;A)|,$

where the supremum is taken over all axis-aligned splits $\mathcal{P}$ of
$[0,1]^{m}$.

For defining Hardy-Krause variation, we need the following notation. For each
$m$-dimensional vector $\beta$ with nonnegative integer components, we let

$\displaystyle\bar{T}^{(\beta)}:=\bar{T}_{1}\times\dots\times\bar{T}_{m}~{}\text{
where }~{}\bar{T}_{k}=\begin{cases}[0,1]&\mbox{if
}\beta_{k}=\max_{j}\beta_{j}\\\ \\{0\\}&\mbox{otherwise}.\end{cases}$ (25)

######  Definition 2.2 (Hardy-Krause variation).

For a real-valued function $g$ on $[0,1]^{m}$, we define the Hardy-Krause
variation of $g$ by

$V_{\mathrm{HK0}}(g)=\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}V^{(|\alpha|)}(g|_{\bar{T}^{(\alpha)}}).$
(26)

Here for each $\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}$,
$g|_{\bar{T}^{(\alpha)}}$ is the restriction of $g$ to $\bar{T}^{(\alpha)}$,
and $V^{(|\alpha|)}(g|_{\bar{T}^{(\alpha)}})$ is the Vitali variation of
$g|_{\bar{T}^{(\alpha)}}$ when treated as a function on $[0,1]^{|\alpha|}$.

When $m=1$, there is only one term in (26) and thus
$V_{\mathrm{HK0}}(g)=V^{(1)}(g)=V(g)$.

The sets $\bar{T}^{(\alpha)}$ for
$\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}$ in Definition 2.2 are the
faces of $[0,1]^{m}$ that are adjacent to the point $\mathbf{0}$. Thus, this
Hardy-Krause variation is often referred to as Hardy-Krause variation anchored
at $\mathbf{0}$. It is also possible to define Hardy-Krause variation using a
different anchor (see, e.g., Aistleitner and Dick [1]), but we only need
Hardy-Krause variation anchored at $\mathbf{0}$ for our purpose. To highlight
the special role of $\mathbf{0}$ in our Hardy-Krause variation, we abbreviate
it to $\mathrm{HK0}$ variation and use the notation $V_{\mathrm{HK0}}(\cdot)$.

There is a well-known close connection between functions with finite
$\mathrm{HK0}$ variation and cumulative distribution functions of finite
signed measures. The following result is a version of this connection that we
shall use. This result follows essentially from Aistleitner and Dick [1,
Theorem 3] with minor modification (we provide the argument for the
modification in Section 8.1.2).

######  Proposition 2.3.

Suppose a real-valued function $g$ defined on $[0,1]^{m}$ has finite
$\mathrm{HK0}$ variation and is coordinate-wise right-continuous on
$[0,1]^{m}$. Assume also that $g$ is coordinate-wise left-continuous at each
point $x\in[0,1]^{m}\setminus[0,1)^{m}$ with respect to all the
$j^{\text{th}}$ coordinates where $x_{j}=1$. Then, there exists a unique
finite signed measure $\nu$ on $[0,1)^{m}$ such that

$g(x)=\nu([\mathbf{0},x]\cap[0,1)^{m})\qquad\text{for every $x\in[0,1]^{m}$}.$
(27)

Also, if a real-valued function $g$ defined on $[0,1]^{m}$ is of the form (27)
for some finite signed measure $\nu$ on $[0,1)^{m}$, then $g$ has finite
$\mathrm{HK0}$ variation and

$|\nu|([0,1)^{m}\setminus\\{\mathbf{0}\\})=V_{\mathrm{HK0}}(g).$

Finally, let us note that when $g$ is sufficiently smooth, its Hardy-Krause
variation $V_{\mathrm{HK0}}(g)$ can be written as the sum of the $L_{1}$ norms
of mixed partial derivatives of $g$ of orders up to 1. The precise statement
is given in the following lemma (proved in Section 8.1.3). In the next lemma
and in the rest of the paper, we use the following notation for partial
derivatives. For an $m$-dimensional vector $\beta$ with nonnegative integer
components,

$g^{(\beta)}=\frac{\partial^{\beta_{1}+\dots+\beta_{m}}g}{\partial
x_{1}^{\beta_{1}}\dots\partial x_{m}^{\beta_{m}}}.$

######  Lemma 2.4.

Suppose a real-valued function $g$ defined on $[0,1]^{m}$ is smooth in the
sense that $g^{(\alpha)}$ exists, is irrespective of the order of
differentiation, and is continuous on $[0,1]^{m}$ for every
$\alpha\in\\{0,1\\}^{m}$. Then, $g$ has finite $\mathrm{HK0}$ variation and

$V_{\mathrm{HK0}}(g)=\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\bar{T}^{(\alpha)}}|g^{(\alpha)}|.$
(28)

###  2.3 Smoothness Characterization for all $d\geq 1$

The following result provides an alternative characterization of
$\mathcal{F}_{\infty-\text{mars}}^{d}$ and $V_{\mathrm{mars}}(\cdot)$ in terms
of smoothness. Proposition 2.3 plays a crucial role in deriving this result
(proved in Section 8.1.4). Recall that we use the notation $t^{(\alpha)}$ to
indicate the vector $(t_{j},j\in S(\alpha))$ for each
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$.

######  Proposition 2.5.

The function class $\mathcal{F}_{\infty-\text{mars}}^{d}$ consists precisely
of all the functions of the form

$f(x_{1},\cdots,x_{d})=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{j\in
S(\alpha)}[0,x_{j}]}g_{\alpha}(t^{(\alpha)})\,dt^{(\alpha)}$ (29)

for some $a_{\mathbf{0}}\in\mathbb{R}$ and for some collection of functions
$\\{g_{\alpha}:\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}\\}$, where for
each $\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$,

  1. (i)

$g_{\alpha}$ is a real-valued function on $[0,1]^{|\alpha|}$,

  2. (ii)

$V_{\mathrm{HK0}}(g_{\alpha})<\infty$,

  3. (iii)

$g_{\alpha}$ is coordinate-wise right-continuous on $[0,1]^{|\alpha|}$, and

  4. (iv)

$g_{\alpha}$ is coordinate-wise left-continuous at each point
$x^{(\alpha)}\in[0,1]^{|\alpha|}\setminus[0,1)^{|\alpha|}$ with respect to all
the $j^{\text{th}}$ coordinates where $x_{j}=1$.

Further, the variation of $f$ in (29) can be written in terms of the Hardy-
Krause variations of $g_{\alpha}$ as

$V_{\mathrm{mars}}(f)=\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}V_{\mathrm{HK0}}(g_{\alpha}).$
(30)

Proposition 2.5 is completely analogous to (22) for $d=1$. Specifically, the
condition (29) is analogous to the univariate absolute continuity condition
(19). The condition $V_{\mathrm{HK0}}(g_{\alpha})<\infty$ for each
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$ is analogous to the
univariate condition $V(g)<\infty$. Coordinate-wise right-continuity of each
$g_{\alpha}$ on $[0,1]^{|\alpha|}$ is analogous to univariate right-continuity
on $[0,1]$. Lastly, coordinate-wise left-continuity of each $g_{\alpha}$ at
each $x^{(\alpha)}\in[0,1]^{|\alpha|}\setminus[0,1)^{|\alpha|}$ (with respect
to all the $j^{\text{th}}$ coordinates where $x_{j}=1$) is analogous to
univariate left-continuity at 1. It is also interesting to note that
$V_{\mathrm{mars}}(f)$ equals the sum of the Hardy-Krause variations of
$g_{\alpha}$ over $\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$.

The condition (29) can be interpreted as a multivariate absolute continuity
condition that imposes a certain kind of smoothness on $f$ (note that there is
no appearance of the basis functions $(x_{j}-t_{j})_{+}$ nor any signed
measures in (29)). We have remarked that the function $g$ appearing in the
univariate absolute continuity condition (19) is identical to the derivative
$f^{\prime}$ if $f$ is sufficiently smooth. One can similarly match each
function $g_{\alpha}$ appearing in (29) with the partial derivative
$f^{(\alpha)}$ when $f$ is sufficiently smooth. Indeed, if $f$ is smooth in
the sense of Lemma 2.4, then it can be shown that for each
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$,

$g_{\alpha}(x^{(\alpha)})=f^{(\alpha)}(\tilde{x}^{(\alpha)})\qquad\text{for
every $x^{(\alpha)}\in[0,1]^{|\alpha|}$,}$

where $x^{(\alpha)}=(x_{j},j\in S(\alpha))$ and $\tilde{x}^{(\alpha)}$ is its
$d$-dimensional extension defined by

$\displaystyle\tilde{x}^{(\alpha)}_{j}=\begin{cases}x_{j}&\mbox{if }j\in
S(\alpha)\\\ 0&\mbox{otherwise}.\end{cases}$

Thus, in this case, we can write the equation (30) as

$V_{\mathrm{mars}}(f)=\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}V_{\mathrm{HK0}}(f^{(\alpha)}|_{\bar{T}^{(\alpha)}}),$

where $f^{(\alpha)}|_{\bar{T}^{(\alpha)}}$ is the restriction of
$f^{(\alpha)}$ to $\bar{T}^{(\alpha)}$ and
$V_{\mathrm{HK0}}(f^{(\alpha)}|_{\bar{T}^{(\alpha)}})$ refers to the Hardy-
Krause variation of $f^{(\alpha)}|_{\bar{T}^{(\alpha)}}$ viewed as a function
on $[0,1]^{|\alpha|}$. With this understanding, somewhat loosely, we can write
the estimator $\hat{f}_{n,V}^{d}$ as

$\hat{f}_{n,V}^{d}\in\mathop{\rm
argmin}_{f}\bigg{\\{}\sum_{i=1}^{n}\big{(}y_{i}-f(x^{(i)})\big{)}^{2}:\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}V_{\mathrm{HK0}}(f^{(\alpha)}|_{\bar{T}^{(\alpha)}})\leq
V\bigg{\\}}.$ (31)

Thus, our estimator $\hat{f}_{n,V}^{d}$ can be viewed as a least squares
estimator under a specific smoothness constraint involving the sum of the
Hardy-Krause variations of mixed partial derivatives of orders up to 1.

Proposition 2.5 implies that every sufficiently smooth function belongs to
$\mathcal{F}_{\infty-\text{mars}}^{d}$. This is proved in the next result
(proved in Section 8.1.5), which also gives an expression for
$V_{\mathrm{mars}}(f)$ in terms of the $L_{1}$ norms of mixed partial
derivatives of $f$ of orders up to 2 for such sufficiently smooth functions.
The following notation will be used below. For each
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$, we let

$J_{\alpha}:=\big{\\{}\beta\in\\{0,1,2\\}^{d}:\beta_{j}=0\text{ if
}\alpha_{j}=0,~{}\beta_{j}\in\\{1,2\\}\text{ if }\alpha_{j}=1,\text{ and
}\max_{j}\beta_{j}=2\big{\\}}.$ (32)

Also, recall the notation $\bar{T}^{(\beta)}$ from (25).

######  Lemma 2.6.

Suppose $f:[0,1]^{d}\rightarrow\mathbb{R}$ is smooth in the sense that

  1. (i)

$f^{(\alpha)}$ exists, is irrespective of the order of differentiation, and is
continuous on $[0,1]^{d}$ for every $\alpha\in\\{0,1\\}^{d}$, and

  2. (ii)

$f^{(\beta)}$ exists, is irrespective of the order of differentiation, and is
continuous on $\bar{T}^{(\alpha)}$ for every $\beta\in J_{\alpha}$, for every
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$.

Then, $f\in\mathcal{F}_{\infty-\text{mars}}^{d}$ and

$V_{\mathrm{mars}}(f)=\sum_{\begin{subarray}{c}\beta\in\\{0,1,2\\}^{d}\\\
\max_{j}\beta_{j}=2\end{subarray}}\int_{\bar{T}^{(\beta)}}|f^{(\beta)}|.$ (33)

It should be noted that the integral on the right-hand side of (33) is only
over those coordinates $x_{l}$ for which $\beta_{l}=\max_{j}\beta_{j}$ (the
remaining coordinates are set to zero).

The formula (33) states that $V_{\mathrm{mars}}(f)$ is the sum of the $L_{1}$
norms of the partial derivatives of $f$, where in each partial derivative, we
take at most two derivatives with respect to every variable and exactly two
derivatives with respect to at least one variable. This is a multivariate
generalization of the univariate formula (24) and highlights the role of the
second derivatives in the determination of $V_{\mathrm{mars}}(f)$. Note that
partial derivatives of total order up to $2d$ appear in (33) and more than two
derivatives are never taken with respect to any variable.

There is an interesting connection between (28) and (33). Specifically, if the
condition $\max_{j}\beta_{j}=2$ in (33) is replaced with
$\max_{j}\beta_{j}=1$, then the resulting formula is identical to (28), which
represents the Hardy-Krause variation of $f$. As a result, we can view
$V_{\mathrm{mars}}(f)$ as second-order Hardy-Krause variation.

In the following example, we specialize the results of this section to the
case $d=2$.

######  Example 2.7 ($d=2$).

According to the definition given in Section 1,
$\mathcal{F}_{\infty-\text{mars}}^{2}$ consists of all the functions of the
form

$\displaystyle f(x_{1},x_{2})$
$\displaystyle=a_{0,0}+\int_{[0,1)}(x_{1}-t_{1})_{+}\,d\nu_{1,0}(t_{1})+\int_{[0,1)}(x_{2}-t_{2})_{+}\,d\nu_{0,1}(t_{2})$
$\displaystyle\qquad+\int_{[0,1)^{2}}(x_{1}-t_{1})_{+}(x_{2}-t_{2})_{+}\,d\nu_{1,1}(t_{1},t_{2})$

for $a_{0,0}\in\mathbb{R}$ and signed measures $\nu_{1,0}$, $\nu_{0,1}$ and
$\nu_{1,1}$ (on $[0,1),[0,1)$ and $[0,1)^{2}$ respectively).

On the other hand, according to the characterization of Proposition 2.5,
$\mathcal{F}_{\infty-\text{mars}}^{2}$ consists of all functions of the form

$f(x_{1},x_{2})=a_{0,0}+\int_{0}^{x_{1}}g_{1,0}(t_{1})\,dt_{1}+\int_{0}^{x_{2}}g_{0,1}(t_{2})\,dt_{2}+\int_{0}^{x_{1}}\int_{0}^{x_{2}}g_{1,1}(t_{1},t_{2})\,dt_{2}\,dt_{1}$
(34)

for a real number $a_{0,0}$ and real-valued functions $g_{1,0},g_{0,1}$, and
$g_{1,1}$ (on $[0,1],[0,1]$, and $[0,1]^{2}$ respectively) which have finite
Hardy-Krause variation (and satisfy the one-sided continuity assumptions
stated in Proposition 2.5). If $f$ is smooth, the functions $g_{1,0},g_{0,1}$,
and $g_{1,1}$ can be matched to the partial derivatives of $f$ as

$g_{1,0}(x_{1})=f^{(1,0)}(x_{1},0)=\frac{\partial f}{\partial
x_{1}}(x_{1},0),~{}~{}~{}~{}~{}g_{0,1}(x_{2})=f^{(0,1)}(0,x_{2})=\frac{\partial
f}{\partial x_{2}}(0,x_{2})$

and

$g_{1,1}(x_{1},x_{2})=f^{(1,1)}(x_{1},x_{2})=\frac{\partial^{2}f}{\partial
x_{1}\partial x_{2}}(x_{1},x_{2}).$

Specializing (31) to $d=2$, we can write the estimator $\hat{f}_{n,V}^{2}$
(somewhat loosely) as

$\displaystyle\hat{f}_{n,V}^{2}\in\mathop{\rm
argmin}_{f}\bigg{\\{}\sum_{i=1}^{n}\big{(}y_{i}-f(x^{(i)})\big{)}^{2}:V_{\mathrm{HK0}}\Big{(}\frac{\partial^{2}f}{\partial
x_{1}\partial x_{2}}(\cdot,\cdot)\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+V\Big{(}\frac{\partial
f}{\partial x_{1}}(\cdot,0)\Big{)}+V\Big{(}\frac{\partial f}{\partial
x_{2}}(0,\cdot)\Big{)}\leq V\bigg{\\}}.$

For smooth functions $f$, the formula (33) for $d=2$ becomes

$\displaystyle
V_{\mathrm{mars}}(f)=\int_{0}^{1}\int_{0}^{1}\big{|}f^{(2,2)}(x_{1},x_{2})\big{|}\,dx_{1}\,dx_{2}$
$\displaystyle+\int_{0}^{1}\big{|}f^{(2,1)}(x_{1},0)\big{|}\,dx_{1}+\int_{0}^{1}\big{|}f^{(2,0)}(x_{1},0)\big{|}\,dx_{1}$
$\displaystyle+\int_{0}^{1}\big{|}f^{(1,2)}(0,x_{2})\big{|}\,dx_{2}+\int_{0}^{1}\big{|}f^{(0,2)}(0,x_{2})\big{|}\,dx_{2}.$

This can be viewed as a second-order version of the Hardy-Krause variation:

$V_{\mathrm{HK0}}(f)=\int_{0}^{1}\int_{0}^{1}\big{|}f^{(1,1)}(x_{1},x_{2})\big{|}\,dx_{1}\,dx_{2}+\int_{0}^{1}\big{|}f^{(1,0)}(x_{1},0)\big{|}\,dx_{1}+\int_{0}^{1}\big{|}f^{(0,1)}(0,x_{2})\big{|}\,dx_{2}.$

###  2.4 Related Work

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 [25] (see also [33, 21, 37]). There are
other ways of generalizing the piecewise linear locally adaptive regression
splines estimator to the multivariate setting as well (see, e.g., Parhi and
Nowak [29, 30] for some recent work).

We have also mentioned that $V_{\mathrm{mars}}(\cdot)$ can be viewed as
second-order Hardy-Krause variation (anchored at $\mathbf{0}$). Our estimation
strategy can thus be seen as second-order Hardy-Krause denoising. In Fang et
al. [14], first-order Hardy-Krause denoising (i.e., least squares estimation
over functions with bounded Hardy-Krause variation) was studied. First-order
Hardy-Krause denoising leads to piecewise constant fits while our method leads
to MARS fits (linear combinations of products of ReLU functions of individual
variables).

The formula (33) implies that our estimator $\hat{f}_{n,V}^{d}$ can be viewed
as a least squares estimator over a smoothness-constrained function class
where smoothness is measured in terms of the $L_{1}$ norms of mixed partial
derivatives of orders up to 2. It is well-known from approximation theory that
$L_{p}$ norm constraints on mixed partial derivatives are very 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.,
Temlyakov [35], Dũng et al. [12], and Bungartz and Griebel [8]). These results
hint that our LASSO variant of the MARS method can be effective in multiple
dimensions, and we state formal rate results in Section 5 demonstrating this.

Lin [24, 23] studied penalized nonparametric regression methods where the
penalty can be seen as a variant of the right-hand side of (33) with the
$L_{1}$ norms replaced by the squared $L_{2}$ norms. Also, the anchoring at
$\mathbf{0}$ in our treatment is handled differently in [24, 23] by means of
averaging. The replacement of the $L_{1}$ norms by the squared $L_{2}$ norms
reduces the class of considered functions significantly, and in particular,
the methods of Lin [24, 23] will not work for fitting functions in
$\mathcal{F}_{\text{mars}}^{d}$. Yang and Tan [44] studied penalized
nonparametric regression for 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, and for functions $f$
that are sufficiently smooth, one of their hierarchical total variations
(corresponding to $m=2$ in their notation) converges to $V_{\mathrm{mars}}(f)$
as the grid resolution becomes arbitrarily small. Ortelli and van de Geer [27]
studied multivariate generalizations of trend filtering in the discrete
setting where the goal is to estimate a multivariate tensor. Their method is
related to the discrete formulation of our estimator described in Proposition
4.5 (see Remark 4.7 for more details). van der Laan et al. [40] discussed
function classes similar to $\mathcal{F}_{\infty-\text{mars}}^{d}$ and norms
similar to $V_{\mathrm{mars}}(\cdot)$ and used them for estimation in settings
that seem different from our classical nonparametric regression framework.

##  3 Existence and Computation

In this section, we prove the existence of our infinite-dimensional LASSO
estimator $\hat{f}_{n,V}^{d}$ (defined in (7)) and provide a method to compute
it via finite-dimensional LASSO algorithms. 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 $x^{(i)},i\in[n]$.
As proved in the next lemma, this observation allows us to restrict attention
to the finite-dimensional subclass of $\mathcal{F}_{\infty-\text{mars}}^{d}$
consisting of the functions (3) where each $\nu_{\alpha}$ is a discrete signed
measure supported on the lattice generated by the design points. For each
$k\in[d]$, let $\mathcal{U}_{k}$ denote the finite subset of $[0,1]$
consisting of all the points $0,x_{k}^{(1)},\dots,x_{k}^{(n)},1$ (recall here
that $x_{k}^{(i)}$ denotes the $k^{th}$ coordinate of the $i^{th}$ design
point $x^{(i)}=(x_{1}^{(i)},\dots,x_{d}^{(i)})$). As there could be ties among
$0,x_{k}^{(1)},\dots,x_{k}^{(n)},1$, we shall write, for some $n_{k}\in[n+1]$,

$\mathcal{U}_{k}=\big{\\{}u^{(k)}_{0},u^{(k)}_{1},\cdots,u^{(k)}_{n_{k}}\big{\\}}\qquad\text{where
$0=u^{(k)}_{0}<\cdots<u^{(k)}_{n_{k}}=1$}.$

Note specially that the cardinality of $\mathcal{U}_{k}$ is $n_{k}+1$, that
$u_{0}^{(k)}$ is always 0, and that $u_{n_{k}}^{(k)}$ is always 1\. The next
lemma (proved in Section 8.2.1) implies that, for the optimization (7), we can
restrict our attention to the functions of the form (3) where each
$\nu_{\alpha}$ is a discrete measure supported on the finite set $\prod_{k\in
S(\alpha)}\mathcal{U}_{k}$.

######  Lemma 3.1.

Suppose we are given a real number $a_{\mathbf{0}}$ and a collection of finite
signed measures $\\{\nu_{\alpha}\\}$ where $\nu_{\alpha}$ is defined on
$[0,1)^{|\alpha|}$ for each $\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$.
Then, there exists a collection of discrete signed measures
$\\{\mu_{\alpha}\\}$ where $\mu_{\alpha}$ is concentrated on $(\prod_{k\in
S(\alpha)}\mathcal{U}_{k})\cap[0,1)^{|\alpha|}$ for each
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$ such that

  1. (i)

$f_{a_{\mathbf{0}},\\{\mu_{\alpha}\\}}(x^{(i)})=f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}(x^{(i)})$
for all $i\in[n]$, and

  2. (ii)

$V_{\mathrm{mars}}(f_{a_{\mathbf{0}},\\{\mu_{\alpha}\\}})\leq
V_{\mathrm{mars}}(f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}})$.

When $\nu_{\alpha}$ is concentrated on $(\prod_{k\in
S(\alpha)}\mathcal{U}_{k})\cap[0,1)^{|\alpha|}$ for each
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$, the function
$f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}$ can be written as

$a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\sum_{l\in\prod\limits_{k\in
S(\alpha)}[0:(n_{k}-1)]}\nu_{\alpha}\big{(}\big{\\{}\big{(}u^{(k)}_{l_{k}},k\in
S(\alpha)\big{)}\big{\\}}\big{)}\cdot\prod_{k\in
S(\alpha)}\big{(}x_{k}-u^{(k)}_{l_{k}}\big{)}_{+},$ (35)

where we use the notation $[p:q]:=\\{p,p+1,\cdots,q\\}$ for two integers
$p\leq q$. Further, the complexity measure
$V_{\mathrm{mars}}(f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}})$ equals

$V_{\mathrm{mars}}(f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}})=\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\sum_{l\in\prod\limits_{k\in
S(\alpha)}[0:(n_{k}-1)]\setminus\\{\mathbf{0}\\}}\Big{|}\nu_{\alpha}\big{(}\big{\\{}\big{(}u^{(k)}_{l_{k}},k\in
S(\alpha)\big{)}\big{\\}}\big{)}\Big{|}.$

The above function (35) is a linear combination of basis functions from the
finite set ${\mathcal{B}}^{d}_{\text{design}:\text{mars}}$ defined in (10),
and the 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, the problem (7) reduces to a finite-dimensional
LASSO problem. Lemma 3.1 then implies that every solution to this finite-
dimensional LASSO problem is also a solution to the original problem (7). A
precise statement is given in the following result, which we prove in Section
8.2.2.

######  Proposition 3.2.

Let

$J=\bigg{\\{}(\alpha,l):\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}\mbox{
and }l\in\prod_{k\in S(\alpha)}[0:(n_{k}-1)]\bigg{\\}}$

and let $M$ be the $n\times|J|$ matrix with columns indexed by $(\alpha,l)\in
J$ such that

$M_{i,(\alpha,l)}=\prod_{k\in
S(\alpha)}\big{(}x^{(i)}_{k}-u^{(k)}_{l_{k}}\big{)}_{+}\qquad\text{for
$i\in[n]$ and $(\alpha,l)\in J$}.$

Let
$(\hat{a}_{\mathbf{0}},\hat{\gamma}_{n,V}^{d})\in{\mathbb{R}}\times{\mathbb{R}}^{|J|}$
be a solution to the following finite-dimensional LASSO problem

$(\hat{a}_{\mathbf{0}},\hat{\gamma}^{d}_{n,V})\in\underset{a_{\mathbf{0}}\in{\mathbb{R}},\gamma\in{\mathbb{R}}^{|J|}}{\mathop{\rm
argmin}}\Bigg{\\{}\left\lVert
y-a_{\mathbf{0}}\mathbf{1}-M\gamma\right\rVert_{2}^{2}:\sum_{\begin{subarray}{c}(\alpha,l)\in
J\\\ l\neq\mathbf{0}\end{subarray}}|\gamma_{\alpha,l}|\leq V\Bigg{\\}},$ (36)

where $y=(y_{i},i\in[n])$ is the vector of observations. Then, the function
$f$ on $[0,1]^{d}$ defined by

$(x_{1},\dots,x_{d})\mapsto
f(x_{1},\cdots,x_{d})=\hat{a}_{\mathbf{0}}+\sum_{(\alpha,l)\in
J}\big{(}\hat{\gamma}^{d}_{n,V}\big{)}_{\alpha,l}\cdot\prod_{k\in
S(\alpha)}\big{(}x_{k}-u^{(k)}_{l_{k}}\big{)}_{+}$ (37)

is a solution to the problem (7). (7) can have multiple solutions, but every
solution $\hat{f}^{d}_{n,V}$ satisfies

$\hat{f}^{d}_{n,V}(x^{(i)})=\hat{a}_{\mathbf{0}}+\big{(}M\hat{\gamma}^{d}_{n,V}\big{)}_{i}=\hat{a}_{\mathbf{0}}+\sum_{(\alpha,l)\in
J}\big{(}\hat{\gamma}^{d}_{n,V}\big{)}_{\alpha,l}\cdot\prod_{k\in
S(\alpha)}\big{(}x^{(i)}_{k}-u^{(k)}_{l_{k}}\big{)}_{+}$

for every $i\in[n]$.

Since it is well-known that a solution to the finite-dimensional LASSO problem
(36) exists and is computable, our estimation problem (37) has a computable
solution as well. Indeed, once we find a solution to the problem (36) via any
finite-dimensional LASSO algorithm, we can construct a solution
$\hat{f}^{d}_{n,V}$ to the problem (7) through the equation (37).

##  4 Alternative Optimization Problem for Equally Spaced Lattice Designs

In this section, we shall describe an alternative finite-dimensional
optimization problem for solving (7) in the case where the design points
$x^{(1)},\dots,x^{(n)}$ form an equally spaced lattice in $[0,1]^{d}$. Equally
spaced lattice designs are useful in applications such as image denoising.
They are also commonly used for theoretical analysis in nonparametric function
estimation (see, e.g., Nemirovski [26]). Our risk results (to be described in
Sections 5 and 6) are stated only for the equally spaced lattice design as
well.

To motivate the general result, let us first consider the univariate case
$d=1$. Let us write the data as $(x^{(i)},y_{i}),i=0,1,\dots,n-1$ with
$x^{(i)}=i/n$ (note that the index $i$ runs from $0$ to $n-1$ as opposed to
$1$ to $n$). In this case, it is well-known that a solution
$\hat{f}^{1}_{n,V}$ to (7) can be obtained by solving the finite-dimensional
optimization problem

$\hat{\theta}=\underset{\theta\in\mathbb{R}^{n}}{\mathop{\rm
argmin}}\bigg{\\{}\sum_{i=0}^{n-1}(y_{i}-\theta_{i})^{2}:n\sum_{i=2}^{n-1}\big{|}(D^{(2)}\theta)_{i}\big{|}\leq
V\bigg{\\}}$ (38)

where

$(D^{(1)}\theta)_{i}:=\theta_{i}-\theta_{i-1},\ i\in[n-1]~{}\text{ and
}~{}(D^{(2)}\theta)_{i}:=(D^{(1)}\theta)_{i}-(D^{(1)}\theta)_{i-1},\
i\in[2:(n-1)]$

and taking

$\hat{f}_{n,V}^{d}(x):=\hat{\theta}_{0}+n\big{(}\hat{\theta}_{1}-\hat{\theta}_{0}\big{)}x+\sum_{l=1}^{n-2}\big{(}\hat{\theta}_{l+1}-2\hat{\theta}_{l}+\hat{\theta}_{l-1}\big{)}(nx-l)_{+}.$

The connection between (7) and (38) can be intuitively seen by noting that

$V_{\mathrm{mars}}(f)=\int_{0}^{1}|f^{\prime\prime}(x)|\,dx\approx
n\sum_{i=2}^{n-1}|(D^{(2)}\theta)_{i}|$

for $\theta_{i}=f(i/n),i\in[0:(n-1)]$ provided $n$ is sufficiently large and
$f$ is smooth. The second-order difference vector $D^{(2)}\theta$ (multiplied
by $n^{2}$) therefore plays the role of the second derivative in the
continuous case.

Let us now consider the case of general $d\geq 1$ where we assume that the
design points $x^{(1)},\dots,x^{(n)}$ are given by the equally spaced lattice:

$\\{x^{(1)},\dots,x^{(n)}\\}=\prod_{k=1}^{d}\Big{\\{}\frac{i_{k}}{n_{k}}:i_{k}\in[0:(n_{k}-1)]\Big{\\}},$
(39)

where $n_{k}\geq 2$ for all $k\in[d]$. It makes sense then to index the data
with the set

$I_{0}:=[0:(n_{1}-1)]\times\dots\times[0:(n_{d}-1)].$ (40)

In this setting, our next result (Proposition 4.5) describes a finite-
dimensional optimization problem for solving (7) that is analogous to (38). As
is clear from (33), the complexity $V_{\mathrm{mars}}(\cdot)$ now involves the
$L_{1}$ norms of the partial derivatives of orders $\beta\in\\{0,1,2\\}^{d}$
with $\max_{j}\beta_{j}=2$. The finite-dimensional optimization problem in
Proposition 4.5 can be seen as an analogous discrete version in which the
continuous partial derivatives are replaced by appropriate discrete
differences. Let us first define these discrete difference operators.

We first define the basic difference operator $D_{d}$ for each positive
integer $d$ as follows. We use the notation $[p:q]=\\{p,p+1,\cdots,q\\}$ for
two integers $p\leq q$.

######  Definition 4.1.

Let $n_{1},\cdots,n_{d}$ be positive integers and $a_{1},\dots,a_{d}$ be
nonnegative integers. For an $n_{1}\dots n_{d}$-dimensional vector $\theta$
that is indexed by
$i\in[a_{1}:(a_{1}+n_{1}-1)]\times\cdots\times[a_{d}:(a_{d}+n_{d}-1)]$, we
define $D_{d}\theta$ to be the $(n_{1}-1)\cdots(n_{d}-1)$-dimensional vector
that is indexed by
$i\in[(a_{1}+1):(a_{1}+n_{1}-1)]\times\cdots\times[(a_{d}+1):(a_{d}+n_{d}-1)]$
and has entries

$(D_{d}\theta)_{i}=\sum_{\delta\in\\{0,1\\}^{d}}(-1)^{\delta_{1}+\cdots+\delta_{d}}\theta_{i-\delta}$

for
$i\in[(a_{1}+1):(a_{1}+n_{1}-1)]\times\cdots\times[(a_{d}+1):(a_{d}+n_{d}-1)]$.

######  Example 4.2 ($d=1$ and $d=2$).

For $\theta\in\mathbb{R}^{n}$ indexed by $i\in[a:(a+n-1)]$, $D_{1}\theta$ is
defined as the $(n-1)$-dimensional vector indexed by $i\in[(a+1),(a+n-1)]$ for
which

$(D_{1}\theta)_{i}=\theta_{i}-\theta_{i-1}$

for $i\in[(a+1):(a+n-1)]$.

For $\theta\in\mathbb{R}^{n_{1}n_{2}}$ indexed by
$(i_{1},i_{2})\in[a_{1}:(a_{1}+n_{1}-1)]\times[a_{2}:(a_{2}+n_{2}-1)]$,
$D_{2}\theta$ is defined as the $(n_{1}-1)(n_{2}-1)$-dimensional vector
indexed by
$(i_{1},i_{2})\in[(a_{1}+1),(a_{1}+n_{1}-1)]\times[(a_{2}+1),(a_{2}+n_{2}-1)]$
for which

$(D_{2}\theta)_{i_{1},i_{2}}=\theta_{i_{1},i_{2}}-\theta_{i_{1}-1,i_{2}}-\theta_{i_{1},i_{2}-1}+\theta_{i_{1}-1,i_{2}-1}$

for
$(i_{1},i_{2})\in[(a_{1}+1):(a_{1}+n_{1}-1)]\times[(a_{2}+1):(a_{2}+n_{2}-1)]$.

Based on the basic difference operators, we define the general difference
operator $D^{(\alpha)}$ for each vector $\alpha$ with nonnegative integer
components.

######  Definition 4.3.

Let $n_{1},\cdots,n_{d}$ be positive integers and suppose that
$\theta\in\mathbb{R}^{n_{1}\cdots n_{d}}$ is indexed by
$i\in[a_{1}:(a_{1}+n_{1}-1)]\times\cdots\times[a_{d}:(a_{d}+n_{d}-1)]$. Then,
we define $D^{(\alpha)}\theta$ in the following recursive way.

  * •

If $\alpha=\mathbf{0}$, then we define $D^{(\mathbf{0})}\theta=\theta$.

  * •

If $\alpha_{k}>0$ for all $k$, then we recursively define
$D^{(\alpha)}\theta=D^{(\alpha-\mathbf{1})}(D_{d}\theta)$.

  * •

Otherwise, let $\theta^{\prime}$ be the lower dimensional slice of $\theta$
consisting of $\theta_{i}$ where $i_{k}=a_{k}$ for all $k\in
S_{0}(\alpha):=\\{j:\alpha_{j}=0\\}$. Index $\theta^{\prime}$ with the set
$\prod_{k\not\in S_{0}(\alpha)}[a_{k}:(a_{k}+n_{k}-1)]$. Next, let
$\alpha^{\prime}$ be the vector with nonnegative integer components obtained
by removing the zero components of $\alpha$. For example, if
$\alpha=(2,0,1,0)$, then $\alpha^{\prime}$ will be $(2,1)$. Then, we can
consider $D^{(\alpha^{\prime})}\theta^{\prime}$, which can be defined
recursively. Let $I^{\prime}$ be the index set of
$D^{(\alpha^{\prime})}\theta^{\prime}$ and construct $I$ from $I^{\prime}$ by
adding the singleton $\\{a_{k}\\}$ to the component where
$[a_{k}:(a_{k}+n_{k}-1)]$ was previously dropped, for every $k\in
S_{0}(\alpha)$. We then define $D^{(\alpha)}\theta$ as follows:

    * –

$D^{(\alpha)}\theta$ is indexed by $i\in I$.

    * –

For $i\in I$,

$(D^{(\alpha)}\theta)_{i}=(D^{(\alpha^{\prime})}\theta^{\prime})_{i^{\prime}},$

where $i^{\prime}$ is obtained from $i$ by dropping out the $k^{th}$ component
for every $k\in S_{0}(\alpha)$. For instance, if $i=(3,0,1,0)$, $i^{\prime}$
will be $(3,1)$.

######  Example 4.4.

Suppose $\alpha=(2,1)$ and $\theta\in\mathbb{R}^{n_{1}n_{2}}$ is indexed by
$(i_{1},i_{2})\in[0:(n_{1}-1)]\times[0:(n_{2}-1)]$. Since $\alpha_{k}>0$ for
all $k$,
$D^{(\alpha)}\theta=D^{(\alpha-\mathbf{1})}(D_{2}\theta)=D^{(1,0)}(D_{2}\theta)$.
Let $\eta=D_{2}\theta$ and note (from Definition 4.1) that $\eta$ is indexed
by $(i_{1},i_{2})\in[1:(n_{1}-1)]\times[1:(n_{2}-1)]$. Next, let
$\eta^{\prime}$ be the lower dimensional slice of $\eta$, i.e.,
$\eta^{\prime}$ is indexed by $i_{1}\in[1:(n_{1}-1)]$ and
$\eta^{\prime}_{i_{1}}=\eta_{i_{1},1}$ for $i_{1}\in[1:(n_{1}-1)]$. Then,
since $D^{(1)}\eta^{\prime}=D^{(0)}(D_{1}\eta^{\prime})=D_{1}\eta^{\prime}$ is
indexed by $i_{1}\in[2:(n_{1}-1)]$, $D^{(1,0)}\eta$ is indexed by
$(i_{1},i_{2})\in[2:(n_{1}-1)]\times\\{1\\}$ and
$(D^{(1,0)}\eta)_{i_{1},1}=(D^{(1)}\eta^{\prime})_{i_{1}}=\eta^{\prime}_{i_{1}}-\eta^{\prime}_{i_{1}-1}=\eta_{i_{1},1}-\eta_{i_{1}-1,1}$
for $i_{1}\in[2:(n_{1}-1)]$. Consequently, $D^{(\alpha)}\theta$ is indexed by
$(i_{1},i_{2})\in[2:(n_{1}-1)]\times\\{1\\}$ and

$\displaystyle(D^{(\alpha)}\theta)_{i_{1},1}$
$\displaystyle=(D^{(1,0)}\eta)_{i_{1},1}=\eta_{i_{1},1}-\eta_{i_{1}-1,1}=(D_{2}\theta)_{i_{1},1}-(D_{2}\theta)_{i_{1}-1,1}$
$\displaystyle=\theta_{i_{1},1}-\theta_{i_{1}-1,1}-\theta_{i_{1},0}+\theta_{i_{1}-1,0}-(\theta_{i_{1}-1,1}-\theta_{i_{1}-2,1}-\theta_{i_{1}-1,0}+\theta_{i_{1}-2,0})$
$\displaystyle=\theta_{i_{1},1}-\theta_{i_{1},0}-2\theta_{i_{1}-1,1}+2\theta_{i_{1}-1,0}+\theta_{i_{1}-2,1}-\theta_{i_{1}-2,0}$

for $i_{1}\in[2:(n_{1}-1)]$.

Let $\theta$ be an $n$-dimensional vector indexed with the set $I_{0}$ (given
by (40)) where $n=n_{1}\cdots n_{d}$. From the definition of the general
difference operators, we can observe that for a $d$-dimensional vector
$\alpha=(\alpha_{1},\cdots,\alpha_{d})$ with nonnegative integer components,
$D^{(\alpha)}\theta$ is indexed with the set
$I_{0}^{(\alpha)}:=I_{1}^{(\alpha)}\times\cdots\times I_{d}^{(\alpha)}$, where

$\displaystyle I_{k}^{(\alpha)}=\begin{cases}\\{\alpha_{k}\\}&\mbox{if
}\alpha_{k}\neq\max_{j}\alpha_{j}\\\ [\max_{j}\alpha_{j}:(n_{k}-1)]&\mbox{if
}\alpha_{k}=\max_{j}\alpha_{j}\end{cases}$

for $k\in[d]$. Thus,

$I_{0}^{(\alpha)}\cap I_{0}^{(\beta)}\neq\emptyset\quad\mbox{if and only
if}\quad\alpha=\beta$

for all $\alpha,\beta$ with $\max_{j}\alpha_{j}=\max_{j}\beta_{j}$, and

$\biguplus_{\begin{subarray}{c}\alpha\in[0:r]^{d}\\\
\max_{j}\alpha_{j}=r\end{subarray}}I_{0}^{(\alpha)}=I_{0}\setminus[0:(r-1)]^{d}.$

Here $\biguplus$ indicates disjoint union. Also, we can check that for each
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$,

$I_{0}^{(\alpha)}=\\{\alpha\\}\uplus\biguplus_{\beta\in
J_{\alpha}}I_{0}^{(\beta)},$

where $J_{\alpha}$ is as defined in (32).

We are now ready to state Proposition 4.5 (proved in Section 8.3).

######  Proposition 4.5.

Let $\hat{\theta}^{d}_{n,V}\in\mathbb{R}^{|I_{0}|}$ be the unique solution to
the problem

$\hat{\theta}^{d}_{n,V}=\underset{\theta\in\mathbb{R}^{|I_{0}|}}{\mathop{\rm
argmin}}\big{\\{}\left\lVert y-\theta\right\rVert_{2}^{2}:\theta\in
C(V)\big{\\}},$ (41)

where

$C(V):=\\{\theta\in\mathbb{R}^{|I_{0}|}:V_{2}(\theta)\leq V\\}$

and

$V_{2}(\theta):=\sum_{\begin{subarray}{c}\beta\in\\{0,1,2\\}^{d}\\\
\max_{j}\beta_{j}=2\end{subarray}}\bigg{[}\bigg{(}\prod_{k=1}^{d}n_{k}^{\beta_{k}-\mathbf{1}\\{\beta_{k}=2\\}}\bigg{)}\cdot\sum_{i\in
I_{0}^{(\beta)}}\big{|}(D^{(\beta)}\theta)_{i}\big{|}\bigg{]}.$ (42)

Then, a solution to the original optimization problem (7) is the function
$\hat{f}$ defined on $[0,1]^{d}$ by

$\hat{f}(x_{1},\cdots,x_{d})=\big{(}\hat{\theta}^{d}_{n,V}\big{)}_{\mathbf{0}}+\sum_{(\alpha,l)\in\tilde{J}}\big{(}H^{(2)}\hat{\theta}^{d}_{n,V}\big{)}_{\tilde{l}+\alpha}\cdot\prod_{k\in
S(\alpha)}(n_{k}x_{k}-l_{k})_{+},$ (43)

where $H^{(2)}\theta$ (for $\theta\in\mathbb{R}^{|I_{0}|}$) is the
$n$-dimensional vector that is indexed by $i\in I_{0}$ and has components

$\displaystyle(H^{(2)}\theta)_{i}=\begin{cases}(D^{(\beta)}\theta)_{i}&\mbox{if
}i\in I_{0}^{(\beta)}\mbox{ for some }\beta\in\\{0,1,2\\}\mbox{ with
}\max_{j}\beta_{j}=2\\\ (D^{(i)}\theta)_{i}&\mbox{if
}i\in\\{0,1\\}^{d},\end{cases}$

$\tilde{J}$ is the set

$\tilde{J}=\bigg{\\{}(\alpha,l):\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}\mbox{
and }l\in\prod_{k\in S(\alpha)}[0:(n_{k}-2)]\bigg{\\}},$

and $\tilde{l}$ is the $d$-dimensional vector with

$\tilde{l}_{k}=\begin{cases}l_{k}&\mbox{if }k\in S(\alpha)\\\
0&\mbox{otherwise}.\end{cases}$

Further, for every solution $\hat{f}^{d}_{n,V}$ to (7), the values of
$\hat{f}^{d}_{n,V}$ at the design points equal the components of
$\hat{\theta}^{d}_{n,V}$, i.e.,

$\hat{f}^{d}_{n,V}\Big{(}\Big{(}\frac{i_{k}}{n_{k}},k\in[d]\Big{)}\Big{)}=\big{(}\hat{\theta}^{d}_{n,V}\big{)}_{i}\qquad\text{for
all $i\in I_{0}$}.$

The formula (42) for $V_{2}(\theta)$ can be seen as a discrete analogue of the
formula (33) for $V_{\mathrm{mars}}(f)$ for smooth $f$. Indeed, for
$\theta_{i}=f(i_{1}/n_{1},\dots i_{d}/n_{d})$ for $i\in I_{0}$, we have

$\bigg{(}\prod_{k=1}^{d}n_{k}^{\beta_{k}-\mathbf{1}\\{\beta_{k}=2\\}}\bigg{)}\cdot\sum_{i\in
I_{0}^{(\beta)}}\big{|}(D^{(\beta)}\theta)_{i}\big{|}\approx\int_{\bar{T}^{(\beta)}}|f^{(\beta)}|$

when $f$ is smooth and $n_{k},k\in[d]$ are large.

######  Example 4.6 (d = 2).

When $d=2$, $V_{2}(\theta)$ can be written as

$\displaystyle V_{2}(\theta)$
$\displaystyle=n_{1}\sum_{i_{1}=2}^{n_{1}-1}\big{|}(D^{(2,0)}\theta)_{i_{1},0}\big{|}+n_{2}\sum_{i_{2}=2}^{n_{2}-1}\big{|}(D^{(0,2)}\theta)_{0,i_{2}}\big{|}$
$\displaystyle\quad+n_{1}n_{2}\sum_{i_{1}=2}^{n_{1}-1}\big{|}(D^{(2,1)}\theta)_{i_{1},1}\big{|}+n_{1}n_{2}\sum_{i_{2}=2}^{n_{2}-1}\big{|}(D^{(1,2)}\theta)_{1,i_{2}}\big{|}+n_{1}n_{2}\sum_{i_{1}=2}^{n_{1}-1}\sum_{i_{2}=2}^{n_{2}-1}\big{|}(D^{(2,2)}\theta)_{i_{1},i_{2}}\big{|}$
$\displaystyle=n_{1}\sum_{i_{1}=2}^{n_{1}-1}\big{|}(\theta_{i_{1},0}-\theta_{i_{1}-1,0})-(\theta_{i_{1}-1,0}-\theta_{i_{1}-2,0})\big{|}$
$\displaystyle\quad+n_{2}\sum_{i_{2}=2}^{n_{2}-1}\big{|}(\theta_{i_{2},0}-\theta_{i_{2}-1,0})-(\theta_{i_{2}-1,0}-\theta_{i_{2}-2,0})\big{|}$
$\displaystyle\quad+n_{1}n_{2}\sum_{i_{1}=2}^{n_{1}-1}\big{|}(\theta_{i_{1},1}-\theta_{i_{1}-1,1}-\theta_{i_{1},0}+\theta_{i_{1}-1,0})-(\theta_{i_{1}-1,1}-\theta_{i_{1}-2,1}-\theta_{i_{1}-1,0}+\theta_{i_{1}-2,0})\big{|}$
$\displaystyle\quad+n_{1}n_{2}\sum_{i_{2}=2}^{n_{2}-1}\big{|}(\theta_{1,i_{2}}-\theta_{0,i_{2}}-\theta_{1,i_{2}-1}+\theta_{0,i_{2}-1})-(\theta_{1,i_{2}-1}-\theta_{0,i_{2}-1}-\theta_{1,i_{2}-2}+\theta_{0,i_{2}-2})\big{|}$
$\displaystyle\quad+n_{1}n_{2}\sum_{i_{1}=2}^{n_{1}-1}\sum_{i_{2}=2}^{n_{2}-1}\big{|}(\theta_{i_{1},i_{2}}-\theta_{i_{1}-1,i_{2}}-\theta_{i_{1},i_{2}-1}+\theta_{i_{1}-1,i_{2}-1})$
$\displaystyle\qquad\qquad\qquad\qquad\quad-(\theta_{i_{1}-1,i_{2}}-\theta_{i_{1}-2,i_{2}}-\theta_{i_{1}-1,i_{2}-1}+\theta_{i_{1}-2,i_{2}-1})$
$\displaystyle\qquad\qquad\qquad\qquad\quad-(\theta_{i_{1},i_{2}-1}-\theta_{i_{1}-1,i_{2}-1}-\theta_{i_{1},i_{2}-2}+\theta_{i_{1}-1,i_{2}-2})$
$\displaystyle\qquad\qquad\qquad\qquad\quad+(\theta_{i_{1}-1,i_{2}-1}-\theta_{i_{1}-2,i_{2}-1}-\theta_{i_{1}-1,i_{2}-2}+\theta_{i_{1}-2,i_{2}-2})\big{|}.$

######  Remark 4.7.

The second-order Vitali trend filtering estimator of Ortelli and van de Geer
[27] is similar to our estimator $\hat{\theta}^{d}_{n,V}$ defined in (41). In
our notation, Ortelli and van de Geer [27] studied the estimator

$\hat{\theta}_{\mathcal{N}_{2}^{\perp}}=\underset{\theta\in\mathbb{R}^{|I_{0}|}}{\mathop{\rm
argmin}}\bigg{\\{}\frac{1}{n}\left\lVert(y-\theta)_{\mathcal{N}_{2}^{\perp}}\right\rVert_{2}^{2}+2\lambda
n\cdot\sum_{i\in
I_{0}^{(2,2,\cdots,2)}}\big{|}(D^{(2,2,\cdots,2)}\theta)_{i}\big{|}\bigg{\\}},$
(44)

where $\lambda>0$ is a tuning parameter and
$(y-\theta)_{\mathcal{N}_{2}^{\perp}}$ denotes the orthogonal projection of
$(y-\theta)$ onto the orthogonal complement $\mathcal{N}_{2}^{\perp}$ of

$\mathcal{N}_{2}:=\big{\\{}\theta\in\mathbb{R}^{|I_{0}|}:(D^{(2,2,\cdots,2)}\theta)_{i}=0\mbox{
for all }i\in I_{0}^{(2,2,\cdots,2)}\big{\\}}.$

The optimization problem (44) is similar to (41) in the sense that it also
aims to find a solution $\theta$ that is close to the observation $y$ and has
small second-order differences at the same time. However, (44) measures the
closeness to the observation via $(y-\theta)_{\mathcal{N}_{2}^{\perp}}$,
instead of $y-\theta$ as in (41). Also, (44) only penalizes
$D^{(2,2,\cdots,2)}\theta$, while (41) constrains $D^{(\beta)}\theta$ for
every $\beta\in\\{0,1,2\\}^{d}$ with $\max_{j}\beta_{j}=2$.

The estimator $\hat{\theta}_{\mathcal{N}_{2}^{\perp}}$ estimates
$\theta^{*}_{\mathcal{N}_{2}^{\perp}}$ but not the whole $\theta^{*}$ (here
$\theta^{*}$ has components $f^{*}(x^{(i)})$, $i\in I_{0}$). In contrast, our
estimator $\hat{\theta}_{n,V}^{d}$ effectively estimates $\theta^{*}$. Ortelli
and van de Geer [27] proved a theoretical result on the accuracy of
$\hat{\theta}_{\mathcal{N}_{2}^{\perp}}$ (as an estimator of
$\theta^{*}_{\mathcal{N}_{2}^{\perp}}$) that can be compared to our risk
results (see Remark 5.2 for details).

##  5 Risk Analysis for Equally Spaced Lattice Designs

We present risk bounds for our infinite-dimensional LASSO estimator
$\hat{f}^{d}_{n,V}$ in this section. We define the risk of an estimator
$\hat{f}_{n}$ based on the data $(x^{(1)},y_{1}),\dots,(x^{(n)},y_{n})$ by

$\mathcal{R}\big{(}\hat{f}_{n},f^{*}\big{)}:=\mathbb{E}\mathcal{L}\big{(}\hat{f}_{n},f^{*}\big{)}\qquad\text{where
$\mathcal{L}\big{(}\hat{f}_{n},f^{*}\big{)}=\frac{1}{n}\sum_{i=1}^{n}\big{(}\hat{f}_{n}(x^{(i)})-f^{*}(x^{(i)})\big{)}^{2}$}.$
(45)

The expectation above is taken over $y_{1},\cdots,y_{n}$ generated according
to (11) with fixed $x^{(1)},\dots,x^{(n)}$. Throughout this section, we assume
that the design points $x^{(1)},\dots,x^{(n)}$ are given by the equally spaced
lattice as in (39).

Our first result states an upper bound of the risk of $\hat{f}^{d}_{n,V}$
under the assumption $f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{d}$ and
$V_{\mathrm{mars}}(f^{*})\leq V$.

######  Theorem 5.1.

Suppose $f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{d}$ and
$V_{\mathrm{mars}}(f^{*})\leq V$ and assume the equally spaced lattice design
(39). The estimator $\hat{f}^{d}_{n,V}$ then satisfies

$\displaystyle\begin{split}\mathcal{R}\big{(}\hat{f}^{d}_{n,V},f^{*}\big{)}&\leq
C_{d}\Big{(}\frac{\sigma^{2}V^{\frac{1}{2}}}{n}\Big{)}^{\frac{4}{5}}\bigg{[}\log\Big{(}2+\frac{Vn^{\frac{1}{2}}}{\sigma}\Big{)}\bigg{]}^{\frac{3(2d-1)}{5}}\\\
&\qquad\
+C_{d}\frac{\sigma^{2}}{n}\bigg{[}\log\Big{(}2+\frac{Vn^{\frac{1}{2}}}{\sigma}\Big{)}\bigg{]}^{\frac{3(2d-1)}{4}}+C_{d}\frac{\sigma^{2}}{n}[\log
n]^{2}\end{split}$ (46)

for some positive constant $C_{d}$ solely depending on $d$.

Note that for fixed $d,\sigma$, and $V$ and sufficiently large $n$, the first
term is the dominant term on the right-hand side of (46), so that

$\mathcal{R}\big{(}\hat{f}^{d}_{n,V},f^{*}\big{)}=O\big{(}n^{-\frac{4}{5}}(\log
n)^{\frac{3(2d-1)}{5}}\big{)},$ (47)

where the multiplicative constant underlying $O(\cdot)$ depends on $d,\sigma$,
and $V$. In the univariate case, we can deduce from Guntuboyina et al. [19,
Theorem 2.1] that

$\mathcal{R}\big{(}\hat{f}^{1}_{n,V},f^{*}\big{)}\leq
C\Big{(}\frac{\sigma^{2}V^{\frac{1}{2}}}{n}\Big{)}^{\frac{4}{5}}+C\frac{\sigma^{2}}{n}\log
n,$

where $C$ is some universal constant. In other words,

$\mathcal{R}\big{(}\hat{f}^{1}_{n,V},f^{*}\big{)}=O\big{(}n^{-\frac{4}{5}}\big{)},$

where the multiplicative constant underlying $O(\cdot)$ depends on $d,\sigma$,
and $V$, as well. Thus, what Theorem 5.1 (proved in Section 8.4.1) tells us is
that for all $d\geq 2$, $\hat{f}^{d}_{n,V}$ can achieve the same rate
$n^{-4/5}$, although it slightly deteriorates by a logarithmic multiplicative
factor depending on $d$. This suggests that our LASSO method for MARS fitting
can avoid the curse of dimensionality and be a useful estimation technique in
higher dimensions.

######  Remark 5.2.

Ortelli and van de Geer [27, Theorem 3.2] proved that the estimator
$\hat{\theta}_{\mathcal{N}_{2}^{\perp}}$ (defined in (44)) has the rate of
convergence

$n^{-\frac{H(d)+3}{2H(d)+3}}(\log n)^{\frac{H(d)}{2H(d)+3}}$ (48)

for estimating $\theta^{*}_{\mathcal{N}_{2}^{\perp}}$ where $\theta^{*}$ has
components $f^{*}(x^{(i)}),i\in I_{0}$. Here $H(d)$ is the harmonic number
defined by $H(d)=1+1/2+\cdots+1/d$. Note that in (48), the exponent of $n$ in
the rate depends on $d$, and it increases as $d$ does. For example, it is
$-4/5$ if $d=1$, $-3/4$ if $d=2$, and $-29/40$ if $d=3$. The rate (48) is
suboptimal for $d\geq 2$ compared to the rate given by (46) whose exponent of
$n$ equals $-4/5$ regardless of the value of $d$.

The key step of our proof of Theorem 5.1 is to find an upper bound of the
metric entropy of $\mathcal{D}_{m}$ (under the $L_{2}$ metric), which is
defined as the collection of all the functions of the form

$(x_{1},\dots,x_{m})\mapsto\int(x_{1}-t_{1})_{+}\dots(x_{m}-t_{m})_{+}\,d\nu(t),$

where $m\in[d]$ is a fixed integer and $\nu$ is a signed measure on
$[0,1]^{m}$ with the total variation $|\nu|([0,1]^{m})\leq 1$. The following
theorem (proved in Section 8.4.2) contains our result on the metric entropy of
$\mathcal{D}_{m}$.

######  Theorem 5.3.

There exist positive constants $C_{m}$ and $\epsilon_{m}$ solely depending on
$m$ such that

$\log N(\epsilon,\mathcal{D}_{m},\|\cdot\|_{2})\leq
C_{m}\epsilon^{-\frac{1}{2}}\Big{[}\log\frac{1}{\epsilon}\Big{]}^{\frac{3(2m-1)}{4}}$

for every $0<\epsilon<\epsilon_{m}$. The log term can be omitted when $m=1$.

######  Remark 5.4.

If the class $\mathcal{D}_{m}$ is altered by replacing $(x-t)_{+}$ with
$\mathbf{1}\\{x\geq t\\}$ and restricting $\nu$ to probability measures, then
one obtains the collection of all the functions of the form

$(x_{1},\dots,x_{m})\mapsto\int\mathbf{1}\\{x_{1}\geq
t_{1}\\}\dots\mathbf{1}\\{x_{m}\geq t_{m}\\}\,d\nu(t)=\nu([\mathbf{0},x]).$

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 et al. [3]. Thus, we are basically extending the argument in
[3] from $\mathbf{1}\\{x\geq t\\}$ to $(x-t)_{+}$.

Theorem 5.3 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 $\mathcal{D}_{m}$ to lower bounds of the small ball
probability of integrated Brownian sheet based on the ideas from Blei et al.
[3, Section 3] and Gao [18, Section 3] and the results from Li and Linde [22,
Theorem 1.2] and Artstein et al. [2, Theorem 5]. The small ball probability of
integrated Brownian sheet here refers to the quantity

$\mathbb{P}\Big{(}\sup_{t\in[0,1]^{m}}|X_{m}(t)|\leq\epsilon\Big{)},$

where $\epsilon>0$ and $X_{m}$ is an $m$-dimensional integrated Brownian sheet
(a description of integrated Brownian sheet is given in Section 8.4.2).
Required bounds on this small ball probability are then obtained using the
results from Dunker et al. [13, Theorem 6] and Chen and Li [10, Theorem 1.2].
Specifically, we show that there exist positive constants $C_{m}$ and
$\epsilon_{m}$ depending only on $m$ such that

$\log\mathbb{P}\Big{(}\sup_{t\in[0,1]^{m}}|X_{m}(t)|\leq\epsilon\Big{)}\geq-
c_{m}\epsilon^{-\frac{2}{3}}\Big{[}\log\frac{1}{\epsilon}\Big{]}^{2m-1}$

for every $0<\epsilon<\epsilon_{m}$. This result, along with the connection
between the metric entropy and the small ball probability, leads to Theorem
5.3, which is the main ingredient in our proof of Theorem 5.1.

Our next result shows that the logarithmic multiplicative factor in the rate
(47) can not be completely removed in the minimax sense. Consider the minimax
risk

$\mathfrak{M}^{d}_{n,V}:=\inf_{\hat{f}_{n}}\sup_{\begin{subarray}{c}f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{d}\\\
V_{\mathrm{mars}}(f^{*})\leq
V\end{subarray}}\mathbb{E}_{f^{*}}\mathcal{L}\big{(}\hat{f}_{n},f^{*}\big{)},$

where the expectation is taken over $y_{1},\cdots,y_{n}$ of (11) and
$\inf_{\hat{f}_{n}}$ denotes the infimum over all estimators $\hat{f}_{n}$ of
$f^{*}$ based on $y_{1},\cdots,y_{n}$. Then, the following theorem (proved in
Section 8.4.3) shows that the supremum risk of every estimator indeed requires
an additional logarithmic multiplicative factor depending on $d$ in addition
to the $n^{-4/5}$ term. We want to note that there is still the gap between
the power $3(2d-1)/5$ of $(\log n)$ in our risk upper bound for
$\hat{f}_{n,V}^{d}$ and the power $4(d-1)/5$ of $(\log n)$ in the minimax
lower bound.

######  Theorem 5.5.

Suppose $\min_{j}n_{j}\geq bn^{1/d}$ for some $0<b\leq 1$ and assume the
equally spaced lattice design (39). Then, there exist positive constants
$C_{b,d}$ depending on $b$ and $d$ and $c_{b,d,\sigma^{2}/V^{2}}$ depending on
$b,d$, and $\sigma^{2}/V^{2}$ such that

$\mathfrak{M}^{d}_{n,V}\geq
C_{b,d}\Big{(}\frac{\sigma^{2}V^{\frac{1}{2}}}{n}\Big{)}^{\frac{4}{5}}\bigg{[}\log\Big{(}\frac{Vn^{\frac{1}{2}}}{\sigma}\Big{)}\bigg{]}^{\frac{4(d-1)}{5}}$

provided $n\geq c_{b,d,\sigma^{2}/V^{2}}$.

In Theorem 5.5, we assume that $\min_{j}n_{j}\geq bn^{1/d}$ for some $0<b\leq
1$. This type of additional assumption is unavoidable. If some $n_{j}$’s are
extremely small compared to the others, then the dimension of the estimation
problem is inherently smaller than $d$ and we can not expect a lower bound as
in Theorem 5.5, which includes a logarithmic multiplicative factor depending
on $d$. Our proof of Theorem 5.5 is based on Assouad’s lemma with a finite set
of functions in
$\\{f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{d}:V_{\mathrm{mars}}(f^{*})\leq
V\\}$ that is constructed by an extension of the ideas in Blei et al. [3,
Section 4].

##  6 Results for the Estimator $\hat{f}_{n,V}^{d,s}$

In this section, we focus on the study of the estimator $\hat{f}_{n,V}^{d,s}$
that was defined in (16). $\hat{f}_{n,V}^{d,s}$ is a constrained least squares
estimator on the class $\mathcal{F}_{\infty-\text{mars}}^{d,s}$ defined in
(15). Intuitively speaking, $\mathcal{F}_{\infty-\text{mars}}^{d,s}$ is the
subcollection of $\mathcal{F}_{\infty-\text{mars}}^{d}$ consisting of all the
functions that have no interaction in variables of order $>s$. This estimator
is especially attractive (compared to $\hat{f}^{d}_{n,V}$) when $s$ is small
(e.g., $s=2$) and $d$ is moderately large.

Our first result for the estimator $\hat{f}_{n,V}^{d,s}$ is Proposition 6.1
below, which guarantees its existence and describes a method for computing it
based on finite-dimensional LASSO. We use here the same notation as in
Proposition 3.2.

######  Proposition 6.1.

Let $J$ and $M$ be as defined in Proposition 3.2 and let
$(\hat{a}_{\mathbf{0}},\hat{\gamma}_{n,V}^{d,s})\in{\mathbb{R}}\times{\mathbb{R}}^{|J|}$
be a solution to the following finite-dimensional LASSO problem

$\displaystyle\begin{split}&(\hat{a}_{\mathbf{0}},\hat{\gamma}^{d,s}_{n,V})\in\underset{a_{\mathbf{0}}\in{\mathbb{R}},\gamma\in{\mathbb{R}}^{|J|}}{\mathop{\rm
argmin}}\Bigg{\\{}\left\lVert
y-a_{\mathbf{0}}\mathbf{1}-M\gamma\right\rVert_{2}^{2}:\sum_{\begin{subarray}{c}(\alpha,l)\in
J\\\ l\neq\mathbf{0}\end{subarray}}|\gamma_{\alpha,l}|\leq V\\\
&\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\mbox{ and
}\gamma_{\alpha,l}=0\mbox{ for all }(\alpha,l)\in J\mbox{ with
}|\alpha|>s\Bigg{\\}}.\end{split}$ (49)

Then, the function $f$ on $[0,1]^{d}$ defined by

$(x_{1},\dots,x_{d})\mapsto
f(x_{1},\cdots,x_{d})=\hat{a}_{\mathbf{0}}+\sum_{\begin{subarray}{c}(\alpha,l)\in
J\\\ |\alpha|\leq
s\end{subarray}}\big{(}\hat{\gamma}^{d,s}_{n,V}\big{)}_{\alpha,l}\cdot\prod_{k\in
S(\alpha)}\big{(}x_{k}-u^{(k)}_{l_{k}}\big{)}_{+}$ (50)

is a solution to the optimization problem (16). Further, (16) can have
multiple solutions, but every solution $\hat{f}^{d,s}_{n,V}$ satisfies

$\hat{f}^{d,s}_{n,V}(x^{(i)})=\hat{a}_{\mathbf{0}}+\big{(}M\hat{\gamma}^{d,s}_{n,V}\big{)}_{i}=\hat{a}_{\mathbf{0}}+\sum_{\begin{subarray}{c}(\alpha,l)\in
J\\\ |\alpha|\leq
s\end{subarray}}\big{(}\hat{\gamma}^{d,s}_{n,V}\big{)}_{\alpha,l}\cdot\prod_{k\in
S(\alpha)}\big{(}x^{(i)}_{k}-u^{(k)}_{l_{k}}\big{)}_{+}$

for every $i\in[n]$.

When the design points $x^{(1)},\cdots,x^{(n)}$ form the equally spaced
lattice (39), $\hat{f}_{n,V}^{d,s}$ can also be computed from the alternative
optimization problem that is given in Proposition 6.2 below. This result is
analogous to and uses the same notation as Proposition 4.5.

######  Proposition 6.2.

Let $\tilde{J}$ and $H^{(2)}\theta$ (for $\theta\in\mathbb{R}^{|I_{0}|}$) be
as defined in Proposition 4.5 and let
$\hat{\theta}^{d,s}_{n,V}\in\mathbb{R}^{|I_{0}|}$ be the unique solution to
the problem

$\hat{\theta}^{d,s}_{n,V}=\underset{\theta\in\mathbb{R}^{|I_{0}|}}{\mathop{\rm
argmin}}\big{\\{}\left\lVert y-\theta\right\rVert_{2}^{2}:\theta\in
C^{s}(V)\big{\\}},$

where

$\displaystyle C^{s}(V)$
$\displaystyle:=\big{\\{}\theta\in\mathbb{R}^{|I_{0}|}:V_{2}(\theta)\leq
V\mbox{ and }(H^{(2)}\theta)_{i}=0\mbox{ for every }i\in I_{0}^{(\alpha)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\mbox{ for all
}\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}\mbox{ with
}|\alpha|>s\big{\\}}.$

Then, the function $f$ on $[0,1]^{d}$ defined by

$(x_{1},\dots,x_{d})\mapsto
f(x_{1},\cdots,x_{d})=\big{(}\hat{\theta}^{d,s}_{n,V}\big{)}_{\mathbf{0}}+\sum_{\begin{subarray}{c}(\alpha,l)\in\tilde{J}\\\
|\alpha|\leq
s\end{subarray}}\big{(}H^{(2)}\hat{\theta}^{d,s}_{n,V}\big{)}_{\tilde{l}+\alpha}\cdot\prod_{k\in
S(\alpha)}\big{(}n_{k}x_{k}-l_{k}\big{)}_{+}$

is a solution to the problem (16). In addition, for every solution
$\hat{f}^{d,s}_{n,V}$ to (16), the values of $\hat{f}^{d,s}_{n,V}$ at the
design points agree with $\hat{\theta}^{d,s}_{n,V}$, i.e,

$\hat{f}^{d,s}_{n,V}\Big{(}\Big{(}\frac{i_{k}}{n_{k}},k\in[d]\Big{)}\Big{)}=\big{(}\hat{\theta}^{d,s}_{n,V}\big{)}_{i}$

for all $i\in I_{0}$.

Next, as in Theorem 5.1, we provide an upper bound of the risk of
$\hat{f}^{d,s}_{n,V}$ under the assumption that
$f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{d,s}$,
$V_{\mathrm{mars}}(f^{*})\leq V$, and the equally spaced lattice design (39).

######  Theorem 6.3.

Suppose $f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{d,s}$ and
$V_{\mathrm{mars}}(f^{*})\leq V$ and assume the equally spaced lattice design
(39). The estimator $\hat{f}^{d,s}_{n,V}$ then satisfies

$\displaystyle\begin{split}\mathcal{R}\big{(}\hat{f}^{d,s}_{n,V},f^{*}\big{)}&\leq
C_{d}\Big{(}\frac{\sigma^{2}V^{\frac{1}{2}}}{n}\Big{)}^{\frac{4}{5}}\bigg{[}\log\Big{(}2+\frac{Vn^{\frac{1}{2}}}{\sigma}\Big{)}\bigg{]}^{\frac{3(2s-1)}{5}}\\\
&\qquad\
+C_{d}\frac{\sigma^{2}}{n}\bigg{[}\log\Big{(}2+\frac{Vn^{\frac{1}{2}}}{\sigma}\Big{)}\bigg{]}^{\frac{3(2s-1)}{4}}+C_{d}\frac{\sigma^{2}}{n}[\log
n]^{2}\end{split}$ (51)

for some positive constant $C_{d}$ solely depending on $d$.

Since the first term is the dominant term on the right-hand side of (51) when
$d,\sigma$, and $V$ are fixed and $n$ is sufficiently large, we have

$\mathcal{R}\big{(}\hat{f}^{d}_{n,V},f^{*}\big{)}=O\big{(}n^{-\frac{4}{5}}(\log
n)^{\frac{3(2s-1)}{5}}\big{)},$ (52)

where the multiplicative constant underlying $O(\cdot)$ depends on $d,\sigma$,
and $V$. Compared to the result (47) for $\hat{f}^{d}_{n,V}$, the order of the
logarithmic factor now decreases to $3(2s-1)/5$. Consequently, when $d$ is
large and $s\ll d$, $\hat{f}^{d,s}_{n,V}$ can be far more effective than
$\hat{f}^{d}_{n,V}$ if $f^{*}$ does not have higher-order interactions in
variables. It may be noted that the constant multiplicative factor underlying
the rate (52) still depends on $d$.

The results of this section are proved by making appropriate modifications to
our proof of the analogous results for $\hat{f}_{n,V}^{d}$. The modifications
needed for the proofs of Propositions 6.1 and 6.2 are straightforward, so we
omit the proofs of these results. Theorem 6.3 is proved in Section 8.7.

##  7 Numerical Experiments

In this section, we provide the results of some numerical experiments
illustrating the performance of $\hat{f}_{n,V}^{d}$ and $\hat{f}_{n,V}^{d,s}$
and compare them with the performance of the usual MARS estimator, which we
will write as $\hat{f}_{\text{mars}}$. The results for simulated data are
presented first, and those for real data follow next. We utilize the R package
$\mathsf{Rmosek}$ (based on interior point convex optimization) to solve the
problem (36) (and (49)) and use the equation (37) (and (50)) to construct our
estimator. For the usual MARS estimator $\hat{f}_{\text{mars}}$, we use the R
package $\mathsf{earth}$ based on Friedman [16] and Friedman [17].

###  7.1 Simulation Studies

For a target function $f^{*}$ defined on $[0,1]^{d}$, we simulate data
$(x^{(1)},y_{1}),\dots,(x^{(n)},y_{n})$ according to the model (11) and
construct each estimator using the whole data twice: once assuming
$V_{\mathrm{mars}}(f^{*})$ is known and once assuming that it is unknown. When
$V_{\mathrm{mars}}(f^{*})$ is unknown, we perform 10-fold cross-validation
with various values of $V$ and choose the value of $V$ that gives the smallest
average mean squared error on the held-out sets. The performance of each
constructed estimator is measured in terms of the loss defined in (45), and to
compute the risk of each estimator, we repeat the above data generation and
estimator construction processes ten times and average the losses computed.

We first start with the case $d=2$ where we present examples of our estimator
$\hat{f}_{n,V}^{2}$ and the usual MARS estimator using heat maps. We annotate
on the plots the loss of the plotted estimators and present the estimated risk
of each estimator in separate tables.

In Figure 1, we consider a specific function
$f^{*}\in\mathcal{F}_{\text{mars}}^{2}$ (which is explicitly given in the
caption to Figure 1) and estimate it from observations on the equally spaced
lattice (39) with $n_{1}=n_{2}=30$. Recall that not only the usual MARS
estimator but also our estimators are functions in
$\mathcal{F}_{\text{mars}}^{d}$ (more precisely, they are finite linear
combinations of functions in ${\mathcal{B}}^{d}_{\text{design}:\text{mars}}$).

In Figure 2, a function
$f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{2}\setminus\mathcal{F}_{\text{mars}}^{2}$
is considered, and we fit it from data on the same lattice as in Figure 1. We
also fit the same function used for Figure 2 with 100 non-equally spaced
design points chosen from the uniform distribution on $[0,1]^{2}$ in Figure 3
(the design points are marked on the heat map for $f^{*}$).

We also conduct experiments for the case when $d>2$, where we use the
estimator $\hat{f}_{n,V}^{d,s}$ with $s=2$ and present involved basis
functions and their coefficients in tables instead of drawing heat maps. In
Table 5, we consider a function
$f^{*}\in\mathcal{F}_{\infty-\text{mars}}^{3,2}\setminus\mathcal{F}_{\text{mars}}^{3}$
and display $\hat{f}_{\text{mars}}$ (whose order of interaction is restricted
to 2) and $\hat{f}_{n,V}^{3,2}$ constructed from observations on the equally
spaced lattice (39) with $n_{1}=n_{2}=n_{3}=12$. Our simulation settings are
described in Table 1.

Table 1:  Functions considered, estimators used, and the type of design points considered in our simulations. | Function class | Our estimator | Design points
---|---|---|---
Figure 1 | $\mathcal{F}_{\text{mars}}^{2}$ | $\hat{f}_{n,V}^{2}$ | Equally spaced lattice
Figure 2 | $\mathcal{F}_{\infty-\text{mars}}^{2}\setminus\mathcal{F}_{\text{mars}}^{2}$ | $\hat{f}_{n,V}^{2}$ | Equally spaced lattice
Figure 3 | $\mathcal{F}_{\infty-\text{mars}}^{2}\setminus\mathcal{F}_{\text{mars}}^{2}$ | $\hat{f}_{n,V}^{2}$ | Uniformly distributed
Table 5 | $\mathcal{F}_{\infty-\text{mars}}^{3,2}\setminus\mathcal{F}_{\text{mars}}^{3}$ | $\hat{f}_{n,V}^{3,2}$ | Equally spaced lattice
Figure 1: The function $f^{*}=10+x_{1}-x_{2}-0.5x_{1}x_{2}+5(x_{1}-0.3)_{+}(x_{2}-0.3)_{+}-4(x_{1}-0.6)_{+}(x_{2}-0.3)_{+}+6(x_{1}-0.6)_{+}(x_{2}-0.6)_{+}$ is shown in the top left plot. We obtain observations from $f^{*}$ (corrupted with Gaussian noise with $\sigma^{2}=0.25$) at each point in the equally spaced lattice design (39) with $n_{1}=n_{2}=30$. Based on the observations, each estimator is constructed. The usual MARS estimator $\hat{f}_{\text{mars}}$ is shown in the top right plot. Our estimator $\hat{f}_{n,V}^{d}$ with $V=V_{\mathrm{mars}}(f^{*})=15$ is shown in the bottom left plot. The bottom right plot shows our estimator $\hat{f}_{n,V}^{d}$ with $V$ selected by cross-validation. Table 2:  This table describes the estimated risk (with standard error) based on ten repetitions of the experiment corresponding to Figure 1. | $\hat{f}_{\text{mars}}$ | $\hat{f}_{n,V}^{2}(V_{\mathrm{mars}}(f^{*}):\text{known})$ | $\hat{f}_{n,V}^{2}(V_{\mathrm{mars}}(f^{*}):\text{unknown})$
---|---|---|---
Estimated risk | $5.1\cdot 10^{-3}\pm 0.6\cdot 10^{-3}$ | $3.3\cdot 10^{-3}\pm 0.1\cdot 10^{-3}$ | $3.3\cdot 10^{-3}\pm 0.5\cdot 10^{-3}$
Figure 2: The function $f^{*}=3+2\sin(\pi x_{1}x_{2})$ is shown in the top left plot. We obtain observations from $f^{*}$ (corrupted with Gaussian noise with $\sigma^{2}=0.25$) at each point in the equally spaced lattice design (39) with $n_{1}=n_{2}=30$. Based on the observations, each estimator is constructed. The usual MARS estimator $\hat{f}_{\text{mars}}$ is shown in the top right plot. Our estimator $\hat{f}_{n,V}^{d}$ with $V=V_{\mathrm{mars}}(f^{*})=32.5$ is shown in the bottom left plot. The bottom right plot shows our estimator $\hat{f}_{n,V}^{d}$ with $V$ selected by cross-validation. Table 3: This table describes the estimated risk (with standard error) based on ten repetitions of the experiment corresponding to Figure 2. | $\hat{f}_{\text{mars}}$ | $\hat{f}_{n,V}^{2}(V_{\mathrm{mars}}(f^{*}):\text{known})$ | $\hat{f}_{n,V}^{2}(V_{\mathrm{mars}}(f^{*}):\text{unknown})$
---|---|---|---
Estimated risk | $7.6\cdot 10^{-3}\pm 0.7\cdot 10^{-3}$ | $4.5\cdot 10^{-3}\pm 0.3\cdot 10^{-3}$ | $4.2\cdot 10^{-3}\pm 0.6\cdot 10^{-3}$
Figure 3: The function $f^{*}=3+2\sin(\pi x_{1}x_{2})$ is shown in the top left plot. We obtain observations from $f^{*}$ (corrupted with Gaussian noise with $\sigma^{2}=0.25$) at each of 100 design points generated according to the uniform distribution on $[0,1]^{2}$ (the design points are displayed in the top left plot). Based on the observations, each estimator is constructed. The usual MARS estimator $\hat{f}_{\text{mars}}$ is shown in the top right plot. Our estimator $\hat{f}_{n,V}^{d}$ with $V=V_{\mathrm{mars}}(f^{*})=32.5$ is shown in the bottom left plot. The bottom right plot shows our estimator $\hat{f}_{n,V}^{d}$ with $V$ selected by cross-validation.  Table 4: This table describes the estimated risk (with standard error) based on ten repetitions of the experiment corresponding to Figure 3. | $\hat{f}_{\text{mars}}$ | $\hat{f}_{n,V}^{2}(V_{\mathrm{mars}}(f^{*}):\text{known})$ | $\hat{f}_{n,V}^{2}(V_{\mathrm{mars}}(f^{*}):\text{unknown})$
---|---|---|---
Estimated risk | $5.7\cdot 10^{-2}\pm 0.8\cdot 10^{-2}$ | $3.2\cdot 10^{-2}\pm 0.4\cdot 10^{-2}$ | $2.8\cdot 10^{-2}\pm 0.4\cdot 10^{-2}$
Table 5:  We consider $d=3$ and the function $f^{*}=3+2\sin(\pi
x_{1}x_{2})+\exp(x_{3})$. We obtain observations from $f^{*}$ (corrupted with
Gaussian noise with $\sigma^{2}=1$) at each point in the equally spaced
lattice design (39) with $n_{1}=n_{2}=n_{3}=12$. Based on the observations,
each estimator is constructed. The involved basis functions and the
coefficients of the MARS estimator $\hat{f}_{\text{mars}}$ is shown in the
left table. The right table shows the the involved basis functions and the
coefficients of our estimator $\hat{f}_{n,V}^{d}$ with $V$ selected to be 14
by cross-validation.

Intercept | 6.608
---|---
$(x_{1}-0.5)_{+}$ | 0.895
$(0.5-x_{1})_{+}$ | -3.744
$(x_{2}-0.4167)_{+}$ | 1.155
$(0.4167-x_{2})_{+}$ | -3.930
$(x_{3}-0.75)_{+}$ | 4.447
$(0.75-x_{3})_{+}$ | -1.516
$(0.5-x_{1})_{+}\cdot(0.4167-x_{2})_{+}$ | 10.24
$(x_{1}-0.6667)_{+}\cdot(x_{2}-0.4167)_{+}$ | -8.405
$(x_{1}-0.5)_{+}\cdot(x_{2}-0.8333)_{+}$ | -19.13

Intercept | 3.888
---|---
$x_{1}$ | 0.161
$x_{2}$ | 0.121
$x_{3}$ | 1.437
$x_{1}\cdot x_{2}$ | 6.039
$x_{1}\cdot x_{3}$ | 0.374
$x_{2}\cdot x_{3}$ | -0.222
$(x_{1}-0.5)_{+}$ | -0.129
$(x_{3}-0.5833)_{+}$ | 0.564
$(x_{3}-0.75)_{+}$ | 1.373
$(x_{1}-0.5833)_{+}\cdot x_{2}$ | -4.036
$(x_{1}-0.6667)_{+}\cdot x_{2}$ | -0.542
$x_{1}\cdot(x_{2}-0.4167)_{+}$ | -0.301
$(x_{1}-0.1667)_{+}\cdot(x_{2}-0.4167)_{+}$ | -2.540
$(x_{1}-0.3333)_{+}\cdot(x_{2}-0.4167)_{+}$ | -3.636
$(x_{1}-0.1667)_{+}\cdot(x_{2}-0.5)_{+}$ | -0.879

Table 6: This table describes the estimated risk (with standard error) based on ten repetitions of the experiment corresponding to Table 5. | $\hat{f}_{\text{mars}}$ | $\hat{f}_{n,V}^{3,2}(V_{\mathrm{mars}}(f^{*}):\text{known})$ | $\hat{f}_{n,V}^{3,2}(V_{\mathrm{mars}}(f^{*}):\text{unknown})$
---|---|---|---
Estimated risk | $1.9\cdot 10^{-2}\pm 0.2\cdot 10^{-2}$ | $1.6\cdot 10^{-2}\pm 0.2\cdot 10^{-2}$ | $1.3\cdot 10^{-2}\pm 0.2\cdot 10^{-2}$

In Figure 1, our estimator with $V$ selected by cross-validation performs the
best. Due to the errors in the observations, the usual MARS method fits the
function that is not smooth around $x_{2}=0.5$. Our estimator for the case
when $V_{\mathrm{mars}}(f^{*})$ is known is smoother than the usual MARS
estimator but still suffers from the same issue. However, the cross-validation
step picks out $V$ smaller than $V_{\mathrm{mars}}(f^{*})$, and it helps
produce the fit that is as smooth as the true function $f^{*}$. In Figure 2,
our estimator with $V=V_{\mathrm{mars}}(f^{*})$ outperforms the other
estimators in terms of the loss. For our estimator for the case when
$V_{\mathrm{mars}}(f^{*})$ is unknown, $V$ bigger than
$V_{\mathrm{mars}}(f^{*})$ is chosen from cross-validation, and such choice of
$V$ brings about slight decline in performance. However, we can see that it
still performs better than the usual MARS estimator. In Figure 3, it is
interesting to observe that our estimator with the cross-validated $V$ is
quite smooth in the regions with limited design points. For example, in the
top left region, while our estimator with the information on
$V_{\mathrm{mars}}(f^{*})$ fluctuates considerably due to the scarcity of
design points, our estimator without the information shows smooth transition.
Also, contrary to the other two estimators, our estimator without the
information on $V_{\mathrm{mars}}(f^{*})$ succeeds in restoring the smooth
shape of the true function $f^{*}$ in the bottom right region. In Table 5, we
can see that as the usual MARS estimator, our estimator with $V$ selected by
cross-validation figures out that there is no interaction between $x_{3}$ and
the other two variables. Note that the only interactions between $x_{3}$ and
the other two variables that appear in our estimator are $x_{1}x_{3}$ and
$x_{2}x_{3}$, which are not penalized in our model. Further, in Table 2, Table
3, Table 4, and Table 6, we can see that the estimated risks of our estimators
are smaller than that of the usual MARS estimator, and the estimated risk of
our estimator for the case when $V_{\mathrm{mars}}(f^{*})$ is unknown is
almost the same as or even smaller than that of our estimator with
$V=V_{\mathrm{mars}}(f^{*})$.

We remark that in these examples, the information on
$V_{\mathrm{mars}}(f^{*})$ does not necessarily lead to better performance for
our estimator. Our estimator for the case when $V_{\mathrm{mars}}(f^{*})$ is
unknown performs at the same level as our estimator with known
$V_{\mathrm{mars}}(f^{*})$. Often, the cross-validation step even figures out
a value of $V$ which yields a better result compared to
$V_{\mathrm{mars}}(f^{*})$. Also, in these examples, we observe that our
estimator outperforms the usual MARS estimator by noticeable margins
regardless of the existence of information on $V_{\mathrm{mars}}(f^{*})$.
Moreover, our estimators in the examples show the almost same level of
sparsity as the corresponding MARS estimators. The usual MARS estimator is
sparser in terms of the usual MARS representation, which allows the functions
$(t_{j}-x_{j})_{+}$. However, if we decompose each basis function
$(t_{j}-x_{j})_{+}$ into $(x_{j}-t_{j})_{+}-(x_{j}-0)_{+}+t_{j}$, the number
of basis functions each MARS estimator uses increases to almost the same
number as our estimators require. For instance, from the conversion, the
number of basis functions of the MARS estimator in Table 5 increases from ten
to thirteen. Overall, the above simulations suggest that our method can
produce estimators that are nearly as interpretable as the usual MARS
estimators while showing comparable and often better accuracy.

The current implementation of our method is however much slower than the usual
MARS method. This is unsurprising considering that the usual MARS method is a
greedy algorithm while our method solves specific optimization problems to
global optimality. We are also using black-box convex optimization techniques
(based on the R package $\mathsf{Rmosek}$) to solve our optimization problems
for computing the estimators. These techniques involve creating the whole
design matrices for LASSO problems, and the number of columns of the design
matrices in our LASSO formulations are usually large. Our current
implementation consumes most of its running time in constructing these design
matrices. For example, the design matrix $M$ involved in our LASSO formulation
for $\hat{f}_{n,V}^{d}$ (described in Proposition 3.2) has the number of
columns equal to

$|J|=\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\prod_{k\in
S(\alpha)}n_{k}=(1+n_{1})\cdots(1+n_{d})-1,$

which can be of order $O(n^{d})$ in the worst case where each $n_{k}=O(n)$
(note however that in the case of a lattice design, the number of columns of
$M$ is $O(n)$ which is much smaller). The estimator $\hat{f}^{d,s}_{n,V}$ can
be computed more efficiently (compared to $\hat{f}_{n,V}^{d}$) if $s\ll d$,
because we can remove in this case the columns of the design matrix $M$ whose
indices do not belong to

$J^{s}:=\bigg{\\{}(\alpha,l):\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\},\
|\alpha|\leq s,\mbox{ and }l\in\prod_{k\in S(\alpha)}[0:(n_{k}-1)]\bigg{\\}}.$

As a result of this removal, the number of columns decreases to

$|J^{s}|=\sum_{\begin{subarray}{c}\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}\\\
|\alpha|\leq s\end{subarray}}\prod_{k\in S(\alpha)}n_{k},$

which is of order $O(n^{s})$ in the worst case (ignoring a multiplicative
factor in $d$). However, our current implementation of $\hat{f}_{n,V}^{d,s}$
is also quite slower compared to the usual MARS estimator. Other approaches
(such as the method used in Boyd et al. [4] for the case of additive models)
can potentially speed up the computations quite a bit, and we shall explore
them in future work.

###  7.2 Ozone Data

Here we use the standard dataset $\mathsf{ozone}$ from the R library
$\mathsf{faraway}$ to illustrate and compare the performance of our estimator
with that of the usual MARS method. We predict atmospheric ozone concentration
level ($\mathsf{ozone}$ or $\mathsf{O_{3}}$, in ppm) from the four covariates:
temperature ($\mathsf{temp}$, in ∘ F), inversion base temperature
($\mathsf{ibt}$, in ∘ F), wind speed ($\mathsf{wind}$, in mph), and visibility
($\mathsf{vis}$, in miles).

As the domain of individual explanatory variables is not $[0,1]$, we first
linearly transform each one into $[0,1]$. In general, there can be multiple
options for the linear transformation. 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 $M-m$. In case the two extreme values of the
domain, $m$ and $M$, are unknown, we simply set $m$ as the minimum and $M$ as
the maximum among observed values. Here for ozone data, we choose the latter
option. For each explanatory variable, we normalize the variable by
subtracting from it the minimal observed value and then dividing it by the
difference between the maximal and the minimal observed values. After fitting
the regression function in the modified domain, we invert the transformation
to write the regression equation in the original domain.

We use our estimator $\hat{f}_{n,V}^{4,2}$ for estimating the regression
function and compare it to the usual MARS estimator $\hat{f}_{\text{mars}}$
whose order of interaction is restricted to 2. We split the data set into a
training set and a test set in advance and use the mean squared error on the
test set for comparison. Out of $n=330$ observations, we use $20\%$
observations as test data and the remaining $80\%$ observations as training
data. We perform 5-fold cross-validation to choose the value of $V$.

The results are presented in Figure 4, Figure 5, and Table 7. In Figure 4, we
plot the two estimated regression equations as functions of the two
explanatory variables $\mathsf{wind}$ (wind speed) and $\mathsf{vis}$
(visibility) while holding $\mathsf{temp}$ (temperature) and $\mathsf{ibt}$
(inverse base temperature) fixed at their median values. We can observe in
Figure 4 that the two estimated regression equations show similar prediction
trend over the domain. Especially, significant changes in prediction trend are
made at the same points $\mathsf{wind}=6$, $\mathsf{vis}=200$, and
$\mathsf{vis}\approx 120$ for both estimators. Also, we find that the two
estimators have similar mean squared errors (the usual MARS estimator:
$2.4\cdot 10^{1}$, our estimator: $2.0\cdot 10^{1}$) on the test set. In
Figure 5, we plot the two estimated regression equations as functions of the
variable $\mathsf{vis}$ (visibility) after we further fix $\mathsf{wind}$
(wind speed) at its median value. In this figure, we can additionally observe
that while the two estimated regression equations show similar prediction
trend, our equation is smoother than the equation from the usual MARS
estimator. In Table 7, we represent our estimator and the usual MARS estimator
in terms of the involved basis functions and their coefficients. Note that
each basis function is transformed back to the original domain from $[0,1]$,
and its coefficient is scaled accordingly. As we can see in Table 7, our
estimator is composed of twenty five basis functions, while the usual MARS
estimator consists of eleven basis functions. Thus, we can say that the usual
MARS is sparser compared to our estimator. However, if we decompose each basis
function $(t_{j}-x_{j})_{+}$ into $(x_{j}-t_{j})_{+}-(x_{j}-0)_{+}+t_{j}$
(before transforming it back to the original domain), the MARS estimator will
then be composed of twenty basis functions. From this perspective, both
estimators show similar levels of sparsity.

We conclude that our estimator and the usual MARS have similar performance on
the $\mathsf{ozone}$ dataset.

Figure 4: Estimated regression equations as functions of $\mathsf{wind}$ and
$\mathsf{vis}$. The usual MARS estimator is on the left and our estimator with
$V$ selected to be 100 by cross-validation is on the right. The variables
$\mathsf{temp}$ and $\mathsf{ibt}$ are fixed at their median values. Figure 5:
Estimated regression equations as functions of $\mathsf{vis}$. The red line
indicates the usual MARS estimator and the blue one indicates our estimator.
The variables $\mathsf{temp}$, $\mathsf{ibt}$, and $\mathsf{wind}$ are fixed
at their median values. Table 7:  The involved basis functions and their
coefficients of the usual MARS estimator (left table) and our estimator (right
table).

Intercept | $1.154\cdot 10^{1}$
---|---
$(\mathsf{temp}-57)_{+}$ | $5.865\cdot 10^{-1}$
$(\mathsf{ibt}-254)_{+}$ | $-6.127\cdot 10^{-1}$
$(\mathsf{ibt}-275)_{+}$ | $5.040\cdot 10^{-1}$
$(\mathsf{wind}-6)_{+}$ | $-6.899\cdot 10^{-1}$
$(\mathsf{vis}-140)_{+}$ | $-1.246\cdot 10^{-1}$
$(\mathsf{vis}-200)_{+}$ | $1.466\cdot 10^{-1}$
$(\mathsf{temp}-69)_{+}(254-\mathsf{ibt})_{+}$ | $-1.143\cdot 10^{-2}$
$(254-\mathsf{ibt})_{+}(100-\mathsf{vis})_{+}$ | $0.724\cdot 10^{-3}$
$(\mathsf{ibt}-250)_{+}(200-\mathsf{vis})_{+}$ | $2.000\cdot 10^{-3}$
$(250-\mathsf{ibt})_{+}(200-\mathsf{vis})_{+}$ | $-0.433\cdot 10^{-3}$

Intercept | $-2.829\cdot 10^{0}$
---|---
$\mathsf{temp}-25$ | $1.029\cdot 10^{-1}$
$(\mathsf{temp}-54)_{+}$ | $0.244\cdot 10^{-1}$
$(\mathsf{temp}-61)_{+}$ | $1.665\cdot 10^{-1}$
$(\mathsf{temp}-63)_{+}$ | $0.051\cdot 10^{-1}$
$(\mathsf{temp}-74)_{+}$ | $-0.550\cdot 10^{-1}$
$\mathsf{ibt}+25$ | $1.828\cdot 10^{-2}$
$(\mathsf{ibt}-180)_{+}$ | $-2.438\cdot 10^{-2}$
$(\mathsf{ibt}-242)_{+}$ | $-2.004\cdot 10^{-2}$
$(\mathsf{ibt}-254)_{+}$ | $-2.236\cdot 10^{-2}$
$\mathsf{wind}$ | $6.374\cdot 10^{-1}$
$(\mathsf{wind}-4)_{+}$ | $-1.745\cdot 10^{-1}$
$(\mathsf{wind}-6)_{+}$ | $-7.299\cdot 10^{-1}$
$\mathsf{vis}$ | $3.163\cdot 10^{-2}$
$(\mathsf{vis}-30)_{+}$ | $-0.077\cdot 10^{-2}$
$(\mathsf{vis}-120)_{+}$ | $-2.764\cdot 10^{-2}$
$(\mathsf{vis}-200)_{+}$ | $2.329\cdot 10^{-2}$
$(\mathsf{temp}-25)\cdot(\mathsf{ibt}+25)$ | $1.007\cdot 10^{-3}$
$(\mathsf{temp}-25)\cdot\mathsf{wind}$ | $-3.207\cdot 10^{-2}$
$(\mathsf{ibt}+25)\cdot\mathsf{wind}$ | $4.005\cdot 10^{-3}$
$(\mathsf{temp}-25)\cdot\mathsf{vis}$ | $2.249\cdot 10^{-4}$
$(\mathsf{ibt}+25)\cdot\mathsf{vis}$ | $-2.730\cdot 10^{-4}$
$(\mathsf{ibt}+25)\cdot(\mathsf{vis}-200)_{+}$ | $1.811\cdot 10^{-4}$
$\mathsf{wind}\cdot\mathsf{vis}$ | $0.716\cdot 10^{-3}$
$\mathsf{wind}\cdot(\mathsf{vis}-140)_{+}$ | $-2.185\cdot 10^{-3}$

###  7.3 Earnings Data

Here using the dataset $\mathsf{ex1029}$ from the R library
$\mathsf{Sleuth3}$, we again compare our estimator and the usual MARS
estimator. We predict weekly earnings ($\mathsf{Weekly}$ $\mathsf{Earnings}$,
in dollars) from years of education ($\mathsf{Education}$, in yrs) and years
of experience ($\mathsf{Experience}$, in yrs). As in the previous subsection,
we first linearly transform each explanatory variable into $[0,1]$ by
subtracting from it the minimal observed value and then dividing it by the
difference between the maximal and the minimal observed values. Once we fit
the regression function in the modified domain, we invert the linear
transformation to represent it in the original domain. Again, out of $n=25437$
observations, we spare $20\%$ observations for test data and use the rest
$80\%$ as training data. Also, as before, we compare the performance of the
estimators based on the mean squared error on the test data. We use
$\hat{f}^{2}_{n,V}$ as our estimator and compare it to the usual MARS
estimator. The value of $V$ is chosen by performing 5-fold cross-validation
with randomly selected $10\%$ of training data.

The results are presented in Figure 6 and Figure 7. In Figure 6, we plot the
two estimated regression equations on the original domain. We can observe in
Figure 6 that the two estimated regression equations show quite similar trend
over the domain. For example, both equations are monotonically increasing with
respect to $\mathsf{Education}$ (with increasing slope) when
$\mathsf{Experience}\leq 40$. Also, we can see that for each value of
$\mathsf{Education}$, both equations reach their peaks around
$\mathsf{Experience}=30$. Indeed, it turns out that the two estimators have
the same mean squared error $1.1\cdot 10^{5}$ on the test set. However, we can
also note in the figure that our estimator gives smoother prediction compared
to the usual MARS estimator, and this is what we can expect in reality. In
Figure 7, we plot the two estimated regression equations as functions of the
explanatory variable $\mathsf{Experience}$ while holding the other variable
$\mathsf{Education}$ fixed at the three different levels
($\mathsf{Education}=6$, $\mathsf{Education}=12$, $\mathsf{Education}=16$). As
we can see in Figure 7, the usual MARS estimator produces counter-intuitive
results. The usual MARS estimator predicts that when
$\mathsf{Experience}\approx 60$, people who are more educated have lower
weekly earnings on average. Contrary to the usual MARS estimator, our
estimator provides more realistic prediction that among people with the same
amount of experience, more educated ones always earn more on average.

We conclude that our estimator produces more realistic results than the usual
MARS estimator for the $\mathsf{Earnings}$ dataset.

Figure 6: Estimated regression equations. The usual MARS estimator is on the
left and our estimator with $V$ selected to be 5000 by cross-validation is on
the right. Figure 7: Estimated regression equations as functions of
$\mathsf{Experience}$. The usual MARS estimator is on the left and our
estimator is on the right. The variable $\mathsf{Education}$ is fixed at the
three different levels ($\mathsf{Education}=6$, $\mathsf{Education}=12$,
$\mathsf{Education}=16$).

##  8 Proofs

###  8.1 Proofs of Lemmas and Propositions in Section 2

####  8.1.1 Uniqueness of Representation (3)

We shall prove the following lemma in this subsection, which shows that the
constant $a_{\mathbf{0}}$ and the signed measures
$\nu_{\alpha},\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$ appearing in
the representation (3) of $f\in\mathcal{F}_{\infty-\text{mars}}^{d}$ are
uniquely determined by $f$.

######  Lemma 8.1.

Suppose $a_{\mathbf{0}}$ and $b_{\mathbf{0}}$ are real numbers and
$\\{\nu_{\alpha}\\}$ and $\\{\mu_{\alpha}\\}$ are collections of finite signed
measures as in (3). If $f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}\equiv
f_{b_{\mathbf{0}},\\{\mu_{\alpha}\\}}$ on $[0,1]^{d}$, then
$a_{\mathbf{0}}=b_{\mathbf{0}}$ and $\nu_{\alpha}=\mu_{\alpha}$ for all
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$.

###### Proof of Lemma 8.1.

By the assumption,

$\displaystyle
a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{[0,1)^{|\alpha|}}\prod_{j\in
S(\alpha)}(x_{j}-t_{j})_{+}\,d\nu_{\alpha}(t^{(\alpha)})$
$\displaystyle\qquad\qquad\qquad=b_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{[0,1)^{|\alpha|}}\prod_{j\in
S(\alpha)}(x_{j}-t_{j})_{+}\,d\mu_{\alpha}(t^{(\alpha)})$

for all $(x_{1},\cdots,x_{d})\in[0,1]^{d}$. Substituting
$x_{1}=\cdots=x_{d}=0$, we get $a_{\mathbf{0}}=b_{\mathbf{0}}$. Suppose for
some positive integer $m$, we have $\nu_{\alpha}=\mu_{\alpha}$ for all
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$ with $|\alpha|\leq m-1$.
Substituting $x_{m+1}=\cdots=x_{d}=0$ and eliminating from the both sides the
components that are already known to be equal, we obtain

$\int_{[0,1)^{m}}\prod_{j=1}^{m}(x_{j}-t_{j})_{+}\,d\nu_{\mathbf{1}_{m}}(t)=\int_{[0,1)^{m}}\prod_{j=1}^{m}(x_{j}-t_{j})_{+}\,d\mu_{\mathbf{1}_{m}}(t).$
(53)

Here $\mathbf{1}_{m}$ indicates the $d$-dimensional vector which has 1 in the
first $m$ components and 0 in the rest. Note that

$\displaystyle\int_{[0,1)^{m}}\prod_{j=1}^{m}(x_{j}-t_{j})_{+}\,d\nu_{\mathbf{1}_{m}}(t)=\int_{\prod\limits_{j=1}^{m}[0,x_{j})}\prod_{j=1}^{m}(x_{j}-t_{j})\,d\nu_{\mathbf{1}_{m}}(t)$
$\displaystyle\qquad\qquad=\int_{\prod\limits_{j=1}^{m}[0,x_{j})}\int_{\prod\limits_{j=1}^{m}(t_{j},x_{j})}1\,ds\,d\nu_{\mathbf{1}_{m}}(t)=\int_{\prod\limits_{j=1}^{m}(0,x_{j})}\int_{\prod\limits_{j=1}^{m}[0,s_{j})}1\,d\nu_{\mathbf{1}_{m}}(t)\,ds.$

Thus, the equation (53) can be written as

$\int_{\prod\limits_{j=1}^{m}(0,x_{j})}\int_{\prod\limits_{j=1}^{m}[0,s_{j})}\,d\nu_{\mathbf{1}_{m}}(t)\,ds=\int_{\prod\limits_{j=1}^{m}(0,x_{j})}\int_{\prod\limits_{j=1}^{m}[0,s_{j})}\,d\mu_{\mathbf{1}_{m}}(t)\,ds.$

By the fundamental theorem of calculus (see, e.g., Rudin [32, Theorem 7.10]),
it follows that

$\int_{\prod\limits_{j=1}^{m}[0,x_{j})}\,d\nu_{\mathbf{1}_{m}}(t)=\int_{\prod\limits_{j=1}^{m}[0,x_{j})}\,d\mu_{\mathbf{1}_{m}}(t)$
(54)

almost everywhere with respect to the Lebesgue measure.

For a subset $R$ of $[0,1)^{m}$, we call it rectangle if it is of the form

$R=\prod_{j=1}^{m}I_{j},$

where for each $j$, $I_{j}=(u_{j},v_{j}),[u_{j},v_{j}],[u_{j},v_{j})$, or
$(u_{j},v_{j}]$ for some real numbers $u_{j}$ and $v_{j}$. Now, let
$\mathcal{P}$ be the collection of all rectangles in $[0,1)^{m}$ and
$\mathcal{Q}$ be the collection of all subsets $E$ of $[0,1)^{m}$ for which

$\nu_{\mathbf{1}_{m}}(E)=\mu_{\mathbf{1}_{m}}(E).$

Using the equation (54), we can show that $\mathcal{P}\subseteq\mathcal{Q}$.
Also, it is clear that $\mathcal{P}$ is a $\pi$-system, and $\mathcal{Q}$ is a
$\lambda$-system (or Dynkin system) on $[0,1)^{m}$. Therefore, by Dynkin’s
$\pi$-$\lambda$ theorem, we have $\sigma(\mathcal{P})\subseteq\mathcal{Q}$.
Since $\sigma(\mathcal{P})$ is indeed the collection of all Borel sets in
$[0,1)^{m}$, we can conclude that $\nu_{\mathbf{1}_{m}}=\mu_{\mathbf{1}_{m}}$.
Consequently, induction argument completes the proof. ∎

####  8.1.2 Proof of Proposition 2.3

Our proof of Proposition 2.3 is based on Aistleitner and Dick [1, Theorem 3],
which we restate as Theorem 8.2 below. Aistleitner and Dick [1, Theorem 3]
connects functions on $[0,1]^{m}$ that has finite $\mathrm{HK0}$ variation and
is coordinate-wise right-continuous to cumulative distribution functions of
finite signed measures on $[0,1]^{m}$. We shall use this theorem again for
proving Lemma 2.4.

######  Theorem 8.2 (Aistleitner and Dick [1], Theorem 3).

Suppose a real-valued function $g$ defined on $[0,1]^{m}$ has finite
$\mathrm{HK0}$ variation and is coordinate-wise right-continuous. Then, there
exists a unique finite signed measure $\nu$ on $[0,1]^{m}$ such that

$g(x)=\nu([\mathbf{0},x])\qquad\text{for every $x\in[0,1]^{m}$}.$ (55)

Also, if a real-valued function $g$ defined on $[0,1]^{m}$ is of the form (55)
for some finite signed measure $\nu$ on $[0,1]^{m}$, then $g$ has finite
$\mathrm{HK0}$ variation and

$|\nu|([0,1]^{m})=V_{\mathrm{HK0}}(g)+|g(\mathbf{0})|.$

###### Proof of Proposition 2.3.

Assume that a real-valued function $g$ defined on $[0,1]^{m}$ has finite
$\mathrm{HK0}$ variation, is coordinate-wise right-continuous, and is
coordinate-wise left-continuous at each point
$x\in[0,1]^{m}\setminus[0,1)^{m}$ with respect to all the $j^{\text{th}}$
coordinates where $x_{j}=1$. Since $g$ has finite $\mathrm{HK0}$ variation and
is coordinate-wise right-continuous, by Theorem 8.2, there exists a finite
signed measure $\bar{\nu}$ on $[0,1]^{m}$ such that

$g(x)=\bar{\nu}([\mathbf{0},x])$

for all $x\in[0,1]^{m}$. Let $\mu$ be the finite signed measure on $[0,1]^{m}$
for which

$\mu(E)=\bar{\nu}(E\cap[0,1)^{m}).$

Then, it follows that

$g(x)=\mu([\mathbf{0},x])$

for all $x\in[0,1]^{m}$. Indeed, if $x\in[0,1)^{m}$, then

$g(x)=\bar{\nu}([\mathbf{0},x])=\bar{\nu}([\mathbf{0},x]\cap[0,1)^{m})=\mu([\mathbf{0},x]).$

Otherwise, without loss of generality, let
$x=(1,\cdots,1,x_{k+1},\cdots,x_{m})$ where $x_{k+1},\cdots,x_{m}$ $<1$ for
some $k\in[m]$. Then, by the left-continuity of $g$,

$\displaystyle g(x)$ $\displaystyle=\lim_{y_{1}\uparrow
1}\cdots\lim_{y_{k}\uparrow 1}g(y_{1},\cdots,y_{k},x_{k+1},\cdots,x_{m})$
$\displaystyle=\lim_{y_{1}\uparrow 1}\cdots\lim_{y_{k}\uparrow
1}\bar{\nu}([\mathbf{0},(y_{1},\cdots,y_{k},x_{k+1},\cdots,x_{m})])$
$\displaystyle=\bar{\nu}([0,1)^{k}\times[0,x_{k+1}]\times\cdots\times[0,x_{m}])$
$\displaystyle=\bar{\nu}([\mathbf{0},x]\cap[0,1)^{m})$
$\displaystyle=\mu([\mathbf{0},x]),$

as well. We thus have the two signed measure $\bar{\nu}$ and $\mu$ satisfying
(55), and the uniqueness in Theorem 8.2 implies that $\bar{\nu}=\mu$, i.e.,

$\bar{\nu}(E)=\mu(E)=\bar{\nu}(E\cap[0,1)^{m}).$

Therefore, if we let $\nu$ be the finite signed measure on $[0,1)^{m}$ for
which

$\nu(E)=\bar{\nu}(E),$

then for every $x\in[0,1]^{m}$,

$g(x)=\bar{\nu}([\mathbf{0},x])=\bar{\nu}([\mathbf{0},x]\cap[0,1)^{m})=\nu([\mathbf{0},x]\cap[0,1)^{m})$

as desired. The uniqueness of such signed measure directly follows from the
uniqueness of the signed measure in Theorem 8.2.

Next, suppose a real-valued function $g$ on $[0,1]^{m}$ is of the form (27)
for some signed measure $\nu$ on $[0,1)^{m}$. Let $\bar{\nu}$ be the finite
signed measure on $[0,1]^{m}$ for which

$\bar{\nu}(E)=\nu(E\cap[0,1)^{m}).$

Then, clearly,

$g(x)=\nu([\mathbf{0},x]\cap[0,1)^{m})=\bar{\nu}([\mathbf{0},x])$

for every $x\in[0,1]^{m}$. By Theorem 8.2, $g$ thus has finite $\mathrm{HK0}$
variation and

$|\bar{\nu}|([0,1]^{m})=V_{\mathrm{HK0}}(g)+|g(\mathbf{0})|.$

Also, since $|g(\mathbf{0})|=|\nu(\\{0\\})|=|\nu|(\\{0\\})$ and
$|\bar{\nu}|([0,1]^{m})=|\nu|([0,1)^{m})$, it follows that

$V_{\mathrm{HK0}}(g)=|\nu|([0,1)^{m}\setminus\\{\mathbf{0}\\}).$

∎

####  8.1.3 Proof of Lemma 2.4

We use the following lemma to prove Lemma 2.4. This lemma also plays a crucial
role in proving Lemma 2.6. We provide our proof of Lemma 8.3 right after the
proof of Lemma 2.4.

######  Lemma 8.3.

If a real-valued function $g$ defined on $[0,1]^{m}$ is smooth in the sense
that $g^{(\alpha)}$ exists, is irrespective of the order of differentiation,
and is continuous on $[0,1]^{m}$ for every $\alpha\in\\{0,1\\}^{m}$, then $g$
satisfies

$g(x_{1},\cdots,x_{m})=g(\mathbf{0})+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{k\in
S(\alpha)}[0,x_{k}]}g^{(\alpha)}(\tilde{t}^{(\alpha)})\,dt^{(\alpha)}$

for all $x=(x_{1},\cdots,x_{m})\in[0,1]^{m}$. Here $t^{(\alpha)}=(t_{k},k\in
S(\alpha))$ and $\tilde{t}^{(\alpha)}$ is its $m$-dimensional extension
defined by

$\displaystyle\tilde{t}^{(\alpha)}_{k}=\begin{cases}t_{k}&\mbox{if }k\in
S(\alpha)\\\ 0&\mbox{otherwise}.\end{cases}$

###### Proof of Lemma 2.4.

First, note that since $g$ is smooth in the sense of Lemma 8.3, we have

$\begin{split}g(x_{1},\cdots,x_{m})&=g(\mathbf{0})+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{k\in
S(\alpha)}[0,x_{k}]}g^{(\alpha)}(\tilde{t}^{(\alpha)})\,dt^{(\alpha)}\end{split}$
(56)

for all $x\in[0,1]^{m}$. For each
$\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}$, let $\mu_{\alpha}$ be the
signed measure on $(0,1]^{|\alpha|}$ defined by

$d\mu_{\alpha}=g^{(\alpha)}(\tilde{t}^{(\alpha)})\,dt^{(\alpha)}.$

Also, for each $\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}$, consider the
projection $\pi_{\alpha}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{|\alpha|}$ for
which

$\pi_{\alpha}\big{(}(t_{k},k\in[m])\big{)}=(t_{k},k\in S(\alpha))$

and let $R^{(\alpha)}=\prod_{k=1}^{m}R_{k}$ where

$\displaystyle R_{k}=\begin{cases}(0,1]&\mbox{if
}\alpha_{k}=\max_{j}\alpha_{j}\\\ \\{0\\}&\mbox{otherwise}.\end{cases}$

Next, we patch together these signed measures and form the signed measure
$\mu$ on $[0,1]^{m}$ by

$\mu(E)=g(\mathbf{0})\cdot 1\\{\mathbf{0}\in
E\\}+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\mu_{\alpha}(\pi_{\alpha}(E\cap
R^{(\alpha)})).$

Then, we can note that

$\displaystyle\begin{split}|\mu|([0,1]^{m})&=|g(\mathbf{0})|+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}|\mu_{\alpha}|((0,1]^{|\alpha|})\\\
&=|g(\mathbf{0})|+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{(0,1]^{|\alpha|}}|g^{(\alpha)}(\tilde{t}^{(\alpha)})|\,dt^{(\alpha)}\\\
&=|g(\mathbf{0})|+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\bar{T}^{(\alpha)}}|g^{(\alpha)}|.\end{split}$
(57)

By (56), we also have

$\displaystyle g(x_{1},\cdots,x_{m})$
$\displaystyle=g(\mathbf{0})+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{k\in
S(\alpha)}(0,x_{k}]}g^{(\alpha)}(\tilde{t}^{(\alpha)})\,dt^{(\alpha)}$
$\displaystyle=g(\mathbf{0})+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{k\in
S(\alpha)}(0,x_{k}]}\,d\mu_{\alpha}$
$\displaystyle=g(\mathbf{0})+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\mu_{\alpha}\Big{(}\prod\limits_{k\in
S(\alpha)}(0,x_{k}]\Big{)}$ $\displaystyle=g(\mathbf{0})\cdot
1\\{\mathbf{0}\in[\mathbf{0},x]\\}+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\mu_{\alpha}\Big{(}\pi_{\alpha}\big{(}[\mathbf{0},x]\cap
R^{(\alpha)}\big{)}\Big{)}$ $\displaystyle=\mu([\mathbf{0},x])$

for all $x\in[0,1]^{m}$. Therefore, Theorem 8.2 implies that $g$ has finite
$\mathrm{HK0}$ variation and

$|\mu|([0,1]^{m})=V_{\mathrm{HK0}}(g)+|g(\mathbf{0})|.$

Combining this with the equation (57), we can derive that

$V_{\mathrm{HK0}}(g)=\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\bar{T}^{(\alpha)}}|g^{(\alpha)}|.$

∎

###### Proof of Lemma 8.3.

Note that for $\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}$,

$\int_{\prod\limits_{k\in
S(\alpha)}[0,x_{k}]}g^{(\alpha)}(\tilde{t}^{(\alpha)})\,dt^{(\alpha)}=\sum_{\delta\in\prod\limits_{k\in
S(\alpha)}\\{0,1\\}}(-1)^{\sum\limits_{k\in S(\alpha)}\delta_{k}}\cdot
g\big{(}\tilde{x}_{\delta}^{(\alpha)}\big{)},$

where $\tilde{x}_{\delta}^{(\alpha)}$ is the $m$-dimensional vector with

$\big{(}\tilde{x}_{\delta}^{(\alpha)}\big{)}_{k}=\begin{cases}(1-\delta_{k})x_{k}&\mbox{if
}k\in S(\alpha)\\\ 0&\mbox{otherwise}\end{cases}$

for $k\in[m]$. It can be easily proved by the fundamental theorem of calculus,
and we use the smoothness assumption on $g$ here. For example, if $m=2$ and
$\alpha=(1,1)$,

$\displaystyle\int_{[0,x_{1}]\times[0,x_{2}]}g^{(1,1)}$
$\displaystyle(t_{1},t_{2})\,d(t_{1},t_{2})=\int_{0}^{x_{1}}\int_{0}^{x_{2}}g^{(1,1)}(t_{1},t_{2})\,dt_{2}\,dt_{1}$
$\displaystyle=\int_{0}^{x_{1}}\big{(}g^{(1,0)}(t_{1},x_{2})-g^{(1,0)}(t_{1},0)\big{)}\,dt_{1}$
$\displaystyle=\int_{0}^{x_{1}}g^{(1,0)}(t_{1},x_{2})\,dt_{1}-\int_{0}^{x_{1}}g^{(1,0)}(t_{1},0)\,dt_{1}$
$\displaystyle=g(x_{1},x_{2})-g(0,x_{2})-g(x_{1},0)+g(0,0)$
$\displaystyle=\sum_{\delta\in\\{0,1\\}^{2}}(-1)^{\delta_{1}+\delta_{2}}\cdot
g\big{(}\tilde{x}_{\delta}^{(1,1)}\big{)},$

and we use that $g,g^{(1,0)},g^{(0,1)}$, and $g^{(1,1)}$ are continuous on
$[0,1]^{2}$ for this argument. Thus,

$\displaystyle g(\mathbf{0})+$
$\displaystyle\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{k\in
S(\alpha)}[0,x_{k}]}g^{(\alpha)}(\tilde{t}^{(\alpha)})\,dt^{(\alpha)}$
$\displaystyle=g(\mathbf{0})+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\sum_{\delta\in\prod\limits_{k\in
S(\alpha)}\\{0,1\\}}(-1)^{\sum\limits_{k\in S(\alpha)}\delta_{k}}\cdot
g\big{(}\tilde{x}_{\delta}^{(\alpha)}\big{)}$
$\displaystyle=\sum_{\alpha\in\\{0,1\\}^{m}}\sum_{\delta\in\prod\limits_{k\in
S(\alpha)}\\{0,1\\}}(-1)^{\sum\limits_{k\in S(\alpha)}\delta_{k}}\cdot
g\big{(}\tilde{x}_{\delta}^{(\alpha)}\big{)}$
$\displaystyle=\sum_{\alpha^{\prime}\in\\{0,1\\}^{m}}g(\tilde{x}^{(\alpha^{\prime})})\cdot\sum_{\begin{subarray}{c}\alpha\in\\{0,1\\}^{m}\\\
\alpha\geq\alpha^{\prime}\end{subarray}}(-1)^{|\alpha|-|\alpha^{\prime}|}.$

Here the last equality results from counting the number of
$g(\tilde{x}^{(\alpha^{\prime})})$ for each $\alpha^{\prime}\in\\{0,1\\}^{m}$.
Note that if $\alpha^{\prime}\neq\mathbf{1}$,

$\sum_{\begin{subarray}{c}\alpha\in\\{0,1\\}^{m}\\\
\alpha\geq\alpha^{\prime}\end{subarray}}(-1)^{|\alpha|-|\alpha^{\prime}|}=\sum_{\delta\in\prod\limits_{k\in
S_{0}(\alpha^{\prime})}\\{0,1\\}}(-1)^{\sum\limits_{k\in
S_{0}(\alpha^{\prime})}\delta_{k}}=(1-1)^{|S_{0}(\alpha^{\prime})|}=0,$

where $S_{0}(\alpha^{\prime})=\\{j\in[m]:\alpha^{\prime}_{j}=0\\}$, and that
if $\alpha^{\prime}=\mathbf{1}$,

$\sum_{\begin{subarray}{c}\alpha\in\\{0,1\\}^{m}\\\
\alpha\geq\alpha^{\prime}\end{subarray}}(-1)^{|\alpha|-|\alpha^{\prime}|}=1.$

As a result, we can obtain that

$g(\mathbf{0})+\sum_{\alpha\in\\{0,1\\}^{m}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{k\in
S(\alpha)}[0,x_{k}]}g^{(\alpha)}(\tilde{t}^{(\alpha)})\,dt^{(\alpha)}=g(\tilde{x}^{(\mathbf{1})})=g(x_{1},\cdots,x_{m})$

as desired. ∎

####  8.1.4 Proof of Proposition 2.5

###### Proof of Proposition 2.5.

Suppose $f\in\mathcal{F}_{\infty-\text{mars}}^{d}$. Then, there exist
$a_{\mathbf{0}}\in\mathbb{R}$ and a collection of finite signed measures
$\\{\nu_{\alpha}:\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}\\}$ as in (3)
such that $f=f_{a_{\mathbf{0}},\\{\nu_{\alpha}\\}}$. For each
$\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}$, define the function
$g_{\alpha}$ on $[0,1]^{|\alpha|}$ by

$g_{\alpha}(x^{(\alpha)})=\nu_{\alpha}([\mathbf{0},x^{(\alpha)}]\cap[0,1)^{|\alpha|})$

for $x^{(\alpha)}\in[0,1]^{|\alpha|}$. Clearly, $g_{\alpha}$ is coordinate-
wise right-continuous and coordinate-wise left-continuous at each point
$x^{(\alpha)}\in[0,1]^{|\alpha|}\setminus[0,1)^{|\alpha|}$ with respect to all
the $j^{\text{th}}$ coordinates where $x_{j}=1$. Also, by Proposition 2.3, we
can see that $g_{\alpha}$ has finite $\mathrm{HK0}$ variation. Moreover, we
have

$\displaystyle f(x_{1},\cdots,x_{d})$
$\displaystyle=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{[0,1)^{|\alpha|}}\prod_{j\in
S(\alpha)}(x_{j}-t_{j})_{+}\,d\nu_{\alpha}(t^{(\alpha)})$
$\displaystyle=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{j\in
S(\alpha)}[0,x_{j})}\prod_{j\in
S(\alpha)}(x_{j}-t_{j})\,d\nu_{\alpha}(t^{(\alpha)})$
$\displaystyle=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{j\in
S(\alpha)}[0,x_{j})}\int_{\prod\limits_{j\in
S(\alpha)}[t_{j},x_{j})}1\,ds^{(\alpha)}\,d\nu_{\alpha}(t^{(\alpha)})$
$\displaystyle=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{j\in
S(\alpha)}[0,x_{j})}\int_{\prod\limits_{j\in
S(\alpha)}[0,s_{j}]}1\,d\nu_{\alpha}(t^{(\alpha)})\,ds^{(\alpha)}$
$\displaystyle=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{j\in
S(\alpha)}[0,x_{j})}\nu_{\alpha}([\mathbf{0},s^{(\alpha)}])\,ds^{(\alpha)}$
$\displaystyle=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{j\in
S(\alpha)}[0,x_{j}]}g_{\alpha}(s^{(\alpha)})\,ds^{(\alpha)}$

for all $x\in[0,1]^{d}$. Thus, $f$ can be represented in the form (29) with
$a_{0}\in\mathbb{R}$ and the collection of functions
$\\{g_{\alpha}:\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}\\}$ that has
the desired properties.

Conversely, suppose a real-valued function $f$ on $[0,1]^{d}$ is of the form

$f(x_{1},\cdots,x_{d})=a_{\mathbf{0}}+\sum_{\alpha\in\\{0,1\\}^{d}\setminus\\{\mathbf{0}\\}}\int_{\prod\limits_{j\in
S(\alpha)}[0,x_{j}]}g_{\alpha}(t^{(\alpha)})\,dt^{(\alpha)}$

for some $a_{\mathbf{0}}\in\mathbb{R}$ and for some collection of functions