# Predicting Census Survey Response Rates via Interpretable Nonparametric
Additive Models with Structured Interactions

Shibal Ibrahim MIT Department of Electrical Engineering and Computer Science.
email<EMAIL_ADDRESS>Rahul Mazumder MIT Sloan School of Management,
Operations Research Center and MIT Center for Statistics. email:
<EMAIL_ADDRESS>Peter Radchenko University of Sydney. email:
<EMAIL_ADDRESS>Emanuel Ben-David United States Census Bureau.
email<EMAIL_ADDRESS>
(August, 2021)

###### Abstract

Accurate and interpretable prediction of survey response rates is critically
important from an operational standpoint. The US Census Bureau’s well-known
ROAM application uses a linear regression model trained on the US Census
Planning Database data to identify hard-to-survey areas. A 2012 crowdsourcing
competition revealed that machine learning methods based on ensembles of
regression trees led to the best performance in predicting survey response
rates; however, the corresponding models could not be adopted for the intended
application due to limited interpretability. In this paper, we present new
interpretable statistical methods to predict, with high accuracy, response
rates in surveys. We propose sparse nonparametric additive models with
pairwise interactions via $\ell_{0}$-regularization, as well as hierarchically
structured variants that provide enhanced interpretability. Our framework
provides a promising approach for the Census Bureau application by improving
upon the linear model approach and the black-box machine learning methods from
the nationwide competition. Experiments based on the US Census Planning
Database demonstrate that our methods lead to high-quality predictive models
that permit actionable interpretability for different segments of the
population. Our models automatically identify useful interactions among key
factors that have been reported in earlier work by the Census Bureau. These
interactions lead to improved predictions and offer useful insights into
existing studies by the Census Bureau on understanding different segments of
the population. In addition to being useful from an interpretability
standpoint, our models lead to predictions that are comparable to popular
black-box machine learning methods based on gradient boosting and neural
networks. Despite strong methodological underpinnings, our models are
computationally challenging—our proposed novel scalable algorithms for
learning these computational tasks may be of independent interest. Our code
implementation in python is available at
https://github.com/ShibalIbrahim/Additive-Models-with-Structured-Interactions.

appendixplaceholder6 appendixplaceholder-1

##  1 Introduction

Context and background. Sample surveys and censuses are primary data sources
in social science studies. However, low and often unpredictable response rates
in surveys remain a continual source of concern (Erdman and Bates, 2016;
Smith, 2014; Tourangeau et al., 2014) – see Figure 1 for an illustration on
the American Community Survey. Tourangeau et al. (2014) discuss a multitude of
factors that make parts of the population hard-to-survey – such factors are
often used to improve sampling strategies, questionnaire designs, recruitment
methods, and the language in which the interview is conducted, among others.
Erdman and Bates (2016) emphasize the usefulness of having an indicator for
hard-to-survey areas to guide targeted surveying (including oversampling),
staff recruitment strategies, and targeted follow-ups for non-responders. For
major campaigns such as the decennial US Census, this approach can help guide
resource allocation for advertisements and building community partnership
activities. Eliciting responses from non-self-responding households through
follow-up operations can be very costly and time consuming. The Census Bureau
estimates that a single percentage increase in the self-response rate amounts
to roughly $85$ million dollars saved in personal follow-up costs (US Census
Bureau, 2009; Bates and Mulry, 2011). Apart from the cost, the quality of
proxy enumerations and imputations is typically significantly lower than that
of self-responses (Mule, 2012).

Figure 1: The 2013-2017 American Community Survey self-response rates for all
tracts in the continental United States. The North in general, and the Upper
Midwest and the Northeast in particular, have higher self-response rates than
the rest of the country. Tracts with lower self-response rates are visible in
many states – in particular, in the South and in the Mountain region.

.

Both the UK Office for National Statistics and the US Census Bureau have
created measures that help quantify the difficulty in gathering data across
different geographic areas (Erdman and Bates, 2016). A report by Bruce et al.
(2001) from the US Census Bureau introduced a hard-to-count (HTC) score for
identifying difficult to enumerate segments of the population. The HTC score,
which is based on 12 carefully chosen covariates (for example, housing
variables and socio-demographic/economic indicators) has found its use in the
planning for the 2010 Census and many other national surveys. The HTC score,
however, has some limitations and was later improved to the low-response-
score, or LRS (Erdman and Bates, 2016). The LRS plays a key role in the US
Census Bureau’s Response Outreach Area Mapper (ROAM)
application111https://www.census.gov/library/visualizations/2017/geo/roam.html
to identify hard-to-survey areas. The LRS appears to have been partly
motivated by the 2012 nationwide competition, organized by the Census Bureau
in partnership with Kaggle, on predicting the mail-return rates (Erdman and
Bates, 2016). This competition aimed to solicit models that were highly
predictive of the decennial census response rates, easily replicable,
inherently interpretable, and consistent for use in the field at various
levels of geography (Erdman and Bates, 2016). Even though the winning models
from the competition had good predictive performance, they lacked
reproducibility and interpretability. Some of the winning models included
covariates that are not publicly available; in particular, they were not
chosen from within the US Census Bureau Planning Database. Additionally, as
pointed out by Erdman and Bates (2016), the winning models included covariates
that were good predictors but not actionable222These include covariates such
as the nearest neighboring block group return rates and margins of error for
various estimates. Models that incorporate such features are not actionable..
Interpretability suffered due to a multitude of factors. First, the winning
models were based on complex ensemble methods (for example, random forests and
gradient boosting). Second, these methods employed a large number of
covariates (the best model had nearly 340), which hurt the model parsimony.
After careful analysis, Erdman and Bates (2016) proposed a linear model based
on 25 covariates—LRS is the prediction from this model. While this model
suffered in terms of predictive accuracy compared to the black-box machine
learning methods, it was highly interpretable and led to useful actionable
insights. For example, Erdman and Bates (2016) highlights three different
block-groups (Columbia Heights, Trinidad, Anacostia) in the District of
Columbia that have similar LRS (predictions) but their characteristics vary,
indicating a need for different actions to increase self-response. Columbia
Heights has a large Hispanic population, suggesting they could benefit from
forms and advertising in Spanish. Anacostia has a stable population (low
relocation to the region) and hence could benefit from community partnerships.
On the other hand, Trinidad is characterized as a region in transition with
high mobility of young population. This area could potentially benefit from
web-based advertisements for internet-based response from the tech-savvy
younger population.

Our aim is to predict self-response rates across all Census tracts in the US
using operational, socio-economic, and demographic characteristics from the US
Census Planning Database. Our goal is to propose well-grounded interpretable
statistical models that deliver strong predictive performance and, thus, can
be used to complement the currently used LRS metric. We seek to create
statistical models with strong methodological underpinnings that lend
operational insights into factors that influence self-response rates in
surveys. This can hopefully facilitate the downstream goal of having a census
that is cost effective and results in improved coverage.

Statistical models. Our statistical estimators are based on Additive Models
(Hastie and Tibshirani, 1987; Hastie et al., 2001), or AMs, with smooth
nonlinear components that include nonlinear pairwise interactions between
covariates. Given response $y\in\mathbb{R}$ and feature-vector
$\mathbf{x}:=(x_{1},\ldots,x_{p})\in\mathbb{R}^{p}$, we model the conditional
mean function as

$\mathbb{E}(y|\mathbf{x})=\sum_{j\in[p]}f_{j}(x_{j})+\sum_{j<k}f_{j,k}(x_{j},x_{k}),$
(1)

where $f_{j}$ and $f_{j,k}$ are unknown smooth univariate and bivariate
functions, respectively. Drawing inspiration from linear model settings (Bach
et al., 2012), we propose new methodology to estimate the unknown functional
components via an optimization problem with structural constraints arising
from interpretability considerations. To obtain a sparse model with few
predictors, that is, with many of the components $\\{f_{j}\\}$ and
$\\{f_{j,k}\\}$ estimated as exactly zero, we present a novel
$\ell_{0}$-regularized approach333For recent examples of such approaches in
the context of linear regression, see Bertsimas et al. (2016); Mazumder and
Radchenko (2017); Mazumder et al. (2017); Hazimeh and Mazumder (2020a) and the
references therein., penalizing the number of nonzero components in the model.
In addition, we explore a refined notion of interpretability in the context of
interaction modeling – namely, hierarchical sparsity – inspired by its usage
in the context of linear models (Hazimeh and Mazumder, 2020b; Bien et al.,
2013). To our knowledge, this is the first methodological exploration of
combinatorially structured statistical models in the context of nonparametric
additive modeling with interactions.

Applied contributions. Our empirical analysis on the US Census Bureau Planning
Database demonstrates the strong appeal of our proposed approach as a
potential tool for predicting self-response rates in the context of the US
Census Bureau application studied herein.

Our approach leads to improved prediction of the self-response rates when
compared to the linear regression methods that appear to be used by the US
Census Bureau (Erdman and Bates, 2016). The improvements are due both to the
departure from linearity and the presence of nonlinear interactions.
Importantly, the prediction accuracy of our interpretable models is at par
with that of black-box machine learning methods such as neural networks and
gradient boosted decision trees, which topped the nationwide Kaggle
competition organized by the Census Bureau (Erdman and Bates, 2016).

Our proposed framework leads to simple interpretable models and allows for
automated variable selection. Our models are substantially more compact than
the natural competitors—we use $8$-$20$ times fewer interaction effects than
the corresponding linear models, and about half as many covariates as off-the-
shelf machine learning methods (e.g, gradient boosting, or feed forward neural
networks) or the sparse _linear_ models with interaction effects. Our models
also use a fewer number of covariates (by about 30%) compared to sparse
nonparametric additive models (without interactions). We use our models to
gather useful operational insights into the factors that influence response
rates in surveys across different segments of the population. In particular,
our models automatically identify interactions between many of the key factors
that have been used in prior publications from the Census Bureau (Bates and
Mulry, 2011) to derive meaningful clustering of the population. These
interaction effects appear to result in improved predictions of survey
response rates and help complement the existing studies by the Census Bureau
(Bates and Mulry, 2011; Kulzick et al., 2019) on understanding different
segments of the population. In summary, our work appears to address many of
the challenges of the US Census application by delivering interpretable
parsimonious models with good predictive performance.

Methodological contributions. We propose a new family of estimators based on
nonparametric additive models with interactions under combinatorial
constraints that promote parsimony and interpretability. We establish
statistical guarantees for the resulting estimators and present large-scale
algorithms to obtain high-quality solutions for our proposed learning tasks,
significantly expanding the current computational landscape for nonparametric
additive models with pairwise interactions. Our approach addresses key
computational challenges posed by the large-scale setting of the Census
dataset, with $\approx 10^{5}$ nonparametric interaction components and
$\approx 10^{5}$ observations. The code implementing our algorithms is
provided in the supplementary material.

Related Work. There is an impressive body of methodological and theoretical
work on using convex $\ell_{1}$-based approaches to fit sparse nonparametric
AMs without interactions (see, for example, Meier et al., 2009; Ravikumar et
al., 2009; Huang et al., 2010; Zhao and Liu, 2012; Yuan and Zhou, 2016, and
the references therein). Even with main-effects alone, these convex
optimization-based methods face computational challenges for the problem-
scales we seek to address444Based on our experience, currently available
software (for example, R package SAM) encounters numerical difficulties for
instances where $n$ is on the order of thousands, and $p$ is on the order of
hundreds.. This prevents a potential user from realizing the full potential of
nonparametric AMs in large-scale settings arising in various practical
applications in the modern-day world. In terms of statistical properties, the
$\ell_{0}$-based estimators can offer improvements over their
$\ell_{1}$-counterparts in certain regimes. Recent work by Hazimeh et al.
(2021) studies the $\ell_{0}$-regularized framework for variable selection in
nonparametric AMs without interactions. In the presence of interactions, which
is the setting we focus on, the problem of variable selection in nonparametric
AMs becomes considerably more challenging. Incorporating hierarchically
structured interactions leads to additional computational challenges. To the
best of our knowledge, our algorithmic framework for sparse nonparametric AMs
with interactions is novel and is likely to be of independent interest from a
computational methodology perspective in applications beyond the one studied
herein.

Organization. Section 2 presents the different statistical models pursued in
this paper. Section 2.3 investigates the statistical properties of our
proposed estimators. Section 2.4 discusses how to obtain solutions to the
corresponding large-scale discrete optimization problems. Section 3 discusses
numerical results for the US Census Bureau application that motivates this
study. Additional technical details are provided in the supplementary
material.

##  2 Statistical Models and Methodology

Here we discuss the statistical models pursued in this paper. Section 2.1
gives an overview of AMs with nonlinear main effects and pairwise
interactions, along with optimization formulations associated with the
estimation procedures. Section 2.2 presents our new models to incorporate
sparsity in the main and interaction effects. Section 2.2.2 discusses
obtaining sparse interactions under strong hierarchy.

###  2.1 Smooth additive models with pairwise interactions

Given data $\\{(y_{i},\mathbf{x}_{i})\\}_{1}^{n}$, our key objective is to
learn a multivariate conditional mean function
$f(\mathbf{x}):=\mathbb{E}(y|\mathbf{x})$, where
$f:\mathbb{R}^{p}\mapsto\mathbb{R}$ is an unknown smooth function (Wahba,
1990). It is well known that such functions become difficult to estimate even
for moderately high $p$ due to the curse of dimensionality—therefore we will
focus on a smaller class of models corresponding to additive structures
(Hastie and Tibshirani, 1987; Stone, 1986). A popular approach considered in
the literature estimates a nonparametric additive model containing main-
effects only, with $f(\mathbf{x})=\sum_{j=1}^{p}f_{j}(x_{j})$, where each
$f_{j}$ is an unknown univariate smooth function of the $j$-th coordinate in
$\mathbf{x}$, namely $x_{j}$. In various applications, however, nonparametric
AMs based on main-effects alone may not lead to a sufficiently rich
representation for predicting the outcome of interest: including interaction
terms can lead to better predictive models that still remain interpretable to
a practitioner (Hastie et al., 2001; Hastie and Tibshirani, 1987). AMs with
pairwise interactions are an useful tool in applied statistical modeling with
various applications in medical sciences and healthcare, e-commerce
applications, recommender system problems, and sentiment analysis, among
others (Bien et al., 2013; Hastie et al., 2015).

The need for interactions. Exploratory analysis on US Census Planning Database
Tract dataset suggests that interaction effects across features can lead to
improved predictions of the survey self-response rate. For example, there is
interaction between the percentages of “People who do not speak English well”
and “Renters” in the area. In areas with a relatively high concentration of
poor English speakers (e.g, $\geq 5.4\%$), the self-response rate decreases on
an average by $0.33\%$ for a unit increase in the percentage of renters. On
the other hand, when the concentration of poor English speakers is relatively
low (e.g, $\leq 1.6\%$), the self-response rate decreases at the rate of
$0.21\%$. Similarly, there is a strong interaction effect between covariates
“Single-unit households” and “Household moved in 2010 or later” in terms of
predicting the low self-response rate. The potential importance of
incorporating interaction effects when predicting the self-response rate has
been highlighted in an earlier paper from the Census Bureau (Erdman and Bates,
2016), although the paper does not pursue statistical modeling involving
interactions.

An additive model with nonlinear main effects and pairwise interactions
extends the traditional AM framework with main effects alone (Hastie et al.,
2001; Radchenko and James, 2010), and is given by model (1), where the unknown
components $\\{f_{j}\\}_{j}$ and $\\{f_{j,k}\\}_{j,k}$ need to be estimated
from the data. This leads to two key challenges. The presence of
$O(p^{2})$-many unknown nonparametric functions results in statistical
challenges even for a moderate value of $p$. Additional regularization (for
example, in the form of sparsity in the components) may be necessary to obtain
a reliable statistical model with good generalization properties (cf Section
2.2). Furthermore, as we mention in Section 1, estimating model (1) leads to
severe computational challenges for large problems (for example, those with
$n\approx 10^{5}$ and $p^{2}\approx 10^{5}$, similar to the instances we
consider in our applications) – we discuss how we address these challenges in
Section 2.4.

We assume that the components $\\{f_{j}\\}$ and $\\{f_{j,k}\\}$ are smooth
(for example, twice continuously differentiable). For illustration, consider
Figure 2, which shows that the marginal fits for the self-response rate appear
to be well-approximated by smooth nonlinear functions. Using squared
$\ell_{2}$-loss as the data fidelity term, the task of learning model (1) can
be expressed as the following infinite dimensional optimization problem:

$\displaystyle\min_{\begin{subarray}{c}f_{j}\in{\mathcal{C}}_{1},\forall j\\\
f_{j,k}\in\mathcal{C}_{2},\forall
j<k\end{subarray}}~{}\frac{1}{n}\sum_{i\in[n]}\Big{(}y_{i}-\sum_{j\in[p]}f_{j}(x_{ij})-\sum_{j<k}f_{j,k}(x_{ij},x_{ik})\Big{)}^{2}+\lambda\Big{[}\sum_{j\in[p]}\Omega(f_{j})+\sum_{j<k}\Omega(f_{j,k})\Big{]},$
(2)

where the main-effects $f_{j}\in{\mathcal{C}}_{1}$ for all $j$ and the
interaction effects $f_{j,k}\in\mathcal{C}_{2}$ for all $(j,k)$; with
${\mathcal{C}}_{1}$, ${\mathcal{C}}_{2}$ denoting the sets of smooth functions
in 1D and 2D respectively. Above, $\Omega$ is a roughness penalty measuring
the smoothness level555For example, $\Omega(f_{j})=\int
f_{j}^{\prime\prime}(x_{j})^{2}dx_{j}$ and
$\Omega(f_{ij})=\int_{x_{j}x_{k}}(\partial^{2}f_{j,k}/\partial
x_{j}^{2})^{2}+(\partial^{2}f_{j,k}/\partial x_{j}\partial
x_{k})^{2}+(\partial^{2}f_{j,k}/\partial x_{k}^{2})^{2}dx_{j}dx_{k}$. of
components $f_{j}$ and $f_{j,k}$ and $\lambda\geq 0$ is a regularization
parameter controlling the smoothness of the functions.

|  |
---|---|---
Figure 2: _Panels [Left]-[Right] illustrate marginal nonparametric fits for
the self-response rate output variable versus three covariates. Each marginal
fit, displayed on a scatter plot with a solid blue line, clearly suggests a
nonlinear relationship of the output vs the individual covariate (we note that
the covariates are standardized.) The $x$-axis corresponds to: [Left] “Persons
of Hispanic origin in the ACS”; [Middle] “Number of households that have only
a smartphone and no other computing device”; [Right] “Persons 25 years and
over with college degree or higher in the ACS”._

A finite dimensional quadratic program. Problem (2) can be written as a
finite-dimensional quadratic program. Using splines to model each of the main-
effects and the interaction-effects, one can write the main-effects and
interaction-effects as a linear combination of suitable bases elements.
Additional details can be found in the supplementary material. Let
$\mathbf{f}_{j}=(f_{j}(x_{1j}),\ldots,f_{j}(x_{nj}))$ and
$\mathbf{f}_{j,k}=(f_{j,k}(x_{1j},x_{1k}),\ldots,f_{j,k}(x_{nj},x_{nk}))$
denote the evaluations of the main effect component $f_{j}$ and the
interaction component $f_{j,k}$, respectively, at the $n$ data points. Then,
$\mathbf{f}_{j}=\mathbf{B}_{j}\boldsymbol{\beta}_{j}$, where
$\mathbf{B}_{j}\in\mathbb{R}^{n\times K_{j}}$ is the model matrix and
$\boldsymbol{\beta}_{j}\in\mathbb{R}^{K_{j}}$ is the vector of coefficients
for each main effect component. Similarly,
$\mathbf{f}_{j,k}=\mathbf{B}_{j,k}\boldsymbol{\theta}_{j,k}$, where
$\mathbf{B}_{j,k}\in\mathbb{R}^{n\times K_{j,k}}$ is the model matrix and
$\boldsymbol{\theta}_{j,k}\in\mathbb{R}^{K_{j,k}}$ is the vector of
coefficients for each interaction effect component. Here, $K_{j}$ and
$K_{j,k}$ denote the corresponding dimensions of the bases – in our
implementation, all the values $K_{j}$ are taken to be the same and all the
$K_{j,k}$ are taken to be the same as well. We use a penalty function to
control the smoothness (i.e., integral of the squared second derivative) of
the 1D and 2D components—see supplementary material for details. Writing
$\boldsymbol{\beta}$ for the vector obtained by stacking together the
coefficients $\boldsymbol{\beta}_{j}$, $j\in[p]$ for the main-effects, and
defining the vector $\boldsymbol{\theta}$ for the interaction effects
analogously, we express the optimization problem (2) as follows:

$\displaystyle\small\min_{\boldsymbol{\beta},\boldsymbol{\theta}}g_{\lambda_{1}}(\boldsymbol{\beta},\boldsymbol{\theta})\;\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\;\frac{1}{n}\Big{\|}\boldsymbol{y}-\big{[}\sum_{j\in[p]}\boldsymbol{B}_{j}\boldsymbol{\beta}_{i}+\sum_{j<k}\boldsymbol{B}_{j,k}\boldsymbol{\theta}_{j,k}\big{]}\Big{\|}_{2}^{2}+$
$\displaystyle\lambda_{1}\big{[}\sum_{j\in[p]}\boldsymbol{\beta}_{j}^{T}\boldsymbol{S}_{j}\boldsymbol{\beta}_{j}+\sum_{j<k}\boldsymbol{\theta}_{j,k}^{T}\boldsymbol{S}_{j,k}\boldsymbol{\theta}_{j,k}\big{]}.$
(3)

Here, $\boldsymbol{S}_{j}=\boldsymbol{D}_{j}^{T}\boldsymbol{D}_{j}$ and
$\boldsymbol{S}_{j,k}=(\boldsymbol{D}_{j}^{T}\boldsymbol{D}_{j})\otimes\boldsymbol{I}_{k}+\boldsymbol{I}_{j}\otimes(\boldsymbol{D}_{k}^{T}\boldsymbol{D}_{k})$
are the smoothness penalty matrices for the main effects and the interaction
components, respectively. In (3), for convenience, we use the same smoothness
penalty $\lambda_{1}$ for both the main and the interaction effects, though in
general they may be taken to be different.

###  2.2 Parsimonious models via $\ell_{0}$-penalization

Here we study $\ell_{0}$-type estimators that limit the number of components
in the additive models introduced earlier. As discussed earlier, we seek to
limit the number of main and interaction effects in our model both for
statistical and interpretability reasons.

####  2.2.1 Sparse pairwise interactions

While a significant amount of impressive work has been devoted towards
studying sparsity in the context of linear models, sparsity in nonlinear
models has received less attention. Interestingly, we observe that the notion
of parsimony is linked to the model being used and changes, for example,
depending on whether we use a linear interaction model of the form
$\mathbb{E}(y|\mathbf{x})=\sum_{j}x_{j}\beta_{j}+\sum_{j<k}x_{j}x_{k}\beta_{jk}$
or model (1) that has nonlinear components. For motivation, we consider Figure
3, which presents our findings on a 2019 US Census Bureau Planning Database
dataset with $40$ covariates666These covariates were selected based on
discussions with researchers at the US Census Bureau.. These features include
covariates used for the low-response-score (Erdman and Bates, 2014), important
covariates highlighted in Appendix C of the 2019 US Census Bureau Planning
Database Documentation (US Census Bureau, 2019), plus some additional
covariates capturing internet penetration and urbanization. Specifically, we
observe that:

  * •

Our proposed AM approach with _nonlinear_ main-effects and interactions777This
corresponds to an estimate available from (4). We present the model leading to
the best prediction performance on the validation set – see Section 3 for
details. results in a substantially more compact model than a _linear_ model
with interactions.

  * •

The above nonparametric additive model leads to better test predictions
compared to its linear model counterpart.

  * •

The best (based on validation tuning with prediction error) model with linear
main and interaction effects contains $37$ main and $555$ interaction effects.
On the other hand, using nonparametric AMs, we need a total of 36 nonlinear
main and interaction effects to obtain a similar test prediction performance.
The nonlinear model also uses fewer covariates compared to the linear model
counterpart.

The above observations appear to suggest that the nonlinear model with sparse
interactions has advantages over its linear model counterpart for predicting
survey self-response rates. As discussed in Section 3 the above findings
generally carry over to the expanded dataset with a larger number of features.

|
---|---

Figure 3: _Sparsity pattern of the main and interaction effects presented in a
$p\times p$ matrix: a black square on the diagonal indicates the presence of a
main effect, and an off-diagonal black square indicates the presence of an
interaction effect in the joint model. [Left] Panel illustrates the sparsity
pattern of a Lasso model with main and interaction effects. There are 37 main
and 555 interaction effects in the optimal model. [Right] Panel illustrates
the sparsity pattern of a nonlinear AM with main and interaction effects,
i.e., model (4). There are 8 main and 92 interaction effects in the optimal
model. Model (4) has prediction performance similar to the Lasso model, with
only 3 main and 33 interaction effects. Nonlinear models lead to much more
compact and hence, easier to interpret models compared to linear models. Both
models were trained on a 2019 US Census Bureau Planning Database dataset
(predicting the tract-level self-response rate) with $p=40$ covariates and
$74,000$ observations._

To obtain an additive model with a small number of main and interaction
effects, we consider the following estimator:

$\displaystyle\min_{\boldsymbol{\beta},\boldsymbol{\theta}}G(\boldsymbol{\beta},\boldsymbol{\theta})\;\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\;g_{\lambda_{1}}(\boldsymbol{\beta},\boldsymbol{\theta})+\lambda_{2}\Big{[}\sum_{j\in[p]}\mathbb{1}\bigl{[}\boldsymbol{\beta}_{j}\neq\boldsymbol{0}\bigr{]}+\alpha\sum_{j<k}\mathbb{1}\bigl{[}\boldsymbol{\theta}_{j,k}\neq\boldsymbol{0}\bigr{]}\Big{]},$
(4)

where $\mathbb{1}[\cdot]$ is an indicator function, $\lambda_{2}$ controls the
number of components selected, and
$\alpha\in\mathopen{[}1,\infty\mathclose{)}$ controls the tradeoff between the
number of main effects and the number of interaction effects in the model. We
note that Problem (4) can be formulated as a mixed integer program (MIP)
(Wolsey and Nemhauser, 1999), and hence can be solved to optimality with
modern commercial solvers (for example, Gurobi, Mosek) for small/moderate
scale problems. Recently, tailored nonlinear branch-and-bound techniques
(Hazimeh et al., 2020) have been shown to be promising for solving
$\ell_{0}$-sparse linear regression problems. Since we intend to compute
solutions to (4) for a family of tuning parameters $(\lambda_{1},\lambda_{2})$
at scale, we consider high-quality approximate solutions (cf Section 2.4).
Once a good solution to (4) is available, we can employ MIP techniques to
improve the solution and/or certify the quality of the solution, extending the
local search techniques presented in Hazimeh and Mazumder (2020a) for
$\ell_{0}$-sparse linear regression problems. To our knowledge, this is the
first work to pursue a methodological investigation of estimator (4).

Related Work. Due to the combinatorial nature of (4), several convex
relaxations888In terms of existing implementations, R package SAM presents
specialized algorithms for a convex relaxation of (4) without interactions.
SAM would not run on the dataset we consider here. of this problem have been
studied in the special case when no interactions are present (Ravikumar et
al., 2009; Zhao and Liu, 2012; Meier et al., 2009; Buhlmann and van de Geer,
2011). Hazimeh et al. (2021) considers $\ell_{0}$-formulations for the setting
without interaction effects and demonstrates the merits of using
$\ell_{0}$-based formulations over convex relaxation-based approaches. Despite
the appeal of estimator (4), computational challenges appear to be a key
limiting factor in exploring the usefulness of this model in applications
arising in practice. In Section 2.4 we present a new algorithm for problem
(4).

####  2.2.2 Incorporating strong hierarchy constraints

Problem (4) limits the total number of main and interaction effects and works
well in our experiments (see Section 3 for the details) in terms of obtaining
a sparse model with good predictive performance. In terms of variable
selection properties, however, (4) can lead to the inclusion of an interaction
effect, say, $\\{(j,k)\\}$ where at least one of the corresponding main-
effects $\\{j\\}$ or $\\{k\\}$ is excluded from the model. This may be
somewhat problematic from an interpretation viewpoint. Hence, it may be
desirable to enforce additional constraints in (4), imposing the so-called
_hierarchy constraints_ (McCullagh and Nelder, 1989; Bien et al., 2013). Here,
we consider the well-known _strong hierarchy_ constraint, where an interaction
effect $\\{(j,k)\\}$ is included in the model only if the corresponding main
effects are also included. In addition to improved interpretation, strong
hierarchy can reduce the _effective_ number of features in the model, which in
turn can reduce the operational costs associated with data collection (Bien et
al., 2013; Hazimeh and Mazumder, 2020b).

For example, based on our analysis, covariates “Males” and “Females” appear to
be unnecessary for self-response prediction when the strong hierarchy
constraint is imposed via model (5), while a related covariate “Households
with female without spouse” remains important. In contrast, the LRS uses both
“Males” and “Households with female without spouse” (Erdman and Bates,
2016)—both features appear to be important perhaps due to the use of a linear
regression model. Similarly, covariates pertaining to civilian employment
status across different segments/age-groups of the population, are excluded
from our nonlinear AM when we learn the (sparse) main and interaction effects
under strong hierarchy.

To enforce strong hierarchy into model (4), we consider the following
estimator:

$\displaystyle\min_{\boldsymbol{\beta},\boldsymbol{\theta}}~{}~{}$
$\displaystyle~{}~{}\,g_{\lambda_{1}}(\boldsymbol{\beta},\boldsymbol{\theta})+\lambda_{2}\big{(}\sum_{j\in[p]}\mathbb{1}[\boldsymbol{\beta}_{j}\neq\boldsymbol{0}]+\alpha\sum_{j<k}\mathbb{1}[\boldsymbol{\theta}_{j,k}\neq\boldsymbol{0}]\big{)}$
(5a) s.t.
$\displaystyle~{}~{}\,\boldsymbol{\theta}_{j,k}\neq\boldsymbol{0}~{}~{}\Longrightarrow~{}~{}\boldsymbol{\beta}_{j}\neq\boldsymbol{0}~{}~{}\&~{}~{}\boldsymbol{\beta}_{k}\neq\boldsymbol{0}~{}~{}~{}~{}\forall
j<k,~{}j\in[p],~{}k\in[p].$ (5b)

We note that (5) differs from (4) in the additional strong hierarchy
constraint appearing in (5b). By using binary variables to model sparsity in
the main/interaction effects and to encode the hierarchy constraint (5b),
Problem (5) can be expressed as the following MIP:

$\displaystyle\min_{\boldsymbol{\beta},\boldsymbol{\theta},\boldsymbol{z}}{}{}$
$\displaystyle
g_{\lambda_{1}}(\boldsymbol{\beta},\boldsymbol{\theta})+\lambda_{2}\big{(}\sum_{j\in[p]}{z_{j}}+\alpha\sum_{j<k}z_{j,k}\big{)}$
(6) s.t. $\displaystyle
z_{j},z_{j,k}\in\\{0,1\\},~{}~{}\left\lVert\boldsymbol{\beta}_{j}\right\rVert_{2}\leq
Mz_{j},~{}~{}\left\lVert\boldsymbol{\theta}_{j,k}\right\rVert_{2}\leq
Mz_{j,k}~{}~{}~{}\forall j<k,~{}j,k\in[p],$

where the BigM parameter $M$ is a sufficiently large finite constant such that
an optimal solution to (6) satisfies
$\max_{j}\|\boldsymbol{\beta}_{j}\|_{2}\leq M$ and
$\max_{j,k}\|\boldsymbol{\theta}_{j,k}\|_{2}\leq M$. Binary variable $z_{j}$
(as well as $z_{j,k}$) indicates whether the corresponding main effect
$\boldsymbol{\beta}_{j}$ (respectively, interaction effect
$\boldsymbol{\theta}_{j,k}$) is zero or not; the constraint appearing in the
last line of (6) enforces the hierarchy constraint in (5b).

To our knowledge, ours is the first work to study estimator (6). The
methodology presented here is of independent interest in the context of
structured nonparametric learning with interactions. In Section 2.4, we
discuss computational challenges associated with (5) and propose scalable
algorithms to obtain good solutions for this problem.

Related Work. Methodology for strong hierarchy in linear models has been
studied in the statistics/machine learning literature (Bien et al., 2013; Lim
and Hastie, 2015; Yan and Bien, 2017; Hazimeh and Mazumder, 2020b) – these
works appear to focus on the linear setting, which is more restrictive than
the nonlinear setting we consider here. Radchenko and James (2010) consider
hierarchy constraints in the nonparametric setting, via convex optimization
schemes. As far as we know, current techniques are unable to scale to the
functional learning instances we consider in our paper.

###  2.3 Statistical Theory

In this section, we explore statistical error bounds of our proposed estimator
in the deterministic design setting. Technical details, assumptions, and
proofs of the established results are provided in the supplementary material.

We write $\mathcal{C}_{\text{gr}}$ for the space of functions
$f(\mathbf{x})=\sum_{j\in[p]}f_{j}(x_{j})+\sum_{j<k}f_{j,k}(x_{j},x_{k})$,
where the components belong to $L_{2}$-Sobolev spaces, and the sparsity level
of the representation is bounded. We let $\Omega_{\text{gr}}(f)$ denote the
corresponding composite smoothness penalty, and write $G(f)$ for the number of
components in the above representation. We let $\|\cdot\|_{n}$ denote the
Euclidean norm divided by $\sqrt{n}$ and focus on the estimator that solves
the $\ell_{0}$-penalized version of optimization problem (2), that is

$\widehat{f}=\mathop{\rm arg\,min}_{f\in\mathcal{C}_{\text{gr}}}\;\,\|{\bf
y}-{\bf f}\|_{n}^{2}+\lambda_{n}\Omega_{\text{gr}}(f)+\mu_{n}G(f).$ (7)

We assume that the observed data follows the model $\mathbf{y}={\bf
f}^{*}+\boldsymbol{\epsilon}$, where $f^{*}$ is the unknown true regression
function of an arbitrary form, and the elements of $\boldsymbol{\epsilon}$ are
independent $N(0,\sigma^{2})$. We refer to $\|\widehat{\bf f}-{\bf
f}^{*}\|^{2}_{n}$ as the prediction error for estimator $\widehat{f}$. We
write $r_{n}=n^{-1/3}$, noting that $r_{n}^{2}$ is the optimal prediction
error rate in the bivariate regression setting where
$f^{*}\in\mathcal{C}_{\text{gr}}$. We say that a constant is universal if it
does not depend on other parameters, such as $n$ or $p$. We use the notation
$\lesssim$ to indicate that inequality $\leq$ holds up to a positive universal
multiplicative factor.

Theorem 1 presented below derives a general non-asymptotic oracle prediction
error bound, comparing the performance of our estimator to that of sparse
approximations to the true regression function. Bounds of this type have been
established for sparse nonparametric AMs with main effects; however, we are
not aware of similar existing bounds for nonparametric models with pairwise
interactions. The existing work has focused mainly on the Lasso-based
estimators, which encourage sparsity in the main-effects by penalizing the
magnitudes of the functional components (see Meier et al., 2009; Tan and
Zhang, 2019, and the references therein). The use of $\ell_{1}$-based
relaxations, in lieu of $\ell_{0}$-penalization, to induce sparsity in the
main-effects can lead to unwanted shrinkage, which may interfere with variable
selection. The approach in Hazimeh et al. (2021) controls the number of main-
effects directly, demonstrating the benefits of $\ell_{0}$-regularization both
theoretically and empirically. Our analysis accounts for the interaction
effects and focuses on the $\ell_{0}$-penalized formulation, which poses
additional theoretical challenges relative to the $\ell_{0}$-constrained
formulation. To our knowledge, the error bounds established here for sparse
nonparametric AMs with pairwise interactions are novel.

The following result corresponds to the high-dimensional setting where $p$ is
large (we discuss the classical fixed $p$ case at the end of this section),
hence the stated error bound holds with high probability.

######  Theorem 1.

There exists a universal $c_{1}$, such that if $\lambda_{n}\geq
c_{1}\sigma\big{[}r_{n}^{2}+r_{n}\sqrt{\log(ep)/n}\,\big{]}$, then

$\|\widehat{\bf f}-{\bf f}^{*}\|^{2}_{n}+\lambda_{n}\Omega_{\text{\rm
gr}}(\widehat{f})+\mu_{n}G(\widehat{f})\lesssim\inf_{f\in\mathcal{C}_{\text{gr}}}\Big{[}\|{\bf
f}-{\bf f}^{*}\|^{2}_{n}+\lambda_{n}\Omega_{\text{\rm
gr}}(f)+\mu_{n}G(f)\Big{]}+\sigma^{2}\Big{[}r_{n}^{2}+\frac{\log(ep)}{n}\Big{]}$

with probability at least $1-1/p$.

We note that a popular approach in the literature, even in the linear setting,
has been to assume a sparse underlying model for the data. In contrast,
Theorem 1 does not assume a ground truth sparse additive model, allowing for
model misspecification and incorporating the approximation error into the
bound. This feature is especially appealing for our application of interest
where the ground truth is not available. The following corollary complements
the general result in Theorem 1 by establishing a non-asymptotic prediction
error bound for the proposed approach in the case where the additive model is
correctly specified.

######  Corollary 1.

Suppose that $f^{*}\in\mathcal{C}_{\text{gr}}$. There exist universal
constants $c_{1}$ and $c_{2}$, such that if $\lambda_{n}\geq
c_{1}\sigma\big{[}r_{n}^{2}+r_{n}\sqrt{\log(ep)/n}\big{]}$ and $\mu_{n}\geq
c_{2}\sigma^{2}\big{[}r_{n}^{2}+\log(ep)/n\big{]}+c_{2}\lambda_{n}\Omega_{\text{\rm
gr}}(f^{*})$, then

$\|\widehat{\bf f}-{\bf
f}^{*}\|^{2}_{n}\lesssim\sigma^{2}\Big{[}r_{n}^{2}+\frac{\log(cp)}{n}\Big{]}+\lambda_{n}\Omega_{\text{\rm
gr}}(f^{*})+\mu_{n}G(f^{*})\quad\text{and}\quad G(\widehat{f})\leq G(f^{*})$
(8)

with probability at least $1-1/p$.

We note that inequality (8) yields the following prediction error bound:

$\|\widehat{\bf f}-{\bf
f}^{*}\|^{2}_{n}\lesssim\sigma^{2}\Big{[}n^{-2/3}+\frac{\log(ep)}{n}\Big{]}.$
(9)

In contrast to the Lasso-based estimators for nonparametric AMs (e.g., Meier
et al., 2009; Tan and Zhang, 2019), the above error rate holds without
imposing assumptions on the design. The accompanying sparsity bound
$G(\widehat{f})\leq G(f^{*})$ is exact, and hence stronger than the
corresponding bounds for the Lasso-based estimators (e.g., Lounici et al.,
2011). The latter bounds hold up to a multiplicative constant, which depends
on the design and can be significantly greater than one. When $\log(p)\lesssim
n^{1/3}$, the prediction error rate in (9) matches the optimal bivariate rate
of $n^{-2/3}$.

Theorem 1 and Corollary 1 are stated for the high-dimensional setting, where
$p$ is large. We show in the proof of Corollary 1 that in the classical
asymptotic setting, where $p$ is fixed and $n$ tends to infinity, an
appropriate choice of $\lambda_{n}$ and $\mu_{n}$ leads to the optimal
bivariate rate of convergence, $\|\widehat{\bf f}-{\bf
f}^{*}\|^{2}_{n}=O_{p}\big{(}n^{-2/3}\big{)}$, and the exact sparsity bound
$G(\widehat{f})\leq G(f^{*})$ that holds with probability tending to one.

###  2.4 Efficient Computations at Scale

We present specialized algorithms to obtain solutions to Problems (4) and (5).
Our approach scales to the problem-size of the Census application, with
$n\approx 10^{5}$ and $p\approx 500$, which poses formidable computational
challenges due to the presence of approximately $10^{5}$ interaction
effects999Note that, using 25 knots for every component, this leads to
estimating around 2.5 million basis coefficients.. To obtain good solutions at
scale, we use techniques inspired by first order methods in continuous
optimization Nesterov (2003) and a careful exploitation of the problem-
structure. A high-level summary is presented below, with the details relegated
to the supplement. For Problem (5) with hierarchical interactions, we propose
a relax-and-round based strategy (see Supplement for details) that appears to
work well in our applications.

Current state-of-art. To appreciate the computational challenges of sparse
nonlinear AMs with interactions, and nonparametric AMs in general, we provide
a few examples of the problem instances that can be handled by current state-
of-the-art algorithms with publicly available implementations. Current
implementations based on R package SAM (Zhao and Liu, 2012), the stepwise GAM
function in R package step.gam (which performs greedy variable selection), and
Python package pyGAM, take on the order of days to run and/or face numerical
difficulties for instances with $n\approx 10^{5}$ to obtain a single solution
without interactions. Furthermore, pyGAM is unable to do automated variable
selection. Wood et al. (2017) present an interesting approach for AMs that
scales to large $n$ settings (see R package mgcv) but doesn’t appear to
perform automated variable selection in the presence of a large number of
features (Wood et al. (2017) report instances containing fewer than $20$ pre-
specified main and interaction effects).

Algorithms for sparse nonlinear interactions: Problem (4). As mentioned
earlier, (4) can be expressed as a MIP and solved for small-to-moderate scale
problems using modern MIP solvers. However, given the problem-sizes of
interest and the fact that we seek to compute solutions to (4) for a family of
tuning parameters, we discuss faster alternatives. We note that the objective
in Problem (4) is a sum of a smooth convex loss function and a discontinuous
regularizer separable across the components $\\{\boldsymbol{\beta}_{j}\\}$ and
$\\{\boldsymbol{\theta}_{i,k}\\}$. Motivated by the strong empirical
performance of cyclical coordinate descent (CD) methods (Wright, 2015) in
$\ell_{0}$-penalized linear regression (Hazimeh and Mazumder, 2020a), we
explore block CD methods to obtain fast approximate solutions for the
nonparametric setting with interactions (4). For convergence guarantees of
this procedure, see Hazimeh and Mazumder (2020a) and references therein.

In our block CD method, the blocks correspond to the basis coefficients for
either the main effects $\\{\boldsymbol{\beta}_{j}\\}_{j}$ or the interaction
effects $\\{\boldsymbol{\theta}_{j,k}\\}_{j,k}$. Given an initialization
$(\boldsymbol{\beta}_{1}^{(0)},\cdots,\boldsymbol{\beta}_{p}^{(0)},\boldsymbol{\theta}_{1,2}^{(0)},\\\
\cdots,\boldsymbol{\theta}_{p-1,p}^{(0)})$, at every cycle, we sequentially
sweep across the main effects and the interaction effects. In particular, if
we denote the solution after $t$ cycles by
$(\boldsymbol{\beta}_{1}^{(t)},\cdots,\boldsymbol{\beta}_{p}^{(t)},\boldsymbol{\theta}_{1,2}^{(t)},\cdots,\boldsymbol{\theta}_{p-1,p}^{(t)})$,
then the block of coefficients for $j$-th main effect at cycle $t+1$, namely,
$\boldsymbol{\beta}_{j}^{(t+1)}$, is obtained by optimizing (4) with respect
to $\boldsymbol{\beta}_{j}$, with other variables held fixed:

$\displaystyle\boldsymbol{\beta}_{j}^{(t+1)}\in\underset{\boldsymbol{\beta}_{j}\in\mathbb{R}^{K_{j}}}{\mathop{\rm
arg\,min}}~{}G(\boldsymbol{\beta}_{1}^{(t+1)},\cdots,\boldsymbol{\beta}_{j-1}^{(t+1)},\boldsymbol{\beta}_{j},\boldsymbol{\beta}_{j+1}^{(t)},\cdots,\boldsymbol{\beta}_{p}^{(t)},~{}~{}\boldsymbol{\theta}_{1,2}^{(t)},\cdots,\boldsymbol{\theta}_{p-1,p}^{(t)}).$
(10)

We update $\boldsymbol{\theta}_{j,k}^{(t+1)}$, containing the coefficients for
the $(j,k)$-th interaction at cycle $t+1$, using

$\displaystyle\boldsymbol{\theta}_{j,k}^{(t+1)}\in\underset{\boldsymbol{\theta}_{j,k}\in\mathbb{R}^{K_{j,k}}}{\mathop{\rm
arg\,min}}~{}G(\boldsymbol{\beta}_{1}^{(t+1)},\cdots,\boldsymbol{\beta}_{p}^{(t+1)},~{}~{}\boldsymbol{\theta}_{1,2}^{(t+1)},\cdots,\boldsymbol{\theta}_{(j,k)-1}^{(t+1)},\boldsymbol{\theta}_{j,k},\boldsymbol{\theta}_{(j,k)+1}^{(t)}\cdots,\boldsymbol{\theta}_{p-1,p}^{(t)}).$
(11)

The block minimization problems (10) and (11) can be solved in closed-form as
discussed in Supplement. These block CD updates need to be paired with several
computational devices in the form of active set updates, cached matrix
factorizations, and warm-starts, among others. We draw inspiration from
similar strategies used in CD-based procedures for sparse linear regression
(Hazimeh and Mazumder, 2020a; Friedman et al., 2010), and extend them to our
problem. See Supplement for more details.

##  3 Case study: Predicting the Census Survey Self-Response Rate

We now discuss the Census application, which is the main motivation for this
paper. As mentioned in Section 1, obtaining interpretable models with good
predictive capabilities is critically important in this application context –
such models have the potential to directly influence the planning of outreach
campaigns and enable stakeholders to make well-informed decisions on the
allocation of spending for different communications channels (for example, TV,
radio, digital), advertising messages with properly tailored content, and the
timing of spending during the campaigns (Kulzick et al., 2019). As discussed
earlier, the black-box predictive models that topped the nationwide
competition were not actionable due to limited interpretability. On the other
hand, the linear models that drive the LRS have limited predictive power. In
this section we discuss how our proposed approach appears to strike a nice
balance between interpretability and good prediction performance.

The data we use for this study is publicly available in the US Census Planning
Database, which provides a range of demographic, socioeconomic, housing, and
Census operational data (US Census Bureau, 2019). The data in the Planning
Database includes covariates from the 2010 Census and the 2013-2017 American
Community Survey (ACS), aggregated at both the Census tract level and the
block level. We use tract level data, with approximately $74,000$ observations
and $500$ covariates. The response is the ACS self-response rate. We exclude
the following covariates from our model for interpretability reasons: spatial
covariates (“State”, “County”, “Tract”, “Flag”, “AIAN Land”) and variables
that serve as a proxy to the response (for example, “Low response score”,
“Number of housing units that returned first forms”, “Replacement forms” or
“Bilingual forms” in Census 2010). We also remove the margin of error
variables corresponding to the ACS. After excluding these variables, we are
left with $p=295$ covariates101010To clarify, the results presented in this
section are based on these $295$ features, though our approach scales to
$p=500$ features..

###  3.1 Experimental setup

We randomly split the data into $58$K observations for training, $7.2$K
observations for validation (used for selecting tuning parameters), and $7.2$K
observations for testing. The features were standardized to have zero mean and
unit variances. We use the squared $\ell_{2}$ loss for training and evaluate
the performance of the models in terms of both root mean square (RMSE) and
mean absolute error (MAE). In addition, we consider variable selection
performance – for example, with respect to the number of features retained and
the associated interpretations they offer.

As noted in Erdman and Bates (2016), the current approach for predicting the
self-response rates–i.e., the LRS–is based on a linear regression model with
around $25$ hand-selected features. Erdman and Bates (2016) also discuss using
tailored variable transformations on some of the features to incorporate
nonlinear effects. Such a manual approach is impractical in our situation, as
we are dealing with a large number of features along with the corresponding
pairwise interactions. However, as we use nonparametric models, which seek to
learn nonlinearities in the main and interaction effects, such feature
engineering may not be necessary.

Benchmark Methods. We evaluate and compare our models against existing linear
and nonparametric approaches. We consider the following linear models with
main effects only (i.e., without interactions): (i) Ridge regression; (ii)
Lasso regression (Hastie et al., 2015); (iii) $\ell_{0}$-penalized regression
with additional ridge regularization (denoted as $\ell_{0}-\ell_{2}$),
implemented via L0Learn (Hazimeh and Mazumder, 2020a). In addition, we
consider the following linear models with interaction effects: (iv) Lasso with
main effects and all pairwise linear interactions, (v) hierScale: a convex
optimization framework for learning sparse main-effects and interactions under
a strong hierarchy constraint (Hazimeh and Mazumder, 2020b). For (i), (ii),
and (iv), we use Python’s scikit-learn library (Pedregosa et al., 2011). In
addition to the linear model benchmarks, we also compare against state-of-the-
art black-box machine learning methods such as (vi) Gradient boosted decision
trees (GBDT) from XGBoost (Chen and Guestrin, 2016); and (vii) Feedforward
neural networks111111We tune both of the nonparametric models (vi) and (vii)
with random search procedures using Python package hyperopt (Bergstra et al.,
2015) for hyperparameter optimization. GBDT is tuned with respect to the
maximum depth, number of estimators, learning rate, and L2 regularization.
Neural networks are tuned with respect to the number of dense layers, number
of units in dense layers, choice and strength of regularizers ($\ell_{1}$,
$\ell_{2}$, dropout), learning rates for Adam optimizer (Kingma and Ba, 2015),
and batch sizes. The number of tuning parameters for all the nonparametric
models is capped at 1000..

Proposed Models. We consider the following estimators: (a)
$\ell_{0}$-penalized AMs with nonlinear main effects. We also consider the
following AMs containing both (nonlinear) main and pairwise interaction
effects: (b) $\ell_{0}$-penalized AMs with both main and interaction effects,
i.e., estimator (4); and (c) sparse hierarchical interactions, i.e., estimator
(5). For all our methods, we use 1000 values of the tuning parameters.

All our algorithms for estimators (a)–(c) above are implemented in Python. We
use cubic B-splines with $10$ knots for the main effects. For the interaction
effects, we use tensor spline bases of degree $3$ with 5 knots in each
coordinate, leading to a total of $5\times 5=25$ knots. Because the problem at
hand has $295$ covariates and $43,365$ possible pairwise interactions, we need
to be careful with implementation aspects while generating spline-transformed
representations for all the interaction effects, which can be memory
intensive.

For method (b), i.e., estimator (4), we set $\alpha=1$, so that both the main
and the interaction effects have the same $\ell_{0}$-penalty parameter. We
note, however, that the strong hierarchy constraint in (5) implicitly gives
higher priority to the main effects as compared to the interaction effects.
For estimator (c), i.e., Problem (5), we fix the smoothing parameter
$\lambda_{1}$ to the optimal value available from (b), to reduce tuning.
$\tau$ is chosen based on validation tuning.

###  3.2 Comparing methods: prediction, sparsity and interpretability

Table 1: _Comparisons of our methods with several benchmark models as
discussed in the text. We display the tuned RMSE and MAE values for the
different models, along with the corresponding number of covariates. LMs and
AMs that include pairwise interaction effects are denoted by
“..+Interactions”. The shorthand 16 Mn + 174 Int (for example) refers to 16
main effects and 174 interaction effects in the selected model, this model
contains 160 effective features. Best metrics are highlighted in bold. Numbers
within parentheses provide the standard errors._

Type | Model | RMSE | MAE | #Covariates
---|---|---|---|---
Linear Models (LM) | Ridge (Model (i)) | $6.804$ | $5.254$ | 295
Lasso (Model (ii)) | $6.803$ | $5.254$ | 221
L0Learn ($\ell_{0}-\ell_{2}$) (Model (iii)) | $6.813$ | $5.268$ | 136
| LM+Interactions (Lasso) (Model (iv)) | $6.528$ | $5.026$ | 264
|  |  |  | (76 Mn + 1598 Int)
| LM+Interactions with | $6.621$ | $5.086$ | 276
| Strong Hierarchy (hierScale) (Model (v)) |  |  | (276 Mn + 4885 Int)
Nonparametric Additive Models (AM) (this paper) |  $\ell_{0}$-sparse AM (Model (a)) | $6.593$ | $5.120$ | 182
$\ell_{0}$-sparse AM+Interactions | $6.467$ | $\mathbf{4.973}$ | 160
(Model (b)) |  |  | (16 Mn + 174 Int)
AM+Interactions with | $\mathbf{6.452}$ | $4.995$ | 131
|  Strong Hierarchy (Model (c)) |  |  | (131 Mn + 173 Int)
Nonparametric (Non-interpretable) | GBDT (Model (vi)) | $6.472$ | $4.983$ | 295
Feedforward Neural Networks (Model (vii)) | $6.501$ | $4.996$ | 295

We discuss the performance of different estimators in terms of prediction
accuracy, model parsimony and interpretability.

Additive nonlinear models vs black-box ML methods. Table 1 reports the
prediction errors (RMSE and MAE) on the test-set along with the number of
nonzero features in the model. Importantly, we observe that the test
performance of AMs with sparse interactions (both with and without the
hierarchy constraints) is comparable to that of the best black box predictive
ML models (vi)-(vii) (GBDT, Neural Networks), when one takes into account the
model standard errors. Model (c) i.e., nonparametric AMs with interactions
with hierarchy, appears to be the winner in terms of the raw RMSE values,
closely followed by model (b), i.e., $\ell_{0}$-sparse nonparametric AMs with
interactions without hierarchy (4). In contrast, model (b) appears to be the
winner in terms of raw MAE values, closely followed by GBDT and model (c).
Interpretability (i.e., model sparsity and/or hierarchy) is the key
differentiating factor among the leading prediction methods. Model (b), i.e.,
$\ell_{0}$-sparse AMs with interactions, selects fewer effects: 16 main (Mn)
and 174 interactions (Int) when compared to model (c), i.e., AMs with sparse
hierarchical interactions (131Mn+173Int). However, recall that model (b) does
not obey the hierarchy constraint. On the other hand, model (c), i.e., AMs
with sparse interactions and hierarchy, selects fewer effective features than
model (b), which is important for cost-savings from an operational standpoint.

The black-box ML methods offer limited interpretability when compared to the
AMs proposed here (additional discussions are provided further below). While
Table 1 presents a summary of the best models chosen based on the validation
set prediction performance, it is desirable to prioritize models with fewer
features (perhaps at the cost of a marginal deterioration in predictive
performance)—see Section 3.3 for further discussion.

Nonlinear models vs linear models. Table 1 suggests that
nonlinear/nonparametric models have better predictive performance than their
linear counterparts. Also, linear models with sparse main and interaction-
effects appear to have an edge over sparse linear models without interactions.
Another appealing aspect of sparse nonparametric AMs compared to linear
models, is in model parsimony—we alluded to this aspect in Figure 3 with a
reduced number of features. We observe that the number of nonzero effects is
significantly lower for the nonlinear AMs. For example, $\ell_{0}$-sparse
nonparametric AMs with no interactions (i.e., Model (a)) can achieve the same
level of predictive performance as its linear model counterpart (i.e., Model
(iii)) with significantly fewer covariates: $40$ versus $136$. Similarly, if
we compare hierScale (i.e., linear main and interaction effects with strong
hierarchical restrictions) to its nonparametric counterpart (5), the number of
main effects (and also, covariates) reduces from 276 to 131, and the number of
interaction effects reduces from 4,885 to 173. Interestingly, all the
covariates selected by the nonparametric AMs with interaction models
(with/without strong hierarchy) are contained in the set of covariates
selected by the sparse linear models with both main and interaction effects.

Additive nonlinear models: interactions vs no interactions. Table 1 shows that
the number of effective features can be reduced by using nonlinear pairwise
interactions in the nonparametric AMs: Comparing performance of Model (a) and
Model (c) we observe that the number of features reduce from $182$ to $131$ by
taking into account interaction effects obeying a strong hierarchy principle.
This difference is possibly because a large number of nonlinear main effects
attempt to explain the nonlinear interaction effects. The resulting
compression points to the redundancy of some of the covariates when
interaction effects are directly included in the model. Furthermore, the
prediction performance of nonlinear AMs seems to improve when pairwise
interactions are included into the model.

Insights from visualizations. As illustrated by Figure 2, which displays some
marginal nonlinear fits, an appealing aspect of nonlinear AMs is that they
naturally allow the practitioner to gather insights into the associations
between the response and a feature of interest, for example, by exploring the
map $x_{j}\mapsto f_{j}(x_{j})$. Similarly, the map $(x_{j},x_{k})\mapsto
f_{j,k}(x_{j},x_{k})$ would shed light into how the interaction of
$(x_{j},x_{k})$ influences the output. Such interpretations are a by-product
of our additive model framework; and may not be readily available via black-
box ML methods such as Neural Networks and GBDT. These association plots, for
example, can help stakeholders identify promising factors and patterns
influencing self-response scores—these insights can in turn inform policy
decisions regarding targeted investments and optimizing the operational costs.

Figure 4: Predicted ACS self-response rates for all tracts in the United
States.

.

Figure 5: _ACS self-response rates for all tracts in the District of Columbia.
a) Actual ACS self-response rates. b) Predicted ACS self-response rates for
AMs with interactions ( 4). c) Difference between the actual and predicted
self-response rate: difference = actual - predicted._

To obtain a finer understanding of the performance of our models, we use
visualization tools inspired by earlier works from the US Census Bureau (US
Census Bureau, 2017; Kulzick et al., 2019). Figure 4 illustrates the tract
self-response rates predicted by our $\ell_{0}$-sparse nonparametric AMs with
interactions model (4) on a map of the United States. Different from Figure 1,
which shows actual data with missing values for some of the tracts121212Some
responses are deliberately suppressed by the US Census Bureau for privacy
considerations, to limit the disclosure of information about individual
respondents and to reduce the number of estimates with unacceptable levels of
statistical reliability (US Census Bureau, 2016), Figure 4 provides predicted
self-response rates for all tracts based on our proposed model. Predictive
models in general and our models in particular, allow us to forecast the
response rates for tracts where gathered data is incomplete. It is worth
noting that our regression-based models may be preferable over commonly used
nearest neighbor type methods: the latter models do not generally offer
insights into factors associated with high/low response rates. Nearest
neighbor based methods may also be ill-posed for isolated regions such as
Alaska, Hawaii and Puerto Rico, each having at least $10\%$ of their tracts
censored.

Figure 5 displays, side by side, the actual and the predicted self-response
rates for the tracts in Washington DC, as well as the corresponding
differences. Both the actual and the predicted rates are higher than average
in most parts of the Northwest and a portion of the Northeast DC, and lower in
the Southeast and most of the Northeast DC.

Figure 6 shows how the sorting of the Census tracts into quintiles (US Census
Bureau, 2017; Kulzick et al., 2019) of the actual self-response rates compares
with the corresponding sorting of the predictions made by model (I):
$\ell_{0}$-sparse nonparametric AMs with main and pairwise interaction effects
(4), and model (II): the Lasso (we use an $\ell_{1}$-penalized linear model
with main effects and no interactions). For example, the top panel in Figure 6
shows that among the tracts in the first quintile of the actual self-response
rates, $85.5\%$ are correctly predicted by model (I) to fall in the first
quintile. The corresponding proportion for model (II) is $84.4\%$. The same
panel shows that among the tracts in the first quintile of the actual self-
response rates, $13.9\%$ are incorrectly predicted to fall in the second
quintile; this proportion is smaller than the one for the Lasso ($14.7\%$).
Similarly, in the second quintile of the actual self-response rates, a higher
proportion ($59.4\%$) is correctly identified by model (I) to be in the second
quintile; this is again an improvement over the $56.6\%$ for the Lasso
regression model. For all 5 quintiles of the actual self response rates, model
(I) identifies a higher proportion to be in the correct quintile than model
(II)—this suggests that our approach appears to have an edge over the Lasso.

Figure 6: Comparison of quintile groups based on Lasso i.e., model (I) (US
Census Bureau, 2017; Kulzick et al., 2019) and our proposed $\ell_{0}$-sparse
AMs with main-effects and interactions i.e., model (II).

.

Summary. Our proposed framework – AMs with sparse nonlinear main and
interaction effects – allows for automated variable selection, provides simple
interpretability, and can deliver useful insights into the underlying
associations between the covariates and the response. Given its strong
predictive performance relative to the black-box ML methods, our proposed
approach has the potential to be a useful practical tool for predicting self-
response rates.

###  3.3 Interpreting important features

We now illustrate how our methodology can guide the practitioner in obtaining
a set of interpretable features and deriving associated actionable insights
into the factors contributing towards low response-rates in surveys. An
important aspect of our regularized learning framework is that it provides an
automated scheme to identify a collection of models, balancing out the
complexity of the model and data-fidelity. To see this, Figure 7[left panel]
plots the number of main and interaction effects against the associated
prediction error for model (4). The plot illustrates that by trading off a
little in the predictive performance (in terms of RMSE, which is shown by the
red dashed line), we can limit the number of effects at any level, as desired
by the practitioner who intends to obtain a more parsimonious representation
of the factors impacting survey response rates. Specifically, if we would like
to limit the number of main effects to under $20$, Figure 7(b) shows the top
$19$ main effects in the order they enter the model along the regularization
path. The definitions of the variables in Figure 7(b) are provided in the
Supplement. Our investigation reveals that most (though not all) of these
variables also appear in the models reported in prior studies by the Census
Bureau—see for example Erdman and Bates (2016) and Kulzick et al. (2019).
Interestingly, some of the features discovered by our framework can be
interpreted in the context of clusters defined in Bates and Mulry (2011) or
mindset solutions in Kulzick et al. (2019), as discussed below.

|
---|---
Figure 7: _Estimates available from model ( 4): sparse AMs with main effects
and interactions. [Left] number of nonzero effects vs the corresponding RMSE.
[Right] Variables corresponding to main effects selected by the model along
the regularization path. For visualization purposes, we focus on the top 19
main effects._

Using cluster analysis on the 12 covariates from the HTC score, Bates and
Mulry (2011) categorize the tracts into eight distinct clusters: All Around
Average I and II, Economically Disadvantaged I and II, Ethnic Enclave I and
II, Single Unattached Mobiles, and Advantaged Homeowners. From the cluster
analysis in Bates and Mulry (2011), we see that the reported clusters have a
wide variation in response rate (the response rate is not used to identify the
clusters). Each cluster has some key characteristics (pertaining to a subset
of variables) that seem to drive the response rate. From the description of
these clusters in Bates and Mulry (2011), we see that the covariates
automatically discovered by our framework, as listed in Figure 7[right panel],
broadly cover the important characteristics upon which the clusters are based.
For example, there is high variability across clusters with respect to the
occupancy rate and home ownership status, which are captured by covariates
“Owner occupied housing units” and “Total occupied units”. Racial diversity
also appears to play a role in cluster formation—in particular, the covariates
related to ethnic groups characterize the Ethnic Enclave clusters. Similarly,
covariates related to poverty, for example, “Persons below poverty line”, lack
of education (“No high school diploma”), and lack of health insurance (“No
health insurance”), determine the Economically Disadvantaged clusters.
Covariates related to high mobility, for example, “Different housing unit in
prior year”, high education (“College or higher degree holders”), marital
status (“Single person households”), and nonspousal occupants (“Unmarried
couples”), characterize the Single Unattached Mobiles cluster.

The above discussion suggests that there may be important interaction effects
within the variable groups that characterize each cluster. We observe that our
AMs with interactions model (4) automatically selects many of these effects.
For instance, our model recovers several interaction effects between important
covariates characterizing the Economically Disadvantaged I cluster. Bates and
Mulry (2011) notes that this cluster has high percentages of people in
poverty, people receiving public assistance, and adults without a high school
education. Our model identifies an interaction effect between “No high school
graduation” and “No health insurance”. The same cluster also has a large
African American population and above average number of children under 18. Our
model potentially captures this characteristic by selecting an interaction
effect between “Non-Hispanic blacks” and “Population with ages 5-17”. We also
recover an interaction effect between “Non-Hispanic blacks” and “No high
school graduation”, which seems to be reflective of this cluster. Consider
also the Advantaged Homeowners cluster characterized by stable homeowners who
are predominantly white and live in single-unit spousal-households. Our model
automatically selects a closely related interaction effect between “Non-
Hispanic whites” and “Single units”. Another example is the Single Unattached
Mobiles cluster, which is characterized by young, highly educated population
living in multi-units structures with high mobility (Bates and Mulry, 2011);
the tracts are almost exclusively urban and have an above-average Asian
population relative to the national average. Our model selects multiple
interaction effects within this group of variables, for example: i) “Moved in
2010 or later” and “Multi-unit structures”; ii) “College or higher degree
holders” and “Single person households”; iii) ”People speaking Asian
languages” and “Multi-unit structures”; iv) “People speaking Asian languages”
and “Urban clusters”.

Some of the covariates in Figure 7[right] can also be interpreted in terms of
the mindset solutions in Kulzick et al. (2019). In a mindset solution, survey
respondents are categorized using factor analysis into six mindsets based on
their predicted self-response patterns and demographic characteristics. These
mindsets are determined by the demographic characteristics such as age, race,
income, home-ownership, presence of children in the household, marital status,
internet use, English proficiency, and country of birth. The six mindsets are:
Eager Engagers, Fence Sitters, Individuals with a Confidentiality, Head
Nodders, Wary Skeptics and Disconnected Doubters. The covariates can be also
interpreted in terms of eight segmented tracts: Responsive Suburbia, Main
Street Middle, Country Roads, Downtown Dynamic, Student and Military
Communities, Sparse Spaces, Multicultural Mosaic, and Rural Delta and Urban
Enclaves – see Kulzick et al. (2019) for additional details.

###  3.4 Conclusions

Our goal has been to create a framework to help practitioners identify key
factors that influence survey response rates. The proposed approach leads to
strong predictive performance while simultaneously delivering parsimonious and
interpretable statistical models. Insights gathered from our models may be
used for targeted follow-up surveys and allocation of resources for
advertisements. We hope that our models will facilitate the goal of achieving
improved census coverage in a timely fashion, while optimizing the operational
risks and costs. Our framework can complement and potentially offer an
alternative to the LRS metric that is currently used in the US Census Bureau
ROAM application.

### Acknowledgements

The authors would like to thank Nancy Bates, Tommy Wright and Joanna Fane
Lineback from the U.S. Census Bureau for helpful discussions. Thanks also to
other researchers from the Centers for Statistical Research and Methodology
and Behavioral Science Methods (U.S Census Bureau) for helpful feedback in
earlier versions of this work. Any opinions and conclusions expressed herein
are those of the authors and do not necessarily represent the views of the
U.S. Census Bureau. This research was supported in part, by grants from the
Office of Naval Research: ONR-N000141812298 (YIP) and the National Science
Foundation: NSF-IIS-1718258 awarded to Rahul Mazumder.

## References

  * Erdman and Bates [2016] Chandra Erdman and Nancy Bates.  The low response score (LRS).  _Public Opinion Quarterly_ , 81(1):144–156, December 2016.  doi: 10.1093/poq/nfw040.
  * Smith [2014] Tom W. Smith.  Hard-to-survey populations in comparative perspective.  In _Hard-to-Survey Populations_ , pages 21–36. Cambridge University Press, Cambridge, 2014.  ISBN 9781139381635.  doi: 10.1017/CBO9781139381635.004.
  * Tourangeau et al. [2014] Roger Tourangeau, Brad Edwards, Timothy P Johnson, et al.  _Hard-to-survey populations_.  Cambridge University Press, 2014.
  * US Census Bureau [2009] US Census Bureau.  2010 census integrated communications campaign plan: The success of the census is in our hands.  November 2009.
  * Bates and Mulry [2011] Nancy Bates and Mary H. Mulry.  Using a geographic segmentation to understand, predict, and plan for census and survey mail nonresponse.  _Journal of Official Statistics_ , 27(4):601–618, July 2011.
  * Mule [2012] Tom Mule.  2010 census coverage measurement estimation report: Summary of estimates of coverage for persons in the united states.  Dssd 2010 census coverage measurement memorandum series no. 2010-g-01, US Census Bureau, May 2012.
  * Bruce et al. [2001] Antonio Bruce, J Gregory Robinson, and Monique V Sanders.  Hard-to-count scores and broad demographic groups associated with patterns of response rates in census 2000.  In _Joint Statistical Meeting_ , 2001.
  * Hastie and Tibshirani [1987] Trevor Hastie and Robert Tibshirani.  Generalized additive models: some applications.  _Journal of the American Statistical Association_ , 82(398):371–386, 1987.
  * Hastie et al. [2001] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  _The Elements of Statistical Learning_.  Springer Series in Statistics. Springer New York Inc., New York, NY, USA, 2001.
  * Bach et al. [2012] Francis Bach, Rodolphe Jenatton, Julien Mairal, and Guillaume Obozinski.  Structured sparsity through convex optimization.  _Statistical Science_ , 27(4):450–468, November 2012.  doi: 10.1214/12-sts394.
  * Bertsimas et al. [2016] Dimitris Bertsimas, Angela King, and Rahul Mazumder.  Best subset selection via a modern optimization lens.  _Annals of Statistics_ , 44(2):813–852, 2016.
  * Mazumder and Radchenko [2017] Rahul Mazumder and Peter Radchenko.  The Discrete Dantzig Selector: Estimating sparse linear models via mixed integer linear optimization.  _IEEE Transactions on Information Theory_ , 63 (5):3053 – 3075, 2017.
  * Mazumder et al. [2017] Rahul Mazumder, Peter Radchenko, and Antoine Dedieu.  Subset selection with shrinkage: Sparse linear modeling when the snr is low.  _arXiv preprint arXiv:1708.03288 (Operations Research, to appear)_ , 2017.
  * Hazimeh and Mazumder [2020a] Hussein Hazimeh and Rahul Mazumder.  Fast best subset selection: Coordinate descent and local combinatorial optimization algorithms.  _Oper. Res._ , 68(5):1517–1537, September 2020a.  ISSN 0030-364X.  doi: 10.1287/opre.2019.1919.
  * Hazimeh and Mazumder [2020b] Hussein Hazimeh and Rahul Mazumder.  Learning hierarchical interactions at scale: A convex optimization approach.  In Silvia Chiappa and Roberto Calandra, editors, _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , volume 108 of _Proceedings of Machine Learning Research_ , pages 1833–1843. PMLR, 26–28 Aug 2020b.
  * Bien et al. [2013] Jacob Bien, Jonathan Taylor, and Robert Tibshirani.  A lasso for hierarchical interactions.  _The Annals of Statistics_ , 41(3):1111–1141, June 2013.  doi: 10.1214/13-aos1096.
  * Kulzick et al. [2019] Robert Kulzick, Laura Kail, Shawnna Mullenax, et al.  2020 census predictive models and audience segmentation report.  June 2019.
  * Meier et al. [2009] Lukas Meier, Sara van de Geer, and Peter Buhlmann.  High-dimensional additive modeling.  _The Annals of Statistics_ , 37(6B):3779 – 3821, 2009.  doi: 10.1214/09-AOS692.
  * Ravikumar et al. [2009] Pradeep Ravikumar, John Lafferty, Han Liu, et al.  Sparse additive models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 71(5):1009–1030, November 2009.  doi: 10.1111/j.1467-9868.2009.00718.x.
  * Huang et al. [2010] J. Huang, J.L. Horowitz, and F. Wei.  Variable selection in nonparametric additive models.  _The Annals of Statistics_ , 38:2282–2313, 2010.
  * Zhao and Liu [2012] Tuo Zhao and Han Liu.  Sparse additive machine.  In Neil D. Lawrence and Mark Girolami, editors, _Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics_ , volume 22 of _Proceedings of Machine Learning Research_ , pages 1435–1443, La Palma, Canary Islands, 21–23 Apr 2012. PMLR.
  * Yuan and Zhou [2016] Ming Yuan and Ding-Xuan Zhou.  Minimax optimal rates of estimation in high dimensional additive models.  _The Annals of Statistics_ , 44(6):2564–2593, 2016.
  * Hazimeh et al. [2021] Hussein Hazimeh, Rahul Mazumder, and Peter Radchenko.  Grouped variable selection with discrete optimization: Computational and statistical perspectives.  _arXiv preprint arXiv:2104.07084 (Annals of Statistics, to appear)_ , 2021.
  * Wahba [1990] Grace Wahba.  _Spline Models for Observational Data_.  Society for Industrial and Applied Mathematics, January 1990.  doi: 10.1137/1.9781611970128.
  * Stone [1986] Charles J Stone.  The dimensionality reduction principle for generalized additive models.  _The Annals of Statistics_ , pages 590–606, 1986.
  * Hastie et al. [2015] Trevor Hastie, Robert Tibshirani, and Martin Wainwright.  _Statistical Learning with Sparsity: The Lasso and Generalizations_.  Chapman & Hall/CRC, 2015.  ISBN 1498712169, 9781498712163.
  * Radchenko and James [2010] Peter Radchenko and Gareth M. James.  Variable selection using adaptive nonlinear interaction structures in high dimensions.  _Journal of the American Statistical Association_ , 105(492):1541–1553, December 2010.  doi: 10.1198/jasa.2010.tm10130.
  * Erdman and Bates [2014] C. Erdman and Nancy J. Bates.  The u.s. census bureau mail return rate challenge: Crowdsourcing to develop a hard-to-count score.  2014\.
  * US Census Bureau [2019] US Census Bureau.  Planning database, Jun 2019.
  * Wolsey and Nemhauser [1999] Laurence A Wolsey and George L Nemhauser.  _Integer and combinatorial optimization_ , volume 55.  John Wiley & Sons, 1999.
  * Hazimeh et al. [2020] Hussein Hazimeh, Rahul Mazumder, and Ali Saab.  Sparse regression at scale: Branch-and-bound rooted in first-order optimization.  _arXiv preprint arXiv:2004.06152_ , 2020.
  * Buhlmann and van de Geer [2011] Peter Buhlmann and Sara van de Geer.  _Statistics for High-Dimensional Data: Methods, Theory and Applications_.  Springer Publishing Company, Incorporated, 1st edition, 2011.  ISBN 3642201911.
  * McCullagh and Nelder [1989] Peter McCullagh and John A Nelder.  _Generalized linear models_ , volume 37.  CRC press, 1989.
  * Lim and Hastie [2015] Michael Lim and Trevor Hastie.  Learning interactions via hierarchical group-lasso regularization.  _Journal of Computational and Graphical Statistics_ , 24(3):627–654, July 2015.  doi: 10.1080/10618600.2014.938812.
  * Yan and Bien [2017] Xiaohan Yan and Jacob Bien.  Hierarchical sparse modeling: A choice of two group lasso formulations.  _Statistical Science_ , 32(4):531–560, November 2017.  doi: 10.1214/17-sts622.
  * Tan and Zhang [2019] Zhiqiang Tan and Cun-Hui Zhang.  Doubly penalized estimation in additive regression with high-dimensional data.  _The Annals of Statistics_ , 47(5):2567–2600, 2019.
  * Lounici et al. [2011] K. Lounici, M. Pontil, S. van de Geer, et al.  Oracle inequalities and optimal inference under group sparsity.  _The Annals of Statistics_ , 39(4):2164–2204, 2011.
  * Nesterov [2003] Yurii Nesterov.  _Introductory lectures on convex optimization: A basic course_ , volume 87.  Springer Science & Business Media, 2003.
  * Wood et al. [2017] Simon N. Wood, Zheyuan Li, Gavin Shaddick, et al.  Generalized additive models for gigadata: Modeling the u.k. black smoke network daily data.  _Journal of the American Statistical Association_ , 112(519):1199–1210, April 2017.  doi: 10.1080/01621459.2016.1195744.
  * Wright [2015] Stephen J Wright.  Coordinate descent algorithms.  _Mathematical Programming_ , 151(1):3–34, 2015\.
  * Friedman et al. [2010] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of Statistical Software_ , 33(1), 2010.  doi: 10.18637/jss.v033.i01.
  * Pedregosa et al. [2011] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, et al.  Scikit-learn: Machine learning in python.  _J. Mach. Learn. Res._ , 12(null):2825–2830, November 2011.  ISSN 1532-4435.
  * Chen and Guestrin [2016] Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.  ISBN 978-1-4503-4232-2.  doi: 10.1145/2939672.2939785.
  * Bergstra et al. [2015] James Bergstra, Brent Komer, Chris Eliasmith, et al.  Hyperopt: a python library for model selection and hyperparameter optimization.  _Computational Science & Discovery_, 8(1):014008, 2015.
  * Kingma and Ba [2015] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.
  * US Census Bureau [2017] US Census Bureau.  Foundational communications models 1.0: the updated low response score (lrs) and an internet lrs.  2017\.
  * US Census Bureau [2016] US Census Bureau.  American community survey data suppression.  September 2016.
  * Hall and Opsomer [2005] Peter Hall and J. D. Opsomer.  Theory for penalised spline regression.  _Biometrika_ , 92(1):105–118, March 2005.  doi: 10.1093/biomet/92.1.105.
  * Li and Ruppert [2008] Yingxing Li and David Ruppert.  On the asymptotics of penalized splines.  _Biometrika_ , 95(2):415–436, 2008.  ISSN 00063444.
  * Wood [2006] Simon N. Wood.  Low-rank scale-invariant tensor product smooths for generalized additive mixed models.  _Biometrics_ , 62(4):1025–1036, May 2006.  doi: 10.1111/j.1541-0420.2006.00574.x.
  * Eilers and Marx [1996] Paul H. C. Eilers and Brian D. Marx.  Flexible smoothing with b-splines and penalties.  _Statistical Science_ , 11:89–121, 1996.
  * Kammann and Wand [2003] E. E. Kammann and M. P. Wand.  Geoadditive models.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 52(1):1–18, January 2003.  doi: 10.1111/1467-9876.00385.
  * Eilers and Marx [2003] Paul H.C. Eilers and Brian D. Marx.  Multivariate calibration with temperature interaction using two-dimensional penalized signal regression.  _Chemometrics and Intelligent Laboratory Systems_ , 66(2):159–174, June 2003.  doi: 10.1016/s0169-7439(03)00029-7.
  * Lin and Zhang [2006] Y. Lin and H. H. Zhang.  Component selection and smoothing in multivariate nonparametric regression.  _The Annals of Statistics_ , 34:2272–2297, 2006.
  * Oden and Reddy [1976] John Tinsley Oden and Junuthula Narasimha Reddy.  _An introduction to the mathematical theory of finite elements_.  Wiley, New York, 1976.
  * Van de Geer [2000] Sara Van de Geer.  _Empirical Processes in M-Estimation_.  Cambridge University Press, Cambridge, 2000.
  * Agmon [1965] S Agmon.  _Lectures on Elliptic Boundary Value Problems_.  Van Nostrand, Princeton, NJ, 1965.
  * Birman and Solomjak [1967] M. S. Birman and M. Z. Solomjak.  Piecewise-polynomial approximations of functions of the classes $w^{\alpha}_{p}$.  _Math. USSR-Sbornik_ , 2(3):295–317, 1967.

## Supplementary Material

This supplemental document contains some additional results not included in
the main body of the paper. In particular, we present

  1. -

a detailed discussion of the univariate and bivariate smooth function
estimation with splines. These serve as building blocks for modelling main and
interaction effects in Section 2.1. The discussion is outlined in Section S1.

  2. -

scalability considerations for solving optimization problem (4) in Section
2.2.

  3. -

technical details and proofs for Theorem 1 in Section 2.3.

  4. -

Census variable descriptions for Figure 3 in Section 2.2 and Figure 7b in
Section 3.3.

This supplement is not entirely self-contained, so readers might need to refer
to the main paper.

.tocmtappendix mtchapternone mtappendixsection mtappendixsubsection

##  Appendix S1 Computational details

###  S1.1 Univariate and bivariate smooth function estimation with splines

We discuss how the univariate and bivariate smooth functions are used to model
main effects and interaction effects in Section 2.1. The computational details
for estimating these building blocks via quadratic optimization are outlined
below.

####  S1.1.1 One dimensional nonparametric function estimation.

Suppose that we have observations $\\{({y}_{i},u_{i})\\}_{1}^{n}$
corresponding to a univariate response ${y}_{i}$ and a univariate predictor
$u_{i}\in[0,1]$. If the underlying relationship between $y$ and ${u}$ can be
modeled via a twice continuously differentiable (i.e., smooth) function
$m(u)$, we can estimate $m(u)$ as a function that minimizes the following
objective:
$\sum_{i=1}^{n}(y_{i}-m(u_{i}))^{2}+\lambda\int(m^{\prime\prime}(u))^{2}du$.
This functional optimization problem can be reduced to a finite dimensional
quadratic program by using a suitable basis representation for $m(u)$. For
example, if all of the $u_{i}$’s are distinct and $b_{1}(u),\dots,b_{n}(u)$
are the basis functions [for example, cubic splines with knots at the
observations $u_{i},$ Wahba, 1990, Hastie et al., 2015], then we can write
$m(u)=\sum_{i=1}^{n}\gamma_{i}b_{i}(u)$. In this case,
$\int(m^{\prime\prime}(u))^{2}du=\boldsymbol{\gamma}^{\top}\mathbf{D}\boldsymbol{\gamma}$,
where $\mathbf{D}$ is a positive semidefinite matrix with the $(i,l)$-th entry
given by $\int b^{{}^{\prime\prime}}_{i}(u)b^{{}^{\prime\prime}}_{l}(u)du$.

Instead of using an $n$-dimensional basis representation for $m(u)$, which is
computationally expensive and perhaps statistically redundant, one can choose
a smaller number of basis functions, such as $K=O(n^{1/5})$; this leads to a
low-rank representation of $\mathbf{D}$ [Hall and Opsomer, 2005, Li and
Ruppert, 2008] and hence an improved computational performance. Reducing the
number of basis elements from $n$ to $K$ leads to fewer parameters at the
expense of a marginal increase in bias [Wahba, 1990], which is often
negligible in practice.

The evaluations of the function $m(u)$ at the $n$ data points can be stacked
together in the form of a vector $\mathbf{B}\boldsymbol{\beta}$, where
$\mathbf{B}$ is an $n\times K$ matrix of basis functions evaluated at the
observations, and $\boldsymbol{\beta}\in\mathbb{R}^{K}$ contains the
corresponding basis coefficients that need to be estimated from the data. The
univariate function fitting problem then reduces to a regularized least
squares problem given by
$\min_{\boldsymbol{\beta}\in\mathbb{R}^{K}}\\{\left\lVert\boldsymbol{y}-\boldsymbol{B}\boldsymbol{\beta}\right\rVert_{2}^{2}+\lambda\boldsymbol{\beta}^{T}\boldsymbol{A}\boldsymbol{\beta}\\}.$

While many choices of spline bases are available [Wood, 2006, Hastie et al.,
2001], we use penalized B-splines for their appealing statistical and
computational properties. Following Eilers and Marx [1996], a second-order
finite difference penalty on the coefficients of adjacent B-splines serves as
a good discrete approximation to the integral of the squared second-order
derivative penalty $\Omega$ discussed above. The regularized least squares
problem takes the form

$\displaystyle\min_{\boldsymbol{\beta}}~{}~{}\left\lVert\boldsymbol{y}-\boldsymbol{B}\boldsymbol{\beta}\right\rVert_{2}^{2}+\lambda\sum_{l=1}^{K-2}(\Delta^{2}\beta_{l})^{2},$
(S1)

where $K$ is the number of the B-spline basis functions and
$\Delta^{2}\beta_{l}=\beta_{l}-2\beta_{l+1}+\beta_{l+2}$ is the second-order
finite difference of the coefficients of adjacent B-splines with equidistant
knots. We note that the regularization term can be represented in matrix form
as $\lambda\left\lVert\boldsymbol{D}\boldsymbol{\beta}\right\rVert_{2}^{2}$,
where $\boldsymbol{D}\in\mathbb{R}^{(K-2)\times K}$ is a banded matrix with
nonzero entries given by $d_{l,l}=1$, $d_{l,l+1}=-2$ and $d_{l,l+2}=1$ for
$l\in[K-2]$.

####  S1.1.2 Two dimensional nonparametric function estimation.

To model the response as a function of two covariates one can again use
reduced rank parameterizations, in the form of multivariate splines [Wahba,
1990], thin-plate splines [Kammann and Wand, 2003], or tensor products of
B-splines [Eilers and Marx, 2003, Wood, 2006], for example. We use tensor
products of B-splines to model the two-dimensional smooth functions of the
form $m(u,v)$, with $(u,v)\in[0,1]^{2}$ for concreteness. We start from a low-
rank B-spline basis representation for the marginal smooth functions as
$m_{1}(u)=\sum_{k=1}^{K}\beta_{k}b_{k}(u)$ and
$m_{2}(v)=\sum_{l=1}^{L}\delta_{l}c_{l}(v)$, where $\\{b_{k}\\}$ and
$\\{c_{l}\\}$ are B-spline basis functions, and
$\\{\beta_{k}\\},\\{\delta_{l}\\}$ are unknown basis coefficients. We convert
the marginal smooth function $m_{1}$ into a smooth function of covariates $u$
and $v$ by allowing the coefficients $\beta_{k}$ to vary in a smooth fashion
with respect to $v$. Given that we have an available basis for representing
smooth functions of $v$, we can write
$\beta_{k}(v)=\sum_{l=1}^{L}\delta_{k,l}c_{l}(v)$ and then arrive at the
tensor smooth given by
$m(u,v)=\sum_{k=1}^{K}\sum_{l=1}^{L}\delta_{k,l}b_{k}(u)c_{l}(v)$. We note
that the evaluations of the function $m(u,v)$ at the $n$ data points can be
stacked together in a vector written as $\boldsymbol{R}\boldsymbol{\gamma}$.
We denote the vectors evaluating the two marginal functional bases,
$\\{b_{k}(\cdot)\\}$ and $\\{c_{l}(\cdot)\\}$, at the $i$-th data point as
$\boldsymbol{B}_{i}\in\mathbb{R}^{K}$ and
$\boldsymbol{C}_{i}\in\mathbb{R}^{L}$, respectively. In matrix notation, the
model matrix $\boldsymbol{R}\in\mathbb{R}^{n\times KL}$ can be expressed as
$\boldsymbol{R}=(\boldsymbol{B}\otimes\boldsymbol{1}_{L})\odot(\boldsymbol{1}_{K}\otimes\boldsymbol{C})$,
where operations $\otimes$ and $\odot$ denote Kronecker product and element-
wise multiplication respectively. The basis coefficients $\delta_{k,l}$ are
appropriately ordered into the vector $\boldsymbol{\gamma}\in\mathbb{R}^{KL}$
via the vectorization operation
$\boldsymbol{\gamma}=\text{vec}(\boldsymbol{\delta})=[\delta_{1,1},\cdots,\delta_{K,1},\delta_{1,2},\cdots,\delta_{K,2},\cdots,\delta_{1,L},\cdots,\delta_{K,L}]^{\top}$,
where $\boldsymbol{\delta}=[\delta_{k,l}]_{k\in[K],l\in[L]}$ stores the basis
coefficients in matrix form.

The smooth 2D estimate can be obtained by making use of tensor products and a
discretized version of the smoothness penalty [Eilers and Marx, 2003], given
by:

$\displaystyle\min_{\boldsymbol{\gamma}}\left\lVert\boldsymbol{y}-\boldsymbol{R}\boldsymbol{\gamma}\right\rVert_{2}^{2}+\lambda_{b}\sum_{k=1}^{K}\sum_{l=1}^{L-2}(\Delta_{(1)}^{2}\delta_{k,l})^{2}+\lambda_{c}\sum_{l=1}^{L}\sum_{k=1}^{K-2}(\Delta_{(2)}^{2}\delta_{k,l})^{2},$
(S2)

where
$\Delta_{(1)}^{2}\delta_{k,l}=\delta_{k,l}-2\delta_{k,l+1}+\delta_{k,l+2}$ and
$\Delta_{(2)}^{2}\delta_{k,l}=\delta_{k,l}-2\delta_{k+1,l}+\delta_{k+2,l}$.
The regularization terms above can be compactly represented as quadratic forms
in $\boldsymbol{\gamma}$ as follows:

$\displaystyle\min_{\boldsymbol{\gamma}}\left\lVert\boldsymbol{y}-\boldsymbol{R}\boldsymbol{\gamma}\right\rVert_{2}^{2}+\lambda_{b}\boldsymbol{\gamma}^{T}\boldsymbol{P}_{b}\boldsymbol{\gamma}+\lambda_{c}\boldsymbol{\gamma}^{T}\boldsymbol{P}_{c}\boldsymbol{\gamma},$
(S3)

where
$\boldsymbol{P}_{b}=(\boldsymbol{D}^{T}\boldsymbol{D})\otimes\boldsymbol{I}_{L}$
and
$\boldsymbol{P}_{c}=\boldsymbol{I}_{K}\otimes(\boldsymbol{D}^{T}\boldsymbol{D})$.
The regularizer ensures that the coefficients in the same row (or column) of
$\boldsymbol{\delta}$ vary regularly, leading to a smooth 2D surface.

###  S1.2 Block Coordinate Descent for solving (4)

In our block CD method, the blocks correspond to the basis coefficients for
either the main effects $\\{\boldsymbol{\beta}_{j}\\}_{j}$ or the interaction
effects $\\{\boldsymbol{\theta}_{j,k}\\}_{j,k}$. Given an initialization
$(\boldsymbol{\beta}_{1}^{(0)},\cdots,\boldsymbol{\beta}_{p}^{(0)},\boldsymbol{\theta}_{1,2}^{(0)},\\\
\cdots,\boldsymbol{\theta}_{p-1,p}^{(0)})$, at every cycle, we sequentially
sweep across the main effects and the interaction effects. If we denote the
solution after $t$ cycles by
$(\boldsymbol{\beta}_{1}^{(t)},\cdots,\boldsymbol{\beta}_{p}^{(t)},\boldsymbol{\theta}_{1,2}^{(t)},\cdots,\boldsymbol{\theta}_{p-1,p}^{(t)})$,
then the block of coefficients for $j$-th main effect
$\boldsymbol{\beta}_{j}^{(t+1)}$ at the cycle $t+1$ is obtained by optimizing
(4) with respect to $\boldsymbol{\beta}_{j}$, with other variables held fixed:

$\displaystyle\boldsymbol{\beta}_{j}^{(t+1)}\in\underset{\boldsymbol{\beta}_{j}\in\mathbb{R}^{K_{j}}}{\mathop{\rm
arg\,min}}~{}G(\boldsymbol{\beta}_{1}^{(t+1)},\cdots,\boldsymbol{\beta}_{j-1}^{(t+1)},\boldsymbol{\beta}_{j},\boldsymbol{\beta}_{j+1}^{(t)},\cdots,\boldsymbol{\beta}_{p}^{(t)},~{}~{}\boldsymbol{\theta}_{1,2}^{(t)},\cdots,\boldsymbol{\theta}_{p-1,p}^{(t)}).$
(S4)

Similarly, $\boldsymbol{\theta}_{j,k}^{(t+1)}$, the coefficients for the
$(j,k)$-th interaction effect at cycle $t+1$ are updated as

$\displaystyle\boldsymbol{\theta}_{j,k}^{(t+1)}\in\underset{\boldsymbol{\theta}_{j,k}\in\mathbb{R}^{K_{j,k}}}{\mathop{\rm
arg\,min}}~{}G(\boldsymbol{\beta}_{1}^{(t+1)},\cdots,\boldsymbol{\beta}_{p}^{(t+1)},~{}~{}\boldsymbol{\theta}_{1,2}^{(t+1)},\cdots,\boldsymbol{\theta}_{(j,k)-1}^{(t+1)},\boldsymbol{\theta}_{j,k},\boldsymbol{\theta}_{(j,k)+1}^{(t)}\cdots,\boldsymbol{\theta}_{p-1,p}^{(t)}).$
(S5)

The block minimization problem (S4) reduces to:

$\displaystyle\boldsymbol{\beta}_{j}^{(t+1)}=\underset{\boldsymbol{\beta}_{j}\in\mathbb{R}^{K_{j}}}{\mathop{\rm
arg\,min}}~{}\psi_{j}(\boldsymbol{r}^{(t)};\boldsymbol{\beta}_{j}):=\frac{1}{n}\big{\|}\boldsymbol{r}^{(t)}-\boldsymbol{B}_{j}\boldsymbol{\beta}_{j}\big{\|}_{2}^{2}+\lambda_{1}~{}\boldsymbol{\beta}_{j}^{T}\boldsymbol{S}_{j}\boldsymbol{\beta}_{j}+\lambda_{2}\cdot\mathbb{1}[\boldsymbol{\beta}_{j}\neq\boldsymbol{0}],$
(S6)

where
$\boldsymbol{r}^{(t)}=\boldsymbol{y}-(\sum_{j^{\prime}=1}^{j-1}\boldsymbol{B}_{j^{\prime}}\boldsymbol{\beta}_{j^{\prime}}^{(t+1)}+\sum_{j^{\prime}=j+1}^{p}\boldsymbol{B}_{j^{\prime}}\boldsymbol{\beta}_{j^{\prime}}^{(t)}+\sum_{j^{\prime}<k^{\prime}}\boldsymbol{B}_{j^{\prime},k^{\prime}}\boldsymbol{\theta}_{j^{\prime},k^{\prime}}^{(t)})$
denotes the residual. A solution to (S6) can be computed in closed form via
the following thresholding operator:

$\displaystyle\boldsymbol{\beta}_{j}^{(t+1)}=\begin{cases}\boldsymbol{0}&\text{if}~{}~{}\psi_{j}(\boldsymbol{r}^{(t)};\boldsymbol{0})\leq\min_{{\boldsymbol{\beta}}_{j}\neq\boldsymbol{0}}\psi_{j}(\boldsymbol{r}^{(t)};\boldsymbol{\beta}_{j})\\\
\big{(}\boldsymbol{B}_{j}^{T}\boldsymbol{B}_{j}+n\lambda_{1}\boldsymbol{S}_{j}\big{)}^{-1}\boldsymbol{B}_{j}^{T}\boldsymbol{r}^{(t)}&\text{otherwise}.\end{cases}$

Similarly, the sub-problem for the interaction effects is

$\displaystyle\boldsymbol{\theta}_{j,k}^{(t+1)}=\underset{\boldsymbol{\theta}_{j,k}\in\mathbb{R}^{K_{j,k}}}{\mathop{\rm
arg\,min}}~{}\psi_{j,k}(\boldsymbol{r}^{(t)};\boldsymbol{\theta}_{j,k}):=\frac{1}{n}\big{\|}\boldsymbol{r}^{(t)}-\boldsymbol{B}_{j,k}\boldsymbol{\theta}_{j,k}\big{\|}_{2}^{2}+\lambda_{1}~{}\boldsymbol{\theta}_{j,k}^{T}\boldsymbol{S}_{j,k}\boldsymbol{\theta}_{j,k}+\alpha\lambda_{2}\cdot\mathbb{1}[\boldsymbol{\theta}_{j,k}\neq\boldsymbol{0}],$
(S7)

where
$\boldsymbol{r}^{(t)}=\boldsymbol{y}-(\sum_{j^{\prime}=1}^{p}\boldsymbol{B}_{j^{\prime}}\boldsymbol{\beta}_{j^{\prime}}^{(t+1)}+\sum_{(j^{\prime},k^{\prime})=1,2}^{(j,k)-1}\boldsymbol{B}_{j^{\prime},k^{\prime}}\boldsymbol{\theta}_{j^{\prime},k^{\prime}}^{(t+1)}+\sum_{(j^{\prime},k^{\prime})=(j,k)+1}^{p-1,p}\boldsymbol{B}_{j^{\prime},k^{\prime}}\boldsymbol{\theta}_{j^{\prime},k^{\prime}}^{(t)})$.
A solution to this problem is given by

$\displaystyle\boldsymbol{\theta}_{j,k}^{(t+1)}=\begin{cases}\boldsymbol{0}&\text{if}\;~{}~{}\psi_{j,k}(\boldsymbol{r}^{(t)};\boldsymbol{0})\leq\min_{{\boldsymbol{\theta}}_{j,k}\neq\boldsymbol{0}}\psi_{j,k}(\boldsymbol{r}^{(t)};\boldsymbol{\theta}_{j,k})\\\
\big{(}\boldsymbol{B}_{j,k}^{T}\boldsymbol{B}_{j,k}+N\lambda_{1}\boldsymbol{S}_{j,k}\big{)}^{-1}\boldsymbol{B}_{j,k}^{T}\boldsymbol{r}^{(t)}&\text{otherwise}.\end{cases}$
(S8)

These block CD updates need to be paired with several computational devices in
the form of active set updates, cached matrix factorizations, and warm-starts,
among others. We draw inspiration from similar strategies used in CD-based
procedures for sparse linear regression [Hazimeh and Mazumder, 2020a, Friedman
et al., 2010], and adapt them to our problem. These devices are discussed in
detail in the Supplement Section S1.3.

###  S1.3 Scalability considerations for solving (4)

We draw inspiration from strategies used in CD-based procedures for sparse
linear regression [Hazimeh and Mazumder, 2020a, Friedman et al., 2010], and
adapt them to scale cyclic block coordinate descent in Section 2.4 to large
problem instances of (4) . These include active set updates, cached matrix
factorizations and warm-starts. They are explained in more detail below:

####  S1.3.1 Active set updates

A main computational bottleneck in the block CD approach is the number of
passes across the $O(p^{2})$ blocks. However, as we anticipate a solution that
is sparse (with few nonzero main and interaction effects), we use an active
set strategy. We restrict our block CD procedure to a small subset
$\mathcal{Q}$ of the $O(p^{2})$ variables, with all blocks outside the active
set being set to zero. Once the CD algorithm converges on the active set, we
check if all blocks outside the active set satisfy the coordinate-wise
optimality conditions111That is, we check if optimal solutions to (10) and
(11) are zero for all blocks outside ${\mathcal{Q}}$.. If there are any
violations, we select the corresponding blocks, append them to the current
active set, and then rerun our block CD procedure. As there are finitely many
active sets, the algorithm is guaranteed to converge – in practice, with warm-
start continuation (discussed below), the number of active set updates is
quite small, and the algorithm is found to converge quite quickly.

####  S1.3.2 Cached matrix factorizations

The updates (10) and (11) require computing a matrix inverse. In particular,
if a block is nonzero, we need to compute a linear system solution of the
form: $\boldsymbol{A}_{l}^{-1}\boldsymbol{b}_{l}$ where
$\boldsymbol{A}_{l}=\boldsymbol{B}_{l}^{T}\boldsymbol{B}_{l}+n\lambda_{1}\boldsymbol{S}_{l}$.
We note that matrix $\boldsymbol{A}_{l}$ is fixed throughout the CD updates
and is independent of the choice of the sparsity regularization parameter
$\lambda_{2}$ – hence, we pre-compute a matrix factorization for
$\boldsymbol{A}_{l}$ (for example, an $LU$ decomposition) and use it to
compute the solution to the linear system. As the dimension of
$\boldsymbol{A}_{l}$ equals the number of basis coefficients, which is small,
this can be done quite efficiently.

####  S1.3.3 Warm-starts

We use our CD procedure to compute a path of solutions to (4) for a 2D grid of
tuning parameters
$(\lambda_{1},\lambda_{2})\in\\{\lambda_{1}^{(l)}\\}_{l=0}^{L}\times\\{\lambda_{2}^{(m)}\\}_{m=0}^{M}$,
where $\lambda_{1}$ corresponds to the smoothness parameter and $\lambda_{2}$
the sparsity parameter. Here, $\lambda_{1}^{(l)}>\lambda_{1}^{(l+1)}$ for all
$l$, and $\lambda_{2}^{(m)}>\lambda_{2}^{(m+1)}$ for all $m$. When
$\lambda_{1}=\lambda_{1}^{(0)}$ (most regularized), we compute a sequence of
solutions across the $\lambda_{2}$-values (from large to small values): a
solution obtained at $(\lambda_{1}^{(0)},\lambda_{2}^{(m)})$ is used to
initialize our CD procedure for the value
$(\lambda_{1}^{(0)},\lambda_{2}^{(m+1)})$. As the number of nonzeros in a
solution to (4) generally increases with $\lambda_{2}$-values, our CD
procedure for $(\lambda_{1}^{(0)},\lambda_{2}^{(m+1)})$ uses an active set
that is slightly larger than the active set222This is usually taken to be
1-10% larger than the current active set, and is chosen in a greedy fashion
from among the main and interaction effects lying outside the current active
set. corresponding to the solution at $(\lambda_{1}^{(0)},\lambda_{2}^{(m)})$.
Once we have traced a full path over $\lambda_{2}$, we use warm-starts in the
lateral direction across the space of $\lambda_{1}$. For all $l\geq 0$, to
obtain a solution to (4) at $(\lambda_{1}^{(l+1)},\lambda_{2}^{(m)})$, we use
the solution at $(\lambda_{1}^{(l)},\lambda_{2}^{(m)})$ as a warm-start.

###  S1.4 Algorithms for sparse hierarchical interactions: Problem (5)

Similar to the case of (4), and with the computational speed in mind, we
present approximate methods to obtain good solutions to (6). As mentioned
earlier, these approximate solutions can be used to initialize MIP-based
approaches for (6) to improve the solution, and/or to certify the quality of
these approximate solutions following Bertsimas et al. [2016], Hazimeh and
Mazumder [2020a].

Our first step is to reduce the number of main and interaction effects in (5)
by making use of the family of solutions available from (4). We consider the
union of supports available from the family of solutions obtained from (4),
across the 2D grid of tuning parameters $\lambda_{1},\lambda_{2}$. Let
$\mathcal{M}\subset[p]$ and $\mathcal{I}\subset[p(p-1)/2]$ denote the sets of
all nonzero main effects and interactions effects, respectively333If
$\mathcal{I}$ includes an interaction $(j,k)$ where main-effect $j$ is not
included in $\mathcal{M}$, we expand $\mathcal{M}$ to include $j$. This way,
we make sure that all interaction effects have the corresponding main effects
included in $\mathcal{M}$., encountered along the 2D path of solutions to (4).
We form a reduced version of problem (6) with $z_{i}=0,i\notin{\mathcal{M}}$
and $z_{j,k}=0$ for all $(j,k)\notin{\mathcal{I}}$. Let us denote the reduced
problem by ${\mathcal{P}}(\mathcal{M},\mathcal{I})$. We consider a convex
relaxation of ${\mathcal{P}}(\mathcal{M},\mathcal{I})$, denoted by
${\mathcal{P}}^{R}(\mathcal{M},\mathcal{I})$, where all binary variables
$\\{z_{i}\\}$, $\\{z_{j,k}\\}$ in ${\mathcal{P}}(\mathcal{M},\mathcal{I})$ are
relaxed to the interval $[0,1]$. As the sizes of $\mathcal{M}$ and
$\mathcal{I}$ are generally small, it is computationally feasible to solve the
relaxation ${\mathcal{P}}^{R}(\mathcal{M},\mathcal{I})$ – we let
$\\{z^{R}_{i}\\}$ and $\\{z^{R}_{j,k}\\}$ denote a solution to this
relaxation. Following Hazimeh and Mazumder [2020b], it can be shown that this
solution satisfies the strong hierarchy constraint (almost surely). To obtain
a feasible solution to Problem (6), we apply a relax-and-round procedure. For
the rounding step, we consider a threshold $\tau\in(0,1)$ and obtain
$\tilde{z}_{i}=\mathbb{1}[z^{R}_{i}>\tau]$ for all $i\in\mathcal{M}$ and
$\tilde{z}_{j,k}=\mathbb{1}[z^{R}_{j,k}>\tau]$ for all $(j,k)\in\mathcal{I}$.
We set $\tilde{z}_{i}=0$, $i\notin{\mathcal{M}}$; and $\tilde{z}_{j,k}=0$ for
all $(j,k)\notin{\mathcal{I}}$. It can be verified that this rounding
procedure maintains strong hierarchy. Finally, we perform a ‘polishing’ step
where we solve (3) restricted to the support defined by
$\\{\tilde{z}_{i}\\}_{i}$ and $\\{\tilde{z}_{j,k}\\}_{j,k}$.

Related Work. In contrast to Problem (4), the regularization penalty in
problem (6) is not separable across the blocks due to the overlapping groups
created by the strong hierarchy constraint. Hence, the CD-based procedures
discussed for (4) do not apply to the hierarchical setting. To the best of our
knowledge, there are no prior specialized algorithms for (5) that apply to the
scale that we consider here. In fact, even in the linear regression setting,
current algorithms for problems with a hierarchy constraint are somewhat
limited in terms of the problem-scales they can address. The sole exceptions
appear to be the convex optimization based approaches of Lim and Hastie
[2015], Hazimeh and Mazumder [2020b], which can address sparse linear
regression problems with a large number of features and a small number of
observations.

##  Appendix S2 Statistics Theory: technical details and proofs

###  S2.1 Notation and assumptions for Theorem 1

Given a candidate regression function $f:[0,1]^{p}\mapsto\mathbb{R}$ of the
form
$f(\mathbf{x})=\sum_{j\in[p]}f_{j}(x_{j})+\sum_{j<k}f_{j,k}(x_{j},x_{k})$, we
define its sparsity level as

$G(f)=\sum_{j\in[p]}\mathbb{1}[{\bf
f}_{j}\neq\mathbf{0}]+\sum_{j<k}\mathbb{1}[{\bf f}_{j,k}\neq\mathbf{0}].$

To ensure identifiability of the additive representation for $f(\mathbf{x})$
additional restrictions need to be imposed on the components $f_{j}$ and
$f_{j,k}$. For the simplicity of the presentation, we avoid specifying a
particular set of restrictions and treat every representation of $f$ as
equivalent, with the understanding that one particular representation is used
when evaluating quantities such as $G(f)$.

We write $\|\cdot\|_{L_{2}}$ for the $L_{2}$ norm of a real-valued function on
$[0,1]^{2}$ or $[0,1]$. We focus on the case where $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$ are $L_{2}$-Sobolev spaces. More specifically, we let

$\mathcal{C}_{2}=\left\\{g(u,v):[0,1]^{2}\mapsto\mathbb{R},\;\|g\|_{L_{2}}+\left\|\frac{\partial^{2}g}{\partial
u^{2}}\right\|_{L_{2}}+\left\|\frac{\partial^{2}g}{\partial u\partial
v}\right\|_{L_{2}}+\left\|\frac{\partial^{2}g}{\partial
v^{2}}\right\|_{L_{2}}<\infty\right\\}$
$\Omega(g)=\left\|\frac{\partial^{2}g}{\partial
u^{2}}\right\|_{L_{2}}^{2}+\left\|\frac{\partial^{2}g}{\partial u\partial
v}\right\|_{L_{2}}^{2}+\left\|\frac{\partial^{2}g}{\partial
v^{2}}\right\|^{2}_{L_{2}},\vskip 3.0pt plus 1.0pt minus 1.0pt$

and also let
${\mathcal{C}_{1}}=\left\\{h:[0,1]\mapsto\mathbb{R},\;\|h\|_{L_{2}}+\|h^{\prime\prime}\|_{L_{2}}<\infty\right\\}$,
${\Omega}(h)=\|h^{\prime\prime}\|^{2}_{L_{2}}$. We define the corresponding
space of additive functions with interactions as

$\mathcal{C}_{\text{gr}}=\\{f:[0,1]^{p}\mapsto\mathbb{R},\;f(\mathbf{x})=\sum_{j\in[p]}f_{j}(x_{j})+\sum_{j<k}f_{j,k}(x_{j},x_{k}),\;f_{j}\in\mathcal{C}_{1},\;f_{j,k}\in\mathcal{C}_{2},\;G(f)\leq
K\\},$

where $K$ is some arbitrarily large universal constant, and let

$\Omega_{\text{gr}}(f)=\sum_{j\in[p]}{\Omega}(f_{j})+\sum_{j<k}\Omega(f_{j,k}).$

We associate each $f\in\mathcal{C}_{\text{gr}}$ with the vector ${\bf
f}=\sum_{j\in[p]}{\bf f}_{j}+\sum_{j<k}{\bf f}_{j,k}$, where ${\bf
f}_{j}=\big{(}f_{j}(x_{1j}),...,\\\ f_{j}(x_{nj})\big{)}$ and ${\bf
f}_{j,k}=\big{(}f_{j,k}(x_{1j},x_{1k}),...,f_{j,k}(x_{nj},x_{nk})\big{)}$.

We assume that the observed data follows the model $\mathbf{y}={\bf
f}^{*}+\boldsymbol{\epsilon}$, where ${\bf
f}^{*}=\big{(}f(\mathbf{x}_{1}),...,\\\ f(\mathbf{x}_{n})\big{)}$ is the
vector representation of a function $f^{*}:[0,1]^{p}\mapsto\mathbb{R}$, and
the elements of $\boldsymbol{\epsilon}$ are independent $N(0,\sigma^{2})$ with
$\sigma>0$.

###  S2.2 Preliminaries and supplementary results

We will prove a more general result by replacing the $2$-nd derivative in the
definition of $\mathcal{C}_{\text{gr}}$ and $\Omega_{\text{gr}}$ with an
$m$-th derivative, and redefining $r_{n}$ as $n^{-m/(2m+2)}$. Theorem 1 will
then follow as a special case corresponding to $m=2$. The more general
prediction error bound extends the result in Theorem 2 of Lin and Zhang
[2006], which is established for the fixed $p$ setting, from the main effects
models to the interaction models. When $\log(p)\lesssim n^{1/(m+1)}$, the
prediction error rate matches the optimal bivariate rate of $n^{-m/(m+1)}$.

We extend the domain of $\|\cdot\|_{n}$ from vectors in $\mathbb{R}^{n}$ to
real-valued functions on $[0,1]^{p}$ by letting $\|\cdot\|_{n}$ be the
empirical $L_{2}$-norm. Thus, given a function $h$, we let
$\|h\|_{n}=[\sum_{i=1}^{n}h(\mathbf{x}_{i})^{2}/n]^{1/2}$. This extension is
consistent in the sense that $\|f\|_{n}=\|{\bf f}\|_{n}$ for
$f\in\mathcal{C}_{\text{\rm gr}}$. By analogy with the $\|\cdot\|_{n}$
notation, we define $(\boldsymbol{\epsilon},{\bf
v})_{n}=(1/n)\sum_{i=1}^{n}\epsilon_{i}v_{i}$ for each ${\bf
v}\in\mathbb{R}^{n}$. For the simplicity of the presentation, we use the
double-index notation $f_{jj}$ for the main effect components $f_{j}$ of
functions $f\in\mathcal{C}_{\text{gr}}$ and then take advantage of the
additive representation

$\sum_{j\in[p]}f_{j}(x_{j})+\sum_{j<k}f_{j,k}(x_{j},x_{k})=\sum_{j\leq
k}f_{j,k}(x_{j},x_{k}).$

We also define $\tilde{\Omega}(f_{j,k})=\sqrt{\Omega(f_{j,k})}$ and
$\tilde{\Omega}_{\text{gr}}(f)=\sum_{j\leq k}\sqrt{\Omega(f_{j,k})}$.

We let $\mathcal{J}=\\{(j,k)\in[p]\times[p],\;\text{s.t.}\;j\leq k\\}$. Given
$J\subseteq\mathcal{J}$, we define functional classes
$\mathcal{F}(J)=\\{f:\;f(\mathbf{x})=\sum_{(j,k)\in
J}f_{j,k}(x_{j},x_{k}),\;f_{j,k}\in\mathcal{C}_{2}\\}$ and
$\mathcal{H}({J})=\\{h:\;h\in\mathcal{F}(J),\;\|h\|_{n}/(r_{n}+\tilde{\Omega}_{\text{\rm
gr}}(h)\leq 1\\}$; we also let
$\mathcal{F}_{s}=\\{f:\;f\in\mathcal{C}_{\text{\rm gr}},\;G(f)\leq s\\}$ for
every natural $s$. Given a positive constant $\delta$ and a metric space
$\mathcal{H}$ endowed with the norm $\|\cdot\|$, we use the standard notation
and write $H(\delta,\mathcal{H},\|\cdot\|)$ for the $\delta$-entropy of
$\mathcal{H}$ with respect to $\|\cdot\|$. More specifically,
$H(\delta,\mathcal{H},\|\cdot\|)$ is the natural logarithm of the smallest
number of balls with radius $\delta$ needed to cover $\mathcal{H}$. The
following result, proved in Section S2.5, bounds the entropy of
$\mathcal{H}({J})$.

######  Lemma 1.

$H(u,\mathcal{H}({J}),\|\cdot\|_{n})\lesssim(1/u)^{2/m}$ for $u\in(0,1)$.

Let $M=K(K+1)/2$. We will need the following maximal inequalities, which are
proved in Sections S2.6 and S2.7, respectively.

######  Lemma 2.

Suppose that $J\subseteq\mathcal{J}$ and $|J|\leq M$. Then, with probability
at least $1-e^{-t}$, inequality

$(\boldsymbol{\epsilon}/\sigma,{\bf
f})_{n}\lesssim\Big{[}r_{n}+\sqrt{t/n}\Big{]}\|{\bf
f}\|_{n}+\Big{[}r_{n}^{2}+r_{n}\sqrt{t/n}\Big{]}\tilde{\Omega}_{\text{\rm
gr}}(f)$

holds uniformly over ${f\in\mathcal{F}(J)}$.

######  Lemma 3.

With probability at least $1-\epsilon$, inequality

$\displaystyle(\boldsymbol{\epsilon}/\sigma,{\bf f})_{n}$
$\displaystyle\lesssim$
$\displaystyle\Big{[}r_{n}+\sqrt{\frac{\log(ep)}{n}}+\sqrt{\frac{\log(1/\epsilon)}{n}}\Big{]}\|{\bf
f}\|_{n}$
$\displaystyle+\Big{[}r_{n}^{2}+r_{n}\sqrt{\frac{\log(ep)}{n}}+r_{n}\sqrt{\frac{\log(1/\epsilon)}{n}}\Big{]}\tilde{\Omega}_{\text{\rm
gr}}(f)$

holds uniformly over ${f\in\mathcal{F}_{M}}$.

###  S2.3 Proof of Theorem 1

Consider an arbitrary function $f\in\mathcal{C}_{\text{gr}}$. For the
remainder of the proof, all universal constants will be chosen independently
of $f$. Because $f$ is feasible for optimization problem (7), we have the
following inequality:

$\|\widehat{\bf f}-{\bf
f}^{*}\|_{n}^{2}+\lambda_{n}\Omega_{\text{gr}}(\widehat{f})+\mu_{n}G(\widehat{f})\leq\|{\bf
f}-{\bf f}^{*}\|_{n}^{2}+2(\boldsymbol{\epsilon},\widehat{\bf f}-{\bf
f})_{n}+\lambda_{n}\Omega_{\text{gr}}(f)+\mu_{n}G(f).$ (S9)

Applying Lemma 3 with $f=\widehat{f}-f$ and $\epsilon=(ep)^{-1}$, we derive
that, with probability at least $1-1/p$,

$\displaystyle(\boldsymbol{\epsilon}/\sigma,\widehat{\bf f}-{\bf f})_{n}$
$\displaystyle\leq$
$\displaystyle{c}_{1}\Big{[}r_{n}+\sqrt{\frac{\log(ep)}{n}}\Big{]}\|\widehat{\bf
f}-{\bf f}\|_{n}$ (S10)
$\displaystyle+a_{2}\Big{[}r_{n}^{2}+r_{n}\sqrt{\frac{\log(ep)}{n}}\Big{]}\tilde{\Omega}_{\text{\rm
gr}}(\widehat{f}-f)$

for some universal constants $a_{1}$ and $a_{2}$. For the remainder of the
proof we restrict our attention to the random event on which (S10) holds.
Multiplying inequality (S9) by two, using (S10), and letting

$\tau_{n}:=2{c}_{1}\sigma\Big{[}r_{n}+\sqrt{\frac{\log(ep)}{n}}\Big{]},\;\tilde{\lambda}_{n}:=4a_{2}\sigma\Big{[}r_{n}^{2}+r_{n}\sqrt{\frac{\log(ep)}{n}}\Big{]},\;\text{and}\;{\lambda}_{n}\geq\tilde{\lambda}_{n},$

we derive

$2\|\widehat{\bf f}-{\bf
f}^{*}\|_{n}^{2}+2\lambda_{n}\Omega_{\text{gr}}(\widehat{f})+2\mu_{n}G(\widehat{f})\leq
2\|{\bf f}-{\bf f}^{*}\|_{n}^{2}+2\tau_{n}\|\widehat{\bf f}-{\bf
f}\|_{n}+\tilde{\lambda}_{n}\tilde{\Omega}_{\text{gr}}(\widehat{f}-f)+2\lambda_{n}\Omega_{\text{gr}}(f)+2\mu_{n}G(f).$

Applying inequalities $2\tau_{n}\|\widehat{\bf f}-{\bf
f}^{*}\|_{n}\leq\tau_{n}^{2}+\|\widehat{\bf f}-{\bf f}^{*}\|_{n}^{2}$ and
$2\tau_{n}\|{\bf f}^{*}-{\bf f}\|_{n}\leq\tau_{n}^{2}+\|{\bf f}-{\bf
f}^{*}\|_{n}^{2}$, we arrive at

$\displaystyle 2\|\widehat{\bf f}-{\bf
f}^{*}\|_{n}^{2}+2\lambda_{n}\Omega_{\text{gr}}(\widehat{f})+2\mu_{n}G(\widehat{f})$
$\displaystyle\leq$ $\displaystyle\|\widehat{\bf f}-{\bf
f}^{*}\|_{n}^{2}+3\|{\bf f}-{\bf f}^{*}\|_{n}^{2}+2\tau_{n}^{2}$
$\displaystyle+\tilde{\lambda}_{n}\big{[}\tilde{\Omega}_{\text{gr}}(\widehat{f})+\tilde{\Omega}_{\text{gr}}(f)\big{]}+2\lambda_{n}\Omega_{\text{gr}}(f)+2\mu_{n}G(f).$

Taking into account inequality $\tilde{\Omega}_{\text{gr}}(\tilde{f})\lesssim
G(\tilde{f})+\Omega_{\text{gr}}(\tilde{f})$, which holds for any
$\tilde{f}\in\mathcal{C}_{\text{gr}}$, we conclude that

$\|\widehat{\bf f}-{\bf
f}^{*}\|_{n}^{2}+\lambda_{n}\Omega_{\text{gr}}(\widehat{f})+2\mu_{n}G(\widehat{f})\leq
3\|{\bf f}-{\bf
f}^{*}\|_{n}^{2}+a_{3}\tau_{n}^{2}+2\lambda_{n}\Omega_{\text{gr}}(f)+2\mu_{n}G(f)$
(S11)

for some universal constant $a_{3}$. ∎

###  S2.4 Proof of Corollary 1

The stated prediction error bound is a direct consequence of the bound in
Theorem 1, hence we only need to establish the sparsity bound. Taking
$f=f^{*}$, we rewrite inequality (S11) as follows:

$\|\widehat{\bf f}-{\bf
f}^{*}\|_{n}^{2}+\lambda_{n}\Omega_{\text{gr}}(\widehat{f})+2\mu_{n}G(\widehat{f})\leq
a_{3}\tau_{n}^{2}+2\lambda_{n}\Omega_{\text{gr}}(f^{*})+2\mu_{n}G(f^{*}).$
(S12)

Provided that constant $c_{2}$ in the lower bound for $\mu_{n}$ is
sufficiently large, we then have $G(\widehat{f})-G(f^{*})<1$, which implies
$G(\widehat{f})\leq G(f^{*})$.

We now focus on the classical asymptotic setting, where $p$ is fixed and $n$
tends to infinity, under the conditions of Corollary 1. Noting that
$1/n=o(r_{n}^{2})$ and repeating the arguments in the proofs of Theorem 1 and
Corollary 1 with $\epsilon=\exp(-nr_{n}^{2})$, we conclude that inequality
(S12) and the accompanying sparsity bound $G(\widehat{f})\leq G(f^{*})$ hold
with probability tending to one. Imposing an additional requirement that
$\lambda_{n}=O\big{(}r_{n}^{2}+\log(ep)/n\big{)}$ and
$\mu_{n}=O\big{(}r_{n}^{2}+\log(ep)/n\big{)}$ then leads to $\|\widehat{\bf
f}-{\bf f}^{*}\|_{n}^{2}=O\big{(}r_{n}^{2}+\log(ep)/n\big{)}$. ∎

###  S2.5 Proof of Lemma 1

We will establish the stated entropy bound for a somewhat larger functional
space
${\mathcal{H}}_{J}^{\prime}=\\{h:\;h\in\mathcal{F}(J),\;\|h\|_{n}+\tilde{\Omega}_{\text{\rm
gr}}(h)\leq 1\\}$. We treat $m$ as a fixed universal constant in all
inequalities that follow.

Consider an arbitrary $g\in\mathcal{C}_{2}$. By the Sobolev embedding theorem
[for example, Oden and Reddy, 1976, Theorem 3.13], we can write $g$ as a sum
of a polynomial of degree $m-1$ and a function $\tilde{g}$ that satisfies
$\|\tilde{g}\|_{L_{2}}\lesssim\tilde{\Omega}(g)$, where we note that
$\tilde{\Omega}(g)=\tilde{\Omega}(\tilde{g})$. Applying Lemma 10.9 in Van de
Geer [2000], which builds on the interpolation inequality of Agmon [1965], we
derive $\|\tilde{g}\|_{\infty}\lesssim\tilde{\Omega}(\tilde{g})$. Thus,
${\mathcal{H}}_{J}^{\prime}\subseteq\\{p+\tilde{h}:\;p\in\mathcal{P}_{J},\;\tilde{h}\in\tilde{\mathcal{H}}_{J}\\}$,
where

$\displaystyle\mathcal{P}_{J}$ $\displaystyle=$
$\displaystyle\\{p:\;p(\mathbf{x})=\sum_{(j,k)\in
J}\sum_{l=0}^{m-1}\sum_{q=0}^{m-1}\alpha_{j,k,l,q}x_{j}^{l}x_{k}^{q},\;\alpha_{j,k,l,q}\in\mathbb{R},\;\|p\|_{n}\leq
2\\}$ $\displaystyle\tilde{\mathcal{H}}_{J}$ $\displaystyle=$
$\displaystyle\\{\tilde{h}:\;\tilde{h}\in\mathcal{F}(J),\;\tilde{\Omega}_{\text{\rm
gr}}(\tilde{h})\leq
1,\;\|\tilde{h}_{j,k}\|_{\infty}\lesssim\tilde{\Omega}(\tilde{h}_{j,k})\;\forall(j,k)\in
J\\}.$

Bound $\|p\|_{n}\leq 2$ in the definition of $\mathcal{P}_{J}$ holds because
if $h=p+\tilde{h}$ for $h\in{\mathcal{H}}_{J}^{\prime}$ and
$\tilde{h}\in\tilde{\mathcal{H}}_{J}$, then $\|p+\tilde{h}\|_{n}\leq 1$ and
$\|\tilde{h}\|_{n}\leq\tilde{\Omega}_{\text{\rm gr}}(\tilde{h})\leq 1$.
Consequently,

$H(u,\mathcal{H}({J}),\|\cdot\|_{n})\leq
H(u,\mathcal{H}_{J}^{\prime},\|\cdot\|_{n})\leq
H(u/2,\mathcal{P}_{J},\|\cdot\|_{n})+H(u/2,\tilde{\mathcal{H}}_{J},\|\cdot\|_{\infty}),$
(S13)

where we used the fact that the unit ball with respect to the
$\|\cdot\|_{\infty}$-norm is contained within the corresponding ball with
respect to the $\|\cdot\|_{n}$-norm. We note that $\mathcal{P}_{J}$ is a ball
of radis $2$, with respect to the $\|\cdot\|_{n}$-norm, in a linear functional
space of dimension $\lesssim|J|+1$. Hence,
$H(u/2,\mathcal{P}_{J},\|\cdot\|_{n})\lesssim|J|+|J|\log(1/u)$ by, for
example, Corollary 2.6 in Van de Geer [2000]. Thus, the result of Lemma 1
follows from S13 if we also establish that
$H(\delta,\tilde{\mathcal{H}}_{J},\|\cdot\|_{\infty})\lesssim|J|(1/\delta)^{2/m}$
for $\delta\in(0,1)$.

We complete the proof by deriving the stated bound on
$H(\delta,\tilde{\mathcal{H}}_{J},\|\cdot\|_{\infty})$. We represent
$\tilde{\mathcal{H}}_{J}$ as

$\Big{\\{}\tilde{h}:\;\tilde{h}(\mathbf{x})=\sum_{{(j,k)}\in
J}\lambda_{(j,k)}g_{(j,k)}(x_{(j,k)}),\,\sum_{{(j,k)}\in
J}|\lambda_{(j,k)}|\leq
1,\,g_{(j,k)}\in\mathcal{C}_{2},\,\tilde{\Omega}(g_{(j,k)})\leq
1,\,\|g_{(j,k)}\|_{\infty}\leq 1\Big{\\}}.$

Given functions $\tilde{h}(\mathbf{x})=\sum_{{(j,k)}\in
J}\lambda_{(j,k)}g_{(j,k)}(x_{j},x_{k}))$ and
$\tilde{h}^{\prime}(\mathbf{x})=\sum_{{(j,k)}\in
J}\lambda_{(j,k)}^{\prime}g_{(j,k)}^{\prime}(x_{j},x_{k}))$ in
$\tilde{\mathcal{H}}_{J}$, we have

$\displaystyle\|\tilde{h}-\tilde{h}^{\prime}\|_{\infty}$ $\displaystyle\leq$
$\displaystyle\|\sum_{{(j,k)}\in J}\lambda_{(j,k)}g_{(j,k)}-\sum_{{(j,k)}\in
J}\lambda_{(j,k)}g_{(j,k)}^{\prime}\|_{\infty}+\|\sum_{{(j,k)}\in
J}\lambda_{(j,k)}g_{(j,k)}^{\prime}-\sum_{{(j,k)}\in
J}\lambda_{(j,k)}^{\prime}g_{(j,k)}^{\prime}\|_{\infty}$ $\displaystyle\leq$
$\displaystyle\max_{{(j,k)}\in
J}\|g_{(j,k)}-g_{(j,k)}^{\prime}\|_{\infty}\sum_{{(j,k)}\in
J}|\lambda_{(j,k)}|+\max_{{(j,k)}\in
J}\|g_{(j,k)}^{\prime}\|_{\infty}\sum_{{(j,k)}\in
J}|\lambda_{(j,k)}-\lambda_{(j,k)}^{\prime}|$ $\displaystyle\leq$
$\displaystyle\max_{{(j,k)}\in
J}\|g_{(j,k)}-g_{(j,k)}^{\prime}\|_{\infty}+\sum_{{(j,k)}\in
J}|\lambda_{(j,k)}-\lambda_{(j,k)}^{\prime}|.$

Consequently, if we let
$\mathcal{G}=\\{g:\;g\in\mathcal{C}_{2},\;\tilde{\Omega}(g)\leq
1,\;\|g\|_{\infty}\leq 1\\}$, let $\|\cdot\|_{1}$ denote the $\ell_{1}$-norm,
and let $B^{d}_{1}$ denote a unit $\ell_{1}$-ball in $\mathbb{R}^{d}$, then

$H(\delta,\tilde{\mathcal{H}}_{J},\|\cdot\|_{\infty})\leq|J|H(\delta/2,\mathcal{G},\|\cdot\|_{\infty})+H(\delta/2,B^{|J|}_{1},\|\cdot\|_{1}).$

By the results in Birman and Solomjak [1967],
$H(\delta/2,\mathcal{G},\|\cdot\|_{\infty})\lesssim(1/\delta)^{2/m}$. By the
standard bounds on the covering numbers of a norm ball,
$H(\delta/2,B^{|J|}_{1},\|\cdot\|_{1})\lesssim|J|+|J|\log(1/\delta)$.
Consequently,
$H(\delta,\tilde{\mathcal{H}}_{J},\|\cdot\|_{\infty})\lesssim|J|(1/\delta)^{2/m}$
for $\delta\in(0,1)$. ∎

###  S2.6 Proof of Lemma 2

We note that $\|h\|_{n}\leq r_{n}$ and $\tilde{\Omega}_{\text{\rm gr}}(h)\leq
1$ for every $h\in\mathcal{H}({J})$. By Lemma 12 in the supplementary material
for Tan and Zhang [2019] (cf. Corollary 8.3 in Van de Geer, 2000),

$\sup_{h\in\mathcal{H}({J})}(\boldsymbol{\epsilon}/\sigma,{\bf h})_{n}\lesssim
n^{-1/2}\int_{0}^{r_{n}}\sqrt{H(u,\mathcal{H}({J}),\|\cdot\|_{n})}du+r_{n}\sqrt{t/n}$
(S14)

with probability at least $1-e^{-t}$.

We note that $r_{n}=n^{-m/(2m+2)}$ and, thus,
$n^{-1/2}r_{n}^{(m-1)/m}=r_{n}^{2}$. Using Lemma 1 to bound the entropy, we
derive

$\displaystyle
n^{-1/2}\int_{0}^{r_{n}}\sqrt{H(u,\mathcal{H}({J}),\|\cdot\|_{n})}du$
$\displaystyle\lesssim$ $\displaystyle
n^{-1/2}\int_{0}^{r_{n}}u^{-1/m}du\lesssim n^{-1/2}r_{n}^{(m-1)/m}=r_{n}^{2}.$

Applying bound (S14), we conclude that

$\sup_{h\in\mathcal{H}({J})}(\boldsymbol{\epsilon}/\sigma,{\bf h})_{n}\lesssim
r_{n}^{2}+r_{n}\sqrt{t/n}$

with probability at least $1-e^{-t}$. The statement of the lemma is then a
consequence of the fact that for every $f\in\mathcal{F}({J})$, function
$f/\big{[}\|f\|_{n}r_{n}^{-1}+\tilde{\Omega}_{\text{\rm gr}}(f)\big{]}$ falls
in the class $\mathcal{H}({J})$. ∎

###  S2.7 Proof of Lemma 3

Let $M_{s}$ denote the number of distinct subsets of $\mathcal{J}$ that have
size $M$ or smaller. We note that $\log(M_{s})\leq 4M\log(ep)$ and, thus,
$M_{s}e^{-t}\leq e^{4M\log(ep)-t}$. Applying Lemma 2 together with the union
bound, we derive that, with probability at least $1-e^{4M\log(ep)-t}$,
inequality

$(\boldsymbol{\epsilon}/\sigma,{\bf
f})_{n}\lesssim\Big{[}r_{n}+\sqrt{t/n}\Big{]}\|{\bf
f}\|_{n}+\Big{[}r_{n}^{2}+r_{n}\sqrt{t/n}\Big{]}\tilde{\Omega}_{\text{\rm
gr}}(f)$

holds uniformly over ${f\in\mathcal{F}_{M}}$. We complete the proof by noting
that for $t=4M\log(ep)+\log(1/\epsilon)$ the above inequality becomes

$\displaystyle(\boldsymbol{\epsilon}/\sigma,{\bf f})_{n}$
$\displaystyle\lesssim$
$\displaystyle\Big{[}r_{n}+\sqrt{\frac{\log(ep)}{n}}+\sqrt{\frac{\log(1/\epsilon)}{n}}\Big{]}\|{\bf
f}\|_{n}$
$\displaystyle+\Big{[}r_{n}^{2}+r_{n}\sqrt{\frac{\log(ep)}{n}}+r_{n}\sqrt{\frac{\log(1/\epsilon)}{n}}\Big{]}\tilde{\Omega}_{\text{\rm
gr}}(f),$

and the corresponding lower-bound on the probability simplifies to
$1-\epsilon$. ∎

##  Appendix S3 Data Analysis: Additional Details

###  S3.1 Definition of important Census/American Community Survey variables

####  S3.1.1 Definition of the variables in Figure 3

  * •

Tot_Population_ACS_13_17: U.S.resident population includes everyone who meets
the ACS residence rules in the tract at the time of the ACS interview.

  * •

pct_Prs_Blw_Pov_Lev_ACS_13_17: Percentage of the ACS eligible population that
are classified as below the poverty level given their total family or
household income within the last year.

  * •

pct_College_ACS_13_17: The percentage of the ACS population aged 25 years and
over that have a college degree or higher

  * •

pct_Not_HS_Grad_ACS_13_17: The percentage of the ACS population aged 25 years
and over that are not high school graduates and have not received a diploma or
the equivalent.

  * •

pct_Pop_5_17_ACS_13_17: The percentage of the ACS population that is between 5
and 17 years old.

  * •

pct_Pop_18_24_ACS_13_17: The percentage of the ACS population that is between
18 and 24 years old

  * •

pct_Pop_25_44_ACS_13_17: The percentage of the ACS population that is between
25 and 44 years old

  * •

pct_Pop_45_64_ACS_13_17: The percentage of the ACS population that is between
45 and 64 years old

  * •

pct_Pop_65plus_ACS_13_17: The percentage of the ACS population that is 65years
old or over

  * •

pct_Hispanic_ACS_13_17: The percentage of the ACS population that identify as
“Mexican”, “Puerto Rican”, “Cuban”, or “another Hispanic, Latino, or Spanish
origin”.

  * •

pct_NH_White_alone_ACS_13_17: The percentage of the ACS population that
indicate no Hispanic origin and their only race as “White” or report entries
such as Irish, German, Italian, Lebanese, Arab, Moroccan, or Caucasian.

  * •

pct_NH_Blk_alone_ACS_13_17: The percentage of the ACS population that indicate
no Hispanic origin and their only race “Black,African Am.,or Negro” or report
entries such as African American, Kenyan, Nigerian, or Haitian.

  * •

pct_ENG_VW_ACS_13_17: The percentage of all ACS occupied housing units where
no one ages 14 years and over speaks English only or speaks English “very
well”.

  * •

pct_Othr_Lang_ACS_13_17: The percentage of the ACS population aged 5 years and
over that speaks a language other than English at home.

  * •

pct_Diff_HU_1yr_Ago_ACS_13_17: The percentage of the ACS population aged 1
year and over that moved from another residence in the U.S.or Puerto Rico
within the last year.

  * •

avg_Tot_Prns_in_HHD_ACS_13_17: The average number of persons per ACS occupied
housing unit. This was calculated by dividing the total household population
in the ACS by the total number of occupied housing units in the ACS.

  * •

pct_Sngl_Prns_HHD_ACS_13_17: The percentage of all ACS occupied housing units
where a householder lives alone.

  * •

pct_Female_No_HB_ACS_13_17: The percentage of all ACS occupied housing units
with a female householder and no spouse of householder present.

  * •

pct_Rel_Under_6_ACS_13_17: The percentage of 2010 ACS family-occupied housing
units with a related child under 6 years old; same-sex couple households with
no relatives of the householder present are not included in the denominator.

  * •

pct_Vacant_Units_ACS_13_17: The percentage of all ACS housing units where no
one is living regularly at the time of interview; units occupied at the time
of interview entirely by persons who are staying two months or less and who
have a more permanent residence elsewhere are classified as vacant.

  * •

pct_Renter_Occp_HU_ACS_13_17: The percentage of ACS occupied housing units
that are not owner occupied, whether they are rented or occupied without
payment of rent.

  * •

pct_Owner_Occp_HU_ACS_13_17: The percentage of ACS occupied housing units with
an owner or co-owner living in it.

  * •

pct_Single_Unit_ACS_13_17: The percentage of all ACS housing units that are in
a structure that contains only that single unit.

  * •

Med_HHD_Inc_ACS_13_17: Median ACS household income for the tract.

  * •

Med_House_Value_ACS_13_17: Median of ACS respondents’ house value estimates
for the tract.

  * •

pct_HHD_Moved_in_ACS_13_17: The percentage of all ACS occupied housing units
where the householder moved into the current unit in the year 2010 or later.

  * •

pct_NO_PH_SRVC_ACS_13_17: The percentage of ACS occupied housing units that do
not have a working telephone and available service.

  * •

pct_HHD_No_Internet_ACS_13_17: Percentage of ACS households that have no
Internet access.

  * •

pct_HHD_w_Broadband_ACS_13_17: Percentage of ACS households that have
broadband Internet access.

  * •

pct_Pop_w_BroadComp_ACS_13_17: Percentage of people that live in households
that have both broadband Internet access and a computing device of any kind in
the ACS.

  * •

pct_URBANIZED_AREA_POP_CEN_2010: The percentage of the 2010 Census total
population that lives in a densely settled area containing 50,000 or more
people.

  * •

pct_MrdCple_HHD_ACS_13_17: The percentage of all ACS occupied housing units
where the householder and his or her spouse are listed as members of the same
household; does include same sex married couples.

  * •

pct_NonFamily_HHD_ACS_13_17: The percentage of all ACS occupied housing units