# Sparse Functional Linear Discriminant Analysis

Juhyun Park
Lancaster University, U.K. and ENSIIE, France
Jeongyoun Ahn
University of Georgia, U.S.A.
and Yongho Jeon
Yonsei University, South Korea Corresponding author: Juhyun Park, Department
of Mathematics and Statistics, Lancaster University, Lancaster LA1 4YF, U.K.
and ENSIIE & LaMME, Université Paris-Saclay, 91025 Évry, France. email:
<EMAIL_ADDRESS>
###### Abstract

Functional linear discriminant analysis offers a simple yet efficient method
for classification, with the possibility of achieving a perfect
classification. Several methods are proposed in the literature that mostly
address the dimensionality of the problem. On the other hand, there is a
growing interest in interpretability of the analysis, which favors a simple
and sparse solution. In this work, we propose a new approach that incorporates
a type of sparsity that identifies non-zero sub-domains in the functional
setting, offering a solution that is easier to interpret without compromising
performance. With the need to embed additional constraints in the solution, we
reformulate the functional linear discriminant analysis as a regularization
problem with an appropriate penalty. Inspired by the success of
$\ell_{1}$-type regularization at inducing zero coefficients for scalar
variables, we develop a new regularization method for functional linear
discriminant analysis that incorporates an $L^{1}$-type penalty, $\int|f|$, to
induce zero regions. We demonstrate that our formulation has a well defined
solution that contains zero regions, achieving a functional sparsity in the
sense of domain selection. In addition, the misclassification probability of
the regularized solution is shown to converge to the Bayes error if the data
are Gaussian. Our method does not presume that the underlying function has
zero regions in the domain, but produces a sparse estimator that consistently
estimates the true function whether or not the latter is sparse. Numerical
comparisons with existing methods demonstrate this property in finite samples
with both simulated and real data examples.

Keywords: Domain selection; Functional classification; Functional sparsity;
Interpretability; $L^{1}$ penalty; Linear discriminant analysis.

##  1 Introduction

We consider a classification problem when data are curves. Denote the
functional predictor by $X\in L^{2}(\mathcal{I})$, observable on a real
interval $\mathcal{I}$, and the class label by $Y$. Assuming that $X$ is a
member of two possible groups ($Y=0$ or $Y=1$), we seek a classification rule
that depends on a linear map, with unknown function $\beta$,

$F_{\beta}(X)=\int_{\mathcal{I}}X(t)\beta(t)\,\mbox{d}t.$

This defines a functional linear classification problem, where we wish to
determine $\beta(\cdot)$ in such a way that the linear map yields a good class
separation. As functional data models are inherently infinite-dimensional,
dimension reduction techniques are essential in constructing a solution.

The standard logistic regression framework for classification can be extended
with functional variables (Müller, 2005) as

$\log\frac{\mbox{pr}(Y=1\mid X)}{\mbox{pr}(Y=0\mid
X)}=a_{0}+\int_{\mathcal{I}}X(t)\beta(t)\,\mbox{d}t\,.$

Alternatively, one can try to directly extend the Bayes classifier. If $X$
were finite-dimensional, the Bayes rule is defined as maximizing the
conditional probability given by

$\mbox{pr}(Y=1\mid
X)=\frac{\pi_{1}f_{1}(X)}{\pi_{0}f_{0}(X)+\pi_{1}f_{1}(X)}\,,$

where $i=0,1,\pi_{i}$ is the probability of $X$ coming from group $i$ and
$f_{i}(X)$ is the marginal density of $X$ if it belongs to group $i$. An
optimal classification can then be achieved when one classifies new variable
$X$ to the class 1 if $\pi_{1}f_{1}(X)>\pi_{0}f_{0}(X)$. Consequently, the
main effort is devoted to estimating $f_{1}(X)$ and $f_{0}(X)$. If $X$ is
Gaussian, further simplifications can be made and such an approach is broadly
known as linear discriminant analysis. The main difficulty with extending the
Bayes rule to functional data is linked to the fact that the marginal
densities do not exist for functional data (Delaigle and Hall, 2010).
Nevertheless, approaches attempting to directly approximate the marginal
densities with dimension reduction techniques have proven useful, even in the
absence of well-defined target densities. For example, James and Hastie (2001)
use a Gaussian framework to develop regularization methods, while Bongiorno
and Goia (2016) and Dai et al. (2017) develop nonparametric approaches without
Gaussian assumptions.

It is well known that functional classification can achieve a perfect
classification, if the infinite-dimensionality is well exploited. This implies
that for the purpose of classification, it is not necessarily advantageous to
have a well-defined finite-dimensional representation. Delaigle and Hall
(2012) demonstrate such a phenomenon with a simple linear centroid classifier
using asymptotic analysis and suggest a practical representation using
components obtained from functional principal component analysis and partial
least squares. Kraus and Stefanucci (2018) propose $L^{2}$ regularization
methods to obtain the representation. Berrendero et al. (2018) further clarify
this phenomenon under a reproducing kernel Hilbert space framework for
Gaussian processes, suggesting an alternative finite dimensional
approximation.

While optimal performance is an important criterion to consider, the
increasing impact of statistical analysis on modern scientific investigations
has created the need to carefully consider interpretability of the outcomes of
the analysis. Some attempts have been made to address interpretability in
functional data, based on the idea that a simpler form of function is easier
to interpret, and thus more useful in practice. As the function is an
infinite-dimensional object, the formulation is often given in terms of a
basis function representation. Under this setting, three different approaches
have been proposed to construct a simpler form of functions. The first one is
to impose a constraint on the coefficients directly with an $\ell_{1}$ norm
(e.g., Zhou et al., 2012). Assuming that $\beta$ can be well approximated by a
finite number of basis functions, say $K$, this can be expressed as
$\beta(t)=\sum_{k=1}^{K}\alpha_{k}B_{k}(t)$ subject to
$\sum_{k=1}^{K}|\alpha_{j}|\leq C$ and thus encourages the coefficients to be
zero. The second approach is to limit the class of the functions $\\{B_{k}\\}$
in terms of their shape such as constant or strictly linear functions only
(e.g., Tian and James, 2013). A difficulty with a standard sparse
regularization for a function is that the penalty that encourages a sparse
representation of the function does not necessarily have a control over domain
selection: i.e., even if $\alpha_{k}=0$ for some $k$, $\beta(t)\neq 0$ for
$t\in I_{k}(t)=\\{t\in\mathcal{I}:B_{k}(t)\neq 0\\}$. The third one is to
limit the support of the function $\beta$ to include zero regions (James et
al., 2009; Zhou et al., 2013; Martin-Barragan et al., 2014; Lin et al., 2016;
Picheny et al., 2019). We wish to incorporate interpretability in the latter
notion of obtaining zero regions in the solution. However, as noted by Kneip
et al. (2016) and Roche (2018), the theoretical framework appropriate to deal
with a discrete notion of sparsity in high dimensions does not necessarily
offer an insight into a problem in an infinite dimensional setting.

A functional formulation on sparsity is relatively scarce. Wang and Kai (2015)
introduce the notion of functional sparsity, distinguishing global sparsity,
which relates to functional variable selection, from local sparsity, which
relates to domain selection with zero regions. Tu et al. (2020) develop a
regularization method to achieve simultaneous estimation of both types of
sparsity in a time varying functional regression setting and Lin et al. (2017)
propose an alternative regularization to achieve local sparsity in functional
linear regression. Both approaches rely on a clever grouping of the sparse
coefficients in the basis function representation. Hall and Hooker (2016)
study the issue of identifiability of domain selection problem in functional
linear regression and suggest a domain search strategy. Kraus and Stefanucci
(2018) follow a similar line.

In this work, we seek an alternative approach to functional linear
classification with a direct estimation method that addresses dimensionality,
optimality and interpretability. With the need to embed additional constraints
on the form of the solution, we reformulate the functional linear discriminant
analysis as a regularization problem with an appropriate choice of penalty
functions. So far, the penalty-based regularization methods have been used
mostly for either smoothness in regression (e.g., Cardot et al., 2003; Crambes
et al., 2009) or invertibility in classification (Kraus and Stefanucci, 2018),
based on an $L^{2}$-type penalty. The idea of sparsity as variable selection
with an $\ell_{1}$-type penalty is actively developed in the high-dimensional
setting with scalar variables, but much less for the infinite-dimensional
functional setting. We develop a new regularization method for functional
linear discriminant analysis with an $L^{1}$-type penalty in the functional
setting for inducing zero regions, as opposed to zero coefficients. We study
the underlying optimization problem in detail, showing that it has a well
defined solution that contains zero regions, achieving a functional sparsity
in the sense of domain selection. Furthermore, we show that our regularized
classifier asymptotically behaves similarly to the optimal classifier. We do
not presume, even in the case where the optimal classifier is well defined,
that the underlying projection function $\beta$ has a local sparse property
with true zero regions in the domain. The role of the $L^{1}$ penalty is to
provide a sparse estimator that consistently estimates the true function
whether or not the latter is sparse, enhancing interpretability without
compromising prediction performance. Unlike the approaches based on direct
dimension reduction techniques such as principal component analysis or partial
least squares, our proposed regularization method does not require assumptions
on the eigenvalue sequences of the covariance function or any other unknown
quantities. Numerical comparisons with existing methods demonstrate the
properties of the proposed estimator in finite samples with both simulated and
real data examples.

##  2 Methodology

###  2.1 Functional data framework

We use the following standard notation: for $f\in L^{p}(\mathcal{I})$, $1\leq
p<\infty$, the norm is defined by
$\|f\|_{p}=(\int_{\mathcal{I}}|f(t)|^{p}\,\mbox{d}t)^{1/p}$. When $p=\infty$,
$\|f\|_{\infty}=\inf\\{C:|f(t)|\leq C\mbox{ a.e. on }\mathcal{I}\\}$.

Assume that $X\in L^{2}(\mathcal{I})$ with the mean $\mu(t)$ and covariance
function $\mbox{cov}\\{X(s),X(t)\\}=\gamma(s,t),s,t,\in\mathcal{I}$. Define
the inner product in $L^{2}(\mathcal{I})$ by $\langle
f,g\rangle=\int_{\mathcal{I}}fg$ for $f,g\in L^{2}(\mathcal{I})$. Assuming
$\int_{\mathcal{I}}\int_{\mathcal{I}}\gamma(s,t)^{2}\,\mbox{d}s\mbox{d}t<\infty$,
define the covariance operator $\Gamma:L^{2}(\mathcal{I})\rightarrow
L^{2}(\mathcal{I})$ as

$\Gamma(\beta)(t)=\int_{\mathcal{I}}\gamma(s,t)\beta(s)\,\mbox{d}s\,,\quad
t\in\mathcal{I}.$

It is known that $\Gamma$ is a compact operator, and is Hilbert-Schmidt, that
admits a spectral representation given by

$\Gamma(\beta)=\sum_{j=1}^{\infty}\theta_{j}\langle\beta,\phi_{j}\rangle\phi_{j}\,,\qquad\theta_{j}\geq
0\,,\quad\theta_{j}\rightarrow 0\,\,\mbox{ as }j\rightarrow\infty\,,$

where $\theta_{j}$ and $\phi_{j}(\cdot)$ correspond to eigenvalues and
eigenfunctions of the covariance operator $\Gamma$, respectively. One of the
important properties of a compact operator is that it is not invertible unless
it has only finitely many distinct eigenvalues. At the same time, since
$\theta_{j}\rightarrow 0$ as $j\rightarrow\infty$, a finite rank approximation
to $\Gamma$ is also well understood, which has been the basis of many
regularization and dimension reduction techniques for functional data, notably
functional principal component analysis and its applications (e.g., Hall et
al., 2001). More details on theoretical foundations of functional data are
found in Hsing and Eubank (2015).

For our classification problem, we denote the variable in each group by
$X_{0}$ and $X_{1}$ and we make the following standard assumptions.

######  Assumption 1.

For $k=0,1$, $X_{k}$ is square integrable function on a compact interval
$\mathcal{I}$ with $\mu_{k}={E}(X_{k})$ and $\mu_{0}\neq\mu_{1}$ with common
covariance function $\gamma$. The mean functions and the covariance function
are continuous.

######  Assumption 2.

For $k=0,1$, ${E}(\|X_{k}\|_{2}^{4})<\infty$.

Without loss of generality, we assume that $\mu_{0}=0$ and $\mu_{1}=\delta$.

###  2.2 An optimal linear classifier

The optimality of the linear methods in the homoskedastic Gaussian scenario is
well established by Delaigle and Hall (2012) based on an asymptotic centroid-
based classifier. We briefly review their framework to motivate our
regularization method, which is introduced in the following section.

For given $X$, the linear classifier with $\beta$ can be defined in multiple
ways. We assume that $\beta\in L^{2}(\mathcal{I})$ is continuous. The class
assignment based on $\beta$ is done via

$T^{0}(X)=\big{(}\langle
X,\beta\rangle-\langle\delta,\beta\rangle\big{)}^{2}-\big{(}\langle
X,\beta\rangle\big{)}^{2}\,,$

which assigns to $X$ the group label $Y=1$ if $T^{0}(X)<0$ and $Y=0$ if
$T^{0}(X)>0$. Then, the misclassification error is
$\pi_{0}\mbox{pr}\\{T_{0}(X)<0\mid Y=0\\}+\pi_{1}\mbox{pr}\\{T_{0}(X)>0\mid
Y=1\\}$. If $X$ is Gaussian, the misclassification probability can be
expressed as

$err_{X}(\beta)=1-\Phi\left(\frac{|\langle\delta,\beta\rangle|}{2\langle\beta,\Gamma\beta\rangle^{1/2}}\right)\,,$
(1)

with its minimal error given by $1-\Phi(\|\Gamma^{-1/2}\delta\|_{2}/2)$
(Delaigle and Hall, 2012; Kraus and Stefanucci, 2018). It can be seen that an
optimal function $\beta$ that minimizes the misclassification error is

$\max_{\beta\neq
0}\frac{\langle\delta,\beta\rangle^{2}}{\langle\Gamma\beta,\beta\rangle}=\max_{\beta\neq
0}\frac{\langle\delta,\beta\rangle^{2}}{\mbox{var}(\langle
X,\beta\rangle)}\,.$

This is equivalent to the extended criterion for Fisher’s discrimination
analysis for functional data (Shin, 2008) defined by

$\max_{\beta\neq 0}\frac{\mbox{var}[{E}\\{F_{\beta}(X)\mid
Y\\}]}{{E}[\mbox{var}\\{F_{\beta}(X)\mid Y\\}]}=\max_{\beta\neq
0}\frac{\pi_{0}\pi_{1}\langle\delta,\beta\rangle^{2}}{\langle\Gamma\beta,\beta\rangle}\,.$
(2)

Attempting to directly solve (2) would lead to a form of generalized
eigenvalue problem. Equivalently, this can be formulated as

$\max_{\beta}\langle\delta,\beta\rangle^{2}\quad\mbox{ subject to
}\quad\langle\Gamma\beta,\,\beta\rangle-1=0.$ (3)

The solution can be derived from the Lagrangian formulation of (3)

$J_{\rho}(\beta)=\langle\delta,\,\beta\rangle^{2}-\rho\\{\langle\Gamma\beta,\,\beta\rangle-1\\}\,,$

for some constant $\rho$, which leads to

$\Gamma\beta=\delta\,.$ (4)

Hence, we see that the optimal linear classifier can be defined as the
solution to (4). However, since a bounded inverse of $\Gamma$ does not exist
in the infinite dimensional setting (e.g., Cardot et al., 2007), this equation
clearly illustrates that linear discrimination analysis for functional data is
an ill-posed inverse problem. When the covariance operator further satisfies
$\|\Gamma^{-1}\delta\|_{2}<\infty$, a unique classifier exists in
$(\mbox{Ker}(\Gamma))^{\bot}$ with the optimal error
$1-\Phi(\|\Gamma^{-1/2}\delta\|_{2}/2)$. When
$\|\Gamma^{-1}\delta\|_{2}=\infty$, an optimal solution does not exist, but
optimal classification can be achieved asymptotically along a non-convergent
path and perfect classification may be possible if
$\|\Gamma^{-1/2}\delta\|_{2}=\infty$. Therefore, unlike regression, it is not
necessary to assume the existence of a unique solution to obtain an
approximate solution, whose performance can mimic the optimal classifier
asymptotically (Delaigle and Hall, 2012; Kraus and Stefanucci, 2018). In the
following, we will denote the solution to (4) by $\beta_{0}$ when it exists
and is finite.

###  2.3 A regularized solution to discrimination

The primary purpose of regularization is to solve an ill-posed inverse
problem. We introduce our regularization method to solve the ill-posed inverse
problem (4) above. A sensible strategy would be to impose some constraints on
the solution set or a penalty term to an objective function. In order to
impose the functional equation (4), we introduce a corresponding minimization
problem. Specifically, viewing the discriminant equation in (4) as a type of
score equation leads to defining an objective function as

$\ell(\beta)=\frac{1}{2}\langle\Gamma\beta,\beta\rangle-\langle\delta,\beta\rangle\,.$
(5)

Hence the functional derivative of (5) with respect to $\beta$ yields (4). In
practice, where $\Gamma$ and $\delta$ are not available, replacing
$(\Gamma,\delta)$ by their empirical counterpart $(\Gamma_{n},\delta_{n})$
gives

$\ell_{n}(\beta)=\frac{1}{2}\langle\Gamma_{n}\beta,\beta\rangle-\langle\delta_{n},\beta\rangle\,.$

For example, the standard sample covariance operator and the sample mean could
be used for $\Gamma_{n}$ and $\delta_{n}$. The standard sample estimators
often lack smoothness (e.g., Cardot et al., 1999) so it is often desirable to
use their smoothed version by taking into account the underlying design
schemes (Zhang and Wang, 2016). Our formulation does not rely on the specific
choice of estimators as long as they are consistent. To incorporate additional
constraints, we consider the following optimization problem

$\min_{\beta\in\mathcal{H}}\\{\ell_{n}(\beta)+P(\beta)\\}\,,$ (6)

where $\mathcal{H}$ is an appropriate function space and $P(\beta)$ is a
penalty corresponding to the constraint. The type of constraint requires some
further consideration. As the underlying problem is defined in
$L^{2}(\mathcal{I})$, it seems natural to add a penalty that depends on the
$L^{2}$ norm of $\beta$ such as $\|\beta\|_{2}^{2}$ as in Kraus and Stefanucci
(2018), while a smoothness constraint on $\beta$ would lead to the $L^{2}$
penalty on the derivatives. Nevertheless, $L^{2}$ regularization in general
does not produce a sparse solution (e.g., Lin et al., 2017; Tu et al., 2020).
This phenomenon is well investigated in multivariate data, and the popularity
of the $\ell_{1}$ or $\ell_{0}$ penalty demonstrates the effectiveness of a
sparse solution in a broader context.

In order to introduce a type of sparsity in the discriminant function, we
propose to directly impose the functional norm constraints
$\|\beta^{\prime}\|_{2}^{2}\leq C_{1}$ and $\|\beta\|_{1}\leq C_{2}$. Although
other constraints are possible, it turns out that this combination of
constraints allows direct control over the $L^{2}$ norm. One might wonder
whether the derivative penalty is necessary in our sparse regularization. Our
experiences suggest that a combined norm give more stable results. We have
assumed compact support to simplify the statement in Assumption 1. Under the
assumption of compact support for $\beta$, the role of the derivative penalty
does not seem so significant with respect to the classification performance.
In general, since we only have $\|\beta\|_{1}\leq\|\beta\|_{2}$, the bound on
the $L^{1}$ norm alone is not sufficient to ensure the $L^{2}$ properties of
$\beta$ and the techniques we have developed below do not necessarily work. On
the other hand, we have (e.g., Gabushin, 1967; Li and Leoni, 2018) that

$\|\beta\|_{2}\leq\|\beta\|_{1}+\|\beta^{\prime}\|_{2}\,,$

which suggests that the proposed method is able to control the $L^{2}$ norm of
the function more effectively by regularizing both the derivative and the
$L^{1}$ norm of the function.

Taking into account the functional constraints in (6) leads to the following
objective function

$J_{n}(\beta)=\frac{1}{2}\langle\Gamma_{n}\beta,\beta\rangle-\langle\delta_{n},\beta\rangle+P(\beta)\,,$
(7)

where $P(\beta)=\lambda\|\beta\|_{1}+(\eta/2)\|\beta^{\prime}\|_{2}^{2}$ and
$\lambda$ and $\eta$ are tuning parameters, chosen from data.

In this formulation, the derivative penalty is a standard smoothness penalty.
The non-standard part is the $L^{1}$ penalty, which is an infinite-dimensional
counterpart of the $\ell_{1}$ sparsity penalty. As is demonstrated below, the
$L^{1}$ norm can be linked to a sparse solution in function space in terms of
domain selection (Li and Leoni, 2018), which is a special case of functional
sparsity (Wang and Kai, 2015). Arguably, the result would be easier to
interpret as the localized effects are automatically identified.

We note that it is possible to use the ridge penalty ($\|\beta\|_{2}^{2}$)
instead of the derivative norm. In fact, this choice makes the analysis of the
optimization problem much easier as it directly controls the $L^{2}$ norm. In
this case, the fact that
$\langle\Gamma_{n}\beta,\beta\rangle+\tau\|\beta\|_{2}^{2}\geq\tau\|\beta\|_{2}^{2}$
for all $\tau>0$ is sufficient to guarantee the existence of a solution.
Nevertheless, we have chosen the derivative norm that has the added
consideration of smoothness of the solution, since a smoother solution is
easier to interpret.

###  2.4 Properties of the regularized classifier in the population model

In this section we study the characteristics of the proposed regularized
solution to the infinite-dimensional convex optimization problem. In order to
simplify the discussion, we first consider the following optimization problem
that replaces the sample quantities by the population counterparts. Let

$J(\beta)=\frac{1}{2}\langle\Gamma\beta,\beta\rangle-\langle\delta,\beta\rangle+\lambda\|\beta\|_{1}+\frac{\eta}{2}\|\beta^{\prime}\|_{2}^{2}\,.$
(8)

The smoothness constraint with $\beta^{\prime}$ reduces the feasible set of
$\beta$ from $L^{2}(\mathcal{I})$ to a differentiable subspace
$H^{1}(\mathcal{I})=\\{\beta\in L^{2}(\mathcal{I}):\beta\mbox{ is absolutely
continuous},\beta^{\prime}\in L^{2}(\mathcal{I})\\}$, which is the Sobolev
space $H^{1}=W^{1,2}$ in standard notation (Brezis, 2011). Since
$L^{2}(\mathcal{I})\subset L^{1}(\mathcal{I})$, we consider
$\mathcal{H}=\\{\beta\in
H^{1}(\mathcal{I}):\|\beta\|_{1}+\|\beta^{\prime}\|_{2}<\infty\\}$ as the
space of feasible solutions, equipped with the norm given by
$\|\beta\|=\|\beta\|_{1}+\|\beta^{\prime}\|_{2}$. Then it is clear that
$\mathcal{H}$ is a convex set.

Under this setting, it can be shown that $J$ is convex, so that the existence
of a solution is guaranteed. Moreover, we can show that $J$ is strictly convex
so the solution is unique. We summarize the result in the following
proposition.

######  Proposition 1.

The function $J:\mathcal{H}\rightarrow\mathbb{R}$ defined in (8) has a global
minimizer $\tilde{\beta}\in\mathcal{H}$ such that

$J(\tilde{\beta})\leq J(\beta)$

for all $\beta\in\mathcal{H}$. Furthermore, the solution is unique.

In order to understand the property of the solution $\tilde{\beta}$,
additional characterizations of the solution beyond its existence are
necessary. Note that if the objective function is differentiable, the first
derivative (along with the second derivative) condition sufficiently
characterizes the solution and a Newton-type algorithm can be used to find it.
Non-differentiable functions require an alternative form to replace the
derivative condition. Following similar arguments in Reyes (2015, ch 6) and
Glowinski (1984, p. 70-71), we derive a set of conditions that $\tilde{\beta}$
satisfies in the following proposition, which could be used for developing an
optimization algorithm.

######  Proposition 2.

If $\tilde{\beta}$ is the minimizer of $J(\beta)$ in (8), there exists a
function $\alpha\in L^{2}(\mathcal{I})$ with $|\alpha(t)|\leq\lambda$ a.e. for
$t\in\mathcal{I}$ such that $\tilde{\beta}$ satisfies the following relation:

$\langle\Gamma\tilde{\beta},\beta\rangle+\eta\langle\tilde{\beta}^{\prime},\beta^{\prime}\rangle-\langle\delta,\beta\rangle+\langle\alpha,\beta\rangle=0\,,\quad\mbox{
for all }\quad\beta\in\mathcal{H}\,,$ (9)

and

$\left\\{\begin{array}[]{ll}\alpha(t)=\lambda&\mbox{ on
}\\{t\in\mathcal{I}:\tilde{\beta}(t)>0\\}\,,\\\ |\alpha(t)|\leq\lambda&\mbox{
on }\\{t\in\mathcal{I}:\tilde{\beta}(t)=0\\}\,,\\\ \alpha(t)=-\lambda&\mbox{
on }\\{t\in\mathcal{I}:\tilde{\beta}(t)<0\\}\,.\end{array}\right.$ (10)

The equations (9) and (10) are necessary conditions for the solution to
satisfy. Equation (9) means that the subgradient of the objective at the
minimizer $\tilde{\beta}$ contains zero. Furthermore, Proposition 2
demonstrates the domain selection property of the solution. In order to make
an analogy to the finite high dimensional case, consider a simplistic example
where $\Gamma=I$, the identity operator, and $\eta=0$ in (8). Note that, since
the identity operator is not compact, this is not a realistic example from the
functional data perspective but serves as a toy example. Then, the optimal
solution $\tilde{\beta}$ can be given explicitly as

$\tilde{\beta}(t)=\left\\{\begin{array}[]{ll}\delta(t)-\lambda&\mbox{ if
}\delta(t)>\lambda\,,\\\ 0&\mbox{ if }\delta(t)\in[-\lambda,\lambda]\,,\\\
\delta(t)+\lambda&\mbox{ if }\delta(t)<-\lambda\,.\end{array}\right.$

Since $\delta$ is continuous, the set
$\\{t\in\mathcal{I}:\tilde{\beta}(t)=0\\}$ will be a union of intervals and
then $\tilde{\beta}$ joins the zero intervals continuously at the boundary.
Hence the solution exhibits a thresholding behaviour, similar to the sparse
estimators in the finite dimensional case, but on the continuous domain, thus
justifying the notion of functional sparsity in the sense of domain selection.
When $\eta=0$, the optimization occurs in $L^{2}(\mathcal{I})$, so the
boundary points of the zero intervals are not necessarily differentiable. For
the general case, it is difficult to express the solution in an analytical
form but the same argument holds and when $\eta>0$ as in (9) and (10), the
function values at zero intervals and non-zero intervals are joined smoothly.

Figure 1: Illustration of the roles of the tuning parameters $\eta$ and
$\lambda$. Each panel plots the true $\beta_{0}(t)$, as well as
$\tilde{\beta}(t)$ with a varying tuning parameter.

(a) Setting 3, $\lambda=2^{-7}$

(b) Setting 3, $\eta=2^{-5}$

(c) Setting 5, $\lambda=2^{-7}$

(d) Setting 5, $\eta=2^{-5}$

Moreover, suppose that $\beta_{c}\in H_{0}^{1}(\mathcal{I})=\\{u\in
H^{1}(\mathcal{I}):u\mbox{ is zero at the boundary of }\mathcal{I}\\}$. Then,
(9) leads to

$\int_{\mathcal{I}}(\Gamma\tilde{\beta}-\delta+\alpha)\beta_{c}+\int_{\mathcal{I}}(\eta\tilde{\beta}^{\prime})\beta_{c}^{\prime}=0\,,\quad\mbox{
for all }\quad\beta_{c}\in H_{0}^{1}(\mathcal{I})\,.$ (11)

This implies (e.g., Brezis, 2011) that $\tilde{\beta}^{\prime}$, the
derivative of the solution, is also absolutely continuous and satisfies

$\eta\tilde{\beta}^{\prime\prime}=\Gamma\tilde{\beta}-\delta+\alpha\,,\quad
a.e.\,,$ (12)

which, together with (10), gives an additional necessary condition for the
solution.

The following observation helps us to determine the upper bound of the tuning
parameter $\lambda$. Suppose that $\delta\in L^{\infty}(\mathcal{I})$. Then,
we have $\langle\delta,u\rangle\leq\|\delta\|_{\infty}\|u\|_{1}$. In addition,
it can be shown that the solution $\tilde{\beta}$ satisfies

$\langle\Gamma\tilde{\beta},\tilde{\beta}\rangle+\eta\langle\tilde{\beta}^{\prime},\tilde{\beta}^{\prime}\rangle+\lambda\|\tilde{\beta}\|_{1}=\langle\delta,\tilde{\beta}\rangle\leq\|\delta\|_{\infty}\|\tilde{\beta}\|_{1}\,.$

It follows that

$\langle\Gamma\tilde{\beta},\tilde{\beta}\rangle+\eta\langle\tilde{\beta}^{\prime},\tilde{\beta}^{\prime}\rangle+(\lambda-\|\delta\|_{\infty})\|\tilde{\beta}\|_{1}\leq
0\,.$

Since the first two terms are non-negative, if
$\|\delta\|_{\infty}\leq\lambda$, then $\tilde{\beta}=0$. This gives a range
of $\lambda$ values for a non-zero solution and this depends on the magnitude
of $\delta$. Figure 1 illustrates the effect of the tuning parameters $\eta$
and $\lambda$ on $\tilde{\beta}$, under two different scenarios for
$\beta_{0}$: sparse and non-sparse. The curves in the top and bottom panels
respectively correspond to the simulation settings 3 and 5 in Section 3. We
observe that for fixed $\lambda$, larger $\eta$ gives smoother estimates and
for fixed $\eta$, larger $\lambda$ gives a more sparse solution.

Suppose that $\mbox{Ker}(\Gamma)=\\{0\\}$. If
$\|\Gamma^{-1}\delta\|_{2}<\infty$, then $\beta_{0}=\Gamma^{-1}\delta\in
L_{2}(\mathcal{I})$ is the unique solution to (4). Denote by $\tilde{\beta}$
the minimizer of $J$ in (8). The following proposition shows that the
regularized solution $\tilde{\beta}$ approximates the underlying model
solution to (4) whether the latter is finite or not. We write
$\tilde{\lambda}=(\lambda,\eta)$, and $\tilde{\lambda}\rightarrow 0$ means
that $\lambda\rightarrow 0$ and $\eta\rightarrow 0$.

######  Proposition 3.

Fix any $u\in L_{2}(\mathcal{I})$. If $\|\Gamma^{-1}\delta\|_{2}<\infty$, then
$\langle\tilde{\beta},u\rangle\rightarrow\langle\beta_{0},u\rangle$ as
$\tilde{\lambda}\rightarrow 0$. If $\|\Gamma^{-1}\delta\|_{2}=\infty$, then
$\|\tilde{\beta}\|\rightarrow\infty$, as $\tilde{\lambda}\rightarrow 0$.

Using Proposition 3 together with (1), we can study the misclassification
error of the regularized classifier. We show below that our regularized
classifier mimics the behaviour of the optimal classifier in terms of the
misclassification rate.

######  Proposition 4.

Assume that $X$ is Gaussian with covariance operator $\Gamma$. The
misclassification probability of the regularized classifier $\tilde{\beta}$,
denoted by $err_{X}(\tilde{\beta})$, converges to
$1-\Phi(\|\Gamma^{-1/2}\delta\|_{2}/2)$ as $\tilde{\lambda}\rightarrow 0$.

######  Remark 1.

It may be suspected that the two classes have different covariance functions,
say $\Gamma_{0}$ and $\Gamma_{1}$, in which case a full generalization of our
linear approach is not trivial. For example, a direct generalization of (2)
replaces the denominator by $\langle\tilde{\Gamma}\beta,\beta\rangle$ where
$\tilde{\Gamma}=\pi_{0}\Gamma_{0}+\pi_{1}\Gamma_{1}$, the pooled covariance
function. Then the same formulation as (8), with $\Gamma$ replaced by
$\tilde{\Gamma}$, follows. Alternatively, we can consider the two-dimensional
projection approach for quadratic discriminant analysis proposed by Gaynanova
and Wang (2019) in the multivariate setting. In our functional framework, we
would optimize the following two objective functions

$J_{i}(\beta)=\frac{1}{2}\langle(\Gamma_{i}\beta,\beta\rangle-\langle\delta,\beta\rangle+\lambda_{i}\|\beta\|_{1}+\frac{\eta_{i}}{2}\|\beta^{\prime}\|_{2}^{2},$

for $i=0,1$. The classification will be then made based on the two linear maps
with the respective solutions. An advantage of this method over standard
quadratic discriminant analysis is that the computational complexity is
essentially the same as the homoscedastic case being considered here. In
either case, the theoretical guarantee on the domain selection and
classification accuracy need to be carefully investigated in the functional
context.

###  2.5 Properties of the regularized classifier in the sample model

The development in the previous section justifies the use of the $L^{1}$
penalty in our formulation, by ensuring that we have a well-defined
regularization problem (Proposition 1) with the desired sparsity in the
solution (Proposition 2). We can show that these properties also hold for the
empirical estimator, denoted by $\hat{\beta}$, the minimizer of $J_{n}$
defined in (7).

In the sample case, our formulation of $J_{n}$ is based on two estimators
$\delta_{n}$ and $\Gamma_{n}$. Provided that the sample covariance operator
$\Gamma_{n}$ is non-negative definite and $\delta_{n}\in L^{2}(\mathcal{I})$,
it can be shown that Propositions 1 and 2 hold for $J_{n}$ and $\hat{\beta}$.
In fact, the proofs are essentially the same, thus we omit the formal
statements of the Corollaries here.

In order to study the limit of the empirical estimator $\hat{\beta}$ and the
corresponding misclassification probability, we make the following
assumptions. Assumption 3 says that some consistent estimators are available
for $\delta_{n}$ and $\Gamma_{n}$ (or $\gamma_{n}$) and Assumption 4 assumes
that the tuning parameters decrease not too rapidly as the sample size gets
large.

######  Assumption 3.

Let $(\delta_{n},\gamma_{n})$ be consistent estimators such that

$\|\gamma_{n}-\gamma\|_{\infty}=O_{p}(a_{n}),\qquad\|\delta_{n}-\delta\|_{\infty}=O_{p}(b_{n})\,,$

for some sequences $a_{n},b_{n}$ converging to 0.

######  Assumption 4.

Let $\lambda=\lambda_{n}\rightarrow 0,\eta=\eta_{n}\rightarrow 0$ with
$\lambda_{n}/\eta_{n}=O(1),a_{n}/\eta_{n}=o(1)$ and $b_{n}/\eta_{n}=o(1)$.

For sufficiently dense data as we consider here, the standard sample mean and
sample covariance function can be used with a root-$n$ convergence.
Alternatively, any preferred smoothing methods could be employed, then the
convergence rates are given as a function of smoothing parameters used, for
example, see Yao et al. (2005); Li and Hsing (2010) for local linear
smoothers. Normally the error for covariance estimation is larger than that
for the mean estimation ($a_{n}\geq b_{n}$). Zhang and Wang (2016) generalize
those results under a general design scheme and show that a root-$n$ rate can
be achieved for dense and ultra-dense data, in which case the uniform
convergence rates can be $a_{n}=b_{n}=(\log n/n)^{1/2}$.

######  Proposition 5.

Fix any $u\in L_{2}(\mathcal{I})$. Suppose that Assumptions 1-4 hold. If
$\|\Gamma^{-1}\delta\|_{2}<\infty$, then
$\langle\hat{\beta},u\rangle\rightarrow\langle\beta_{0},u\rangle$ in
probability as $n\rightarrow\infty$. On the other hand, if
$\|\Gamma^{-1}\delta\|_{2}=\infty$, then $\|\hat{\beta}\|\rightarrow\infty$ in
probability as $n\rightarrow\infty$.

######  Proposition 6.

Assume that $X$ is Gaussian with covariance operator $\Gamma$. Under
Assumptions 1-4, the misclassification probability of the regularized
classifier $\hat{\beta}$, denoted by $err_{X}(\hat{\beta})$, converges to
$1-\Phi(\|\Gamma^{-1}\delta\|_{2}/2)$ in probability as $n\rightarrow\infty$.

###  2.6 Numerical algorithms

A fundamental strategy to solve infinite-dimensional optimization problems
could be said to be either to discretize-and-optimize or to optimize-and-
discretize (Reyes, 2015). For example, Glowinski (1984) analysed numerical
methods based on the former approaches using piecewise linear and quadratic
approximations on $H^{1}(\mathcal{I})$, while Reyes (2015) includes examples
of the latter. Generally, the former is easier, while the latter may give a
more elegant solution. According to Proposition 2, the second approach
requires solving equations (9), or (12), and (10) simultaneously. Based on
these equations, Newton-like update steps for both $\tilde{\beta}$ and
$\alpha$ could be derived. However, as this involves further development of
the testable conditions at the iteration to be adapted to a proper function
space (Reyes, 2015), we leave it for future work.

Here we have used piecewise linear approximations in our implementation with
equidistant grid points. We use estimates from the method by Yao et al. (2005)
for $\Gamma_{n}$ and $\delta_{n}$ and evaluate the derivative by the finite-
difference approximation. This strategy, combined with the rather simple form
of the optimization problem, leads to a lasso-type problem, for which many
efficient optimization algorithms are available. We implement the coordinate
descent algorithm by Fu (2012) to solve our discretized problem.

###  2.7 Implications in functional linear regression

Although we started with a classification problem, our regularization method
with an $L^{1}$-type penalty can be motivated from a regression point of view.
Consider a functional regression problem with scalar response $Y$ and
functional predictor $X\in H=L^{2}(\mathcal{I})$. For simplicity, assume that
${E}(Y)=0$ and ${E}(X)=0$ and consider a linear regression given by

$Y=\int_{\mathcal{I}}\beta(t)X(t)\,\mbox{d}t+\epsilon\,,$

where ${E}(\epsilon)=0$ and $var(\epsilon)=\sigma^{2}$. As before, the
covariance operator of $X$ is denoted by $\Gamma$. In addition, let
$\Lambda:L^{2}(\mathcal{I})\rightarrow\mathbb{R}$ denote the covariance
operator between $X$ and $Y$, defined as

$\Lambda(u)=\int_{\mathcal{I}}{E}\\{X(s)Y\\}u(s)\,\mbox{d}s\,.$

Then, the population least squares criterion can be expressed as

$E\\{(Y-\langle\beta,X\rangle)^{2}\\}=\langle\Gamma\beta,\beta\rangle-2\Lambda(\beta)+const\,,$

from which it follows that if $\tilde{\beta}$ is the minimizer if and only if
it satisfies

$\langle\Gamma\tilde{\beta},u\rangle=\Lambda(u)\qquad\mbox{ for all }u\in
H\,.$ (13)

To be consistent with the earlier notation, write
$\Lambda(u)=\langle\Delta,u\rangle$ where $\Delta(\cdot)=E\\{X(\cdot)Y\\}$.
Then, (13) can be expressed as

$\langle\Gamma\tilde{\beta},u\rangle=\langle\Delta,u\rangle\qquad\mbox{ for
all }u\in H\,.$

Note that when $Y$ is binary with values in $\\{-1,1\\}$, $\Delta$ is equal to
$\delta$ in (4), the difference in the group means in $X$ between two classes.

The existence of the solution requires the so-called Picard condition (Cardot
et al., 2003; Hsing and Eubank, 2015)

$\sum_{j=1}^{\infty}\frac{\langle
E\\{X(\cdot)Y\\},\phi_{j}\rangle^{2}}{\theta_{j}^{2}}=\sum_{j=1}^{\infty}\frac{\langle\Delta,\phi_{j}\rangle^{2}}{\theta_{j}^{2}}<\infty\,,$

which is equivalent to $\|\Gamma^{-1}\delta\|_{2}<\infty$ in the case of
classification problems (c.f., Section 2.2).

The above discussion shows that even though the motivation differs, the
underlying problem for functional linear regression is equivalent to that of
linear classification. Hence, our methodology is equally applicable in a
regression setting and offers a regularization method to estimate the
coefficient function $\beta$. Given a sample of size $n$ of
$(Y_{i},X_{i}),i=1,\ldots,n$, denote the sample version of the operators of
$\Gamma$ and $\Lambda$ by $\Gamma_{n}$ and $\Lambda_{n}$, respectively,
defined as

$\Gamma_{n}(u)(t)=\frac{1}{n}\sum_{i=1}^{n}\langle X_{i},u\rangle
X_{i}(t)\,,\quad\Lambda_{n}(u)=\frac{1}{n}\sum_{i=1}^{n}\langle X_{i},u\rangle
Y_{i}=\langle\Delta_{n},u\rangle\qquad u\in L^{2}(\mathcal{I})\,.$

where $\Delta_{n}(\cdot)=(1/n)\sum_{i=1}^{n}X_{i}(\cdot)Y_{i}$. Assuming that
the true coefficient function $\beta$ is a smooth function, Cardot et al.
(2003) proposed a generalization of a ridge regression with an $L^{2}$
derivative penalty as

$\frac{1}{2}\langle\Gamma_{n}\beta_{K},\beta_{K}\rangle-\langle\Delta_{n},\beta_{K}\rangle+\frac{\eta}{2}\|\beta_{K}^{(m)}\|_{2}^{2}\,,$

where $\beta_{K}=\sum_{j=1}^{K}c_{k}B_{k}$ with $B_{k}$ is a B-spline basis
function and $\beta_{K}^{(m)},m\geq 1$ is the $m$th derivative of $\beta_{K}$.
Our analysis suggests that adding an $L^{1}$ penalty would lead to a sparse
estimator that consistently estimates the true function, whether or not the
latter is sparse.

######  Remark 2.

Although the theoretical framework differs, our main consideration of
interpretability is not limited to this type of problem. For example, semi-
parametric problems of single-index or multiple-index regression or
classification can be expressed as

$Y=F\left(\sum_{k}\langle X,\beta_{k}\rangle\right)+\epsilon$

with a known or unknown link function $F$. Then an $L^{1}$-type functional
regularization for $\beta_{k}$ can be developed with an appropriate loss
function. We leave such a generalization to our future work.

##  3 Numerical studies

###  3.1 Simulation models

We use simulated curves under six different underlying population scenarios to
assess the performance of the proposed method and compare it with existing
approaches. Settings 1 – 5 represent scenarios where the underlying structure
follows model (4) in Section 2.2 and Setting 6 is borrowed from Delaigle and
Hall (2012). The proposed approach, which we name SFLDA (Sparse Functional
Linear Discriminant Analysis), is compared with eight different approaches:
the non-sparse version of the proposed approach (FLDA) that sets the tuning
parameter $\lambda=0$ in (7); ridge FLDA (RFLDA) proposed by Kraus and
Stefanucci (2018); partial least squares (PLS) by Delaigle and Hall (2012);
the three Bayes methods (B-Gauss, B-NPD, B-NPR) in Dai et al. (2017); the
quadratic discriminant analysis (QDA) by Galeano et al. (2015); functional
logistic regression (Logistic) (Müller, 2005).

The first four settings have a sparse underlying discriminant function. The
sparse discriminant functions, $\beta_{1}(t)$ and $\beta_{2}(t)$, are
generated with $33$ cubic B-spline basis functions $B_{i,4}~{}(i=1,\dots,33)$
with the knots $\xi_{j}=j/30,j=0,\dots,30$, where $\xi_{0}$ and $\xi_{30}$ are
the boundary knots. We set $\beta_{1}(t)=0.2B_{5,4}(t)+0.2B_{28,4}(t)$ for
Settings 1 and 2, and $\beta_{2}(t)=0.2B_{5,4}(t)-0.2B_{24,4}(t)$ for Settings
3 and 4. For the non-sparse $\beta_{3}(t)$ for Setting 5, we use 8 cubic
B-spline basis and set

$\beta_{3}(t)=0.1B_{1,4}(t)-0.3B_{3,4}(t)-0.2B_{5,4}(t)+0.2B_{7,4}(t)-0.1B_{8,4}(t).$

For the covariance function $\Gamma$, we used the Matérn covariance function:

$\gamma(s,t)=\sigma^{2}\frac{2^{1-\nu}}{G(\nu)}\left\\{(2\nu)^{1/2}\frac{|s-t|}{\rho}\right\\}^{\nu}K_{\nu}\left\\{(2\nu)^{1/2}\frac{|s-t|}{\rho}\right\\},$

where $G$ is the gamma function, $K_{\nu}$ is the modified Bessel function of
the second kind, for which the R-function besselK was used, and the parameters
were set as $\sigma=1$, $\rho=0.2$, and $\nu=3$. We set the mean of the first
group as $\mu_{0}(t)=5t+\sum_{j=1}^{5}(c_{j}/j)\phi(t)$, where
$\phi(t)=2^{1/2}\sin(\pi jt)$,
$(c_{1},\dots,c_{5})=(2.19,-0.18,-0.19,-2.51,-0.56)$, and the second group
mean as $\mu_{1}(t)=\mu_{0}(t)+\delta(t)$. Here the mean difference function
$\delta(t)$ is determined by
$\delta(t)=\int_{\mathcal{I}}\Gamma(t,s)\beta_{k}(s)ds$, $k=1,2,3$.

The data were assumed to be available on a fine grid. We used $T=100$
equispaced grid points on $[0,1]$. For grid points $t_{j},j=1,\dots,T$ , let
$\tilde{\Gamma}$ be the $T\times T$ matrix with $(i,j)$ entry
$\Gamma(t_{i},t_{j})$, and let $m_{0}$ and $m_{1}$ be the $T\times 1$ vectors
with $i$th entry $\mu_{0}(t_{i})$ and $\mu_{1}(t_{i})$, respectively. Data
were generated by $x_{0}=m_{0}+\tilde{\Gamma}^{1/2}e_{0}$ and
$x_{1}=m_{1}+\tilde{\Gamma}^{1/2}e_{1}$ where $e_{0}$ and $e_{1}$ are the
vectors of independent noise for each class, generated from either a standard
normal for Settings 1, 3, and 5 or a $t$-distribution with $5$ degrees of
freedom for Settings 2 and 4.

As mentioned earlier, Setting 6 is borrowed from Delaigle and Hall (2012).
Each curve from the $i$th group is generated as
$x_{i}=\sum_{j=1}^{40}(j^{-1}Z_{ij}+\mu_{ij})\phi_{j}(t)$, where $Z_{ij}$ are
centred exponential variables,
$(\mu_{01},\ldots,\mu_{06})=(0,-0.5,1,-0.5,1,-0.5)$, and
$(\mu_{11},\ldots,\mu_{16})=(0,-0.75,0.75,-0.15,1.4,0.1)$.

We generated $n=100$ curves for each group to train each classifier and
independently generated 300 curves for each group to evaluate the
misclassification error. The tuning parameters in each method are chosen via
5-fold cross-validation within the training data based on the
misclassification error.

###  3.2 Results

Table 1:  Test error rates (%) for simulated data based on 100 repetitions. Standard errors are in the parentheses. Setting | SFLDA | FLDA | RFLDA | PLS | B-Gauss | B-NPD | B-NPR | Logistic | QDA
---|---|---|---|---|---|---|---|---|---
1 | 33.3 | 33.4 | 33.6 | 33.7 | 33.4 | 34.6 | 35.1 | 34.7 | 34.2
(0.22) | (0.22) | (0.24) | (0.23) | (0.24) | (0.25) | (0.25) | (0.27) | (0.28)
2 | 33.5 | 33.5 | 33.8 | 33.7 | 33.6 | 34.7 | 35.4 | 34.9 | 34.2
(0.21) | (0.20) | (0.20) | (0.20) | (0.20) | (0.22) | (0.24) | (0.25) | (0.22)
3 | 37.2 | 36.9 | 37.5 | 37.1 | 37.2 | 37.8 | 39.0 | 38.5 | 38.0
(0.23) | (0.20) | (0.20) | (0.21) | (0.22) | (0.25) | (0.29) | (0.27) | (0.25)
4 | 37.7 | 37.3 | 38.1 | 37.6 | 37.5 | 39.1 | 40.3 | 39.2 | 38.2
(0.24) | (0.22) | (0.24) | (0.23) | (0.24) | (0.27) | (0.30) | (0.28) | (0.26)
5 | 3.2 | 3.1 | 3.1 | 3.1 | 3.2 | 3.5 | 3.7 | 3.6 | 4.2
(0.08) | (0.07) | (0.07) | (0.07) | (0.08) | (0.10) | (0.10) | (0.14) | (0.08)
6 | 3.3 | 3.2 | 3.3 | 3.1 | 4.0 | 3.4 | 2.5 | 2.5 | 3.6
(0.08) | (0.07) | (0.07) | (0.07) | (0.08) | (0.08) | (0.10) | (0.10) | (0.14)
Table 2:  Average norm difference ($\|\hat{\beta}-\beta_{0}\|_{j}$) between the estimated discriminant function and the true function. Numbers are multiplied by $100$. Setting | $j$ | SFLDA | FLDA | RFLDA | PLS
---|---|---|---|---|---
1 | $1$ | 5.37 (0.19) | 7.97 (0.10) | 8.81 (0.13) | 8.20 (0.09)
$2$ | 9.64 (0.22) | 10.07 (0.11) | 11.29 (0.18) | 10.14 (0.10)
3 | $1$ | 5.84 (0.18) | 8.10 (0.11) | 9.00 (0.12) | 8.30 (0.11)
$2$ | 10.39 (0.22) | 10.40 (0.14) | 11.69 (0.18) | 10.26 (0.11)
5 | $1$ | 6.47 (0.14) | 4.15 (0.11) | 5.78 (0.17) | 4.50 (0.10)
$2$ | 8.21 (0.18) | 5.25 (0.13) | 7.23 (0.20) | 5.64 (0.11)

Table 1 lists test errors based on 100 repetitions. The proposed approach is
generally competitive in all settings, even though the differences among the
methods do not stand out. A main benefit of the proposed sparse estimation can
be seen in its ability to estimate zero regions in the discriminant function.
Table 2 displays the norm differences of the estimated discriminant function
$\hat{\beta}(t)$ by the proposed method with sparse functional penalty, its
non-sparse version, ridge regularization, and partial least squares, in
Settings 1, 3, and 5. The proposed method has the lowest error with respect to
both norms in the two sparse settings (1 and 3), while its non-sparse version
and partial least squares are better in the non-sparse setting (5). Also shown
in Figure 2 are the estimated curves for each of the four methods for Settings
1, 3, and 5, along with the median in solid black and the true discriminant
curve in a dot-dashed black line. The median curve of the estimators from the
proposed method resembles the true function the most in Settings 1 and 3,
successfully identifying the non-signal regions. The estimators with a ridge
penalty show larger variability than other methods, which suggests that a
ridge-type regularization in functional classification does not necessarily
reduce variance, unlike in ridge regression.

Figure 2: Discriminant functions from 100 repetitions of the simulation study
estimated by SFLDA, FLDA, RFLDA, and PLS, for Settings 1, 3, and 5. The
thicker black line in each plot represents the median of the estimated curves,
while the dot-dashed curve is the true discriminant function.

SFLDA

FLDA

RFLDA

PLS

SFLDA

FLDA

RFLDA

PLS

SFLDA

FLDA

RFLDA

PLS

Figure 3: For each real data example, 20 sample curves from each class are
shown along with average curves. Solid and dotted curves represent curves from
different classes.

(a) Tecator

(b) Wine

(c) Growth

(d) dGrowth

(e) Wheat

(f) Phoneme

###  3.3 Real data examples

Figure 4: Estimated discriminant curves for real data examples from 30
repetitions. The median curve is shown as a thicker black solid line.

Tecator - SFLDA

FLDA

RFLDA

PLS

Wine - SFLDA

FLDA

RFLDA

PLS

Growth - SFLDA

FLDA

RFLDA

PLS

Dgrowth - SFLDA

FLDA

RFLDA

PLS

Wheat - SFLDA

FLDA

RFLDA

PLS

Phoneme - SFLDA

FLDA

RFLDA

PLS

Table 3:  Mean test errors (%) from real data examples, based on 30 repetitions. Standard errors are in the parentheses. Data | SFLDA | FLDA | RFLDA | PLS | B-Gauss | B-NPD | B-NPR | Logistic | QDA
---|---|---|---|---|---|---|---|---|---
Tecator | 5.0 | 3.4 | 7.4 | 4.8 | 6.3 | 4.2 | 4.8 | 1.8 | 3.3
(0.44) | (0.41) | (0.68) | (0.46) | (0.48) | (0.42) | (0.59) | (0.27) | (0.42)
Wine | 10.3 | 10.0 | 9.2 | 9.2 | 11.1 | 10.4 | 9.9 | 8.7 | 11.6
(0.81) | (0.74) | (0.72) | (0.81) | (0.91) | (0.89) | (0.57) | (0.65) | (0.80)
Growth | 6.0 | 5.7 | 8.3 | 6.1 | 6.2 | 7.2 | 7.1 | 5.9 | 6.7
(0.39) | (0.39) | (0.43) | (0.32) | (0.44) | (0.36) | (0.40) | (0.48) | (0.40)
dGrowth | 7.0 | 6.4 | 7.4 | 6.3 | 7.4 | 8.1 | 7.6 | 7.1 | 9.5
(0.57) | (0.47) | (0.52) | (0.52) | (0.48) | (0.49) | (0.55) | (0.56) | (0.60)
Wheat | 0.1 | 0.0 | 0.1 | 0.0 | 0.1 | 0.0 | 0.2 | 0.2 | 0.2
(0.09) | (0.00) | (0.33) | (0.00) | (0.09) | (0.00) | (0.13) | (0.17) | (0.13)
Phoneme | 19.7 | 21.6 | 20.9 | 21.4 | 23.7 | 23.8 | 24.2 | 22.1 | 24.4
(0.53) | (0.59) | (0.53) | (0.58) | (0.60) | (0.57) | (0.60) | (0.62) | (0.65)

We tested the same set of the functional classification methods as in the
previous section for real data examples. The Tecator dataset (Delaigle and
Hall, 2012) has 240 curves of near infrared absorbance spectra (850 - 1050 nm)
of finely chopped meat, using a Tecator Infratec Food & Feed Analyser. Here
the groups are defined according to the fat content. The Wine spectra dataset
(Dai et al., 2017) contains 123 samples of mid infrared spectra (4000 - 400
cm-1), which are divided into two groups based on the alcohol content level.
The Growth dataset (Gasser et al., 1984; Sheehy et al., 1999) has height
growth curves of 112 boys and 120 girls from births to the 18th year. We also
analysed the velocity of the Growth curves. The Wheat data (Kalivas, 1997)
have near infrared spectra of 100 wheat samples, divided into two groups
according to the protein content. The Phoneme data are log-periodograms
constructed from digitized speech of two different sounds “aa” and “ao”, as
described in Hastie et al. (2009). Each panel of Figure 3 respectively shows
40 randomly selected curves from each data set, with each class shown with
different colors.

To estimate test errors, we split into 2/3 training and 1/3 testing data, and
tune each method in the same way as the simulation study, and repeat this
process 30 times. Mean test errors are displayed in Table 3. While there are
no meaningful differences among the test errors of these methods, we can see
that the proposed method achieves compatible accuracies with a much better
ability to find localized discriminant regions, as can be seen in Figure 4. In
particular, we have found that the selected regions of Tecator, Wine, and
Phoneme are consistent with findings from the spectroscopy literature.
Specifically, according to the body composition study (Conway et al., 1984),
fat has high absorbance around 930nm, the absorbance region for ethanol is
concentrated in 1200 - 850 cm-1 (Debebe et al., 2017), and it has been found
by Hastie et al. (1995) that the discriminating feature in phonemes are in the
low frequencies about 500 - 1000Hz, corresponding to the frequencies 16–32.
All these known regions are consistent with the estimates from the proposed
method in the figure.

## Acknowledgment

We are grateful to Aurore Delaigle and Xiongtao Dai for sharing the Matlab
codes for compared methods. We are also grateful to two referees, an associate
editor and the editor for careful reading and helpful suggestions. The third
author’s work was supported by Basic Science Research Program of the National
Research Foundation of Korea (NRF-2019R1A2C1005979) funded by the Korean
government.

## Supplementary material

Supplementary material includes additional details on relevant background
material and proofs, as well as an example of matlab codes for implementation.

## References

  * Berrendero et al. (2018) Berrendero, J. R., A. Cuevas, and J. L. Torrecilla (2018).  On the use of reproducing kernel Hilbert spaces in functional classification.  J. Am. Statist. Assoc. 113(523), 1210–1218.
  * Bongiorno and Goia (2016) Bongiorno, E. G. and A. Goia (2016).  Classification methods for Hilbert data based on surrogate density.  Comput. Stat. Data Anal. 99, 204–222.
  * Brezis (2011) Brezis, H. (2011).  Functional analysis, Sobolev spaces and partial differential equations.  Springer.
  * Cardot et al. (1999) Cardot, H., F. Ferraty, and P. Sarda (1999).  Functional linear model.  Stat. Probabil. Lett. 45(5), 11–22.
  * Cardot et al. (2003) Cardot, H., F. Ferraty, and P. Sarda (2003).  Spline estimators for the functional linear model.  Stat. Sinica 13, 571–591.
  * Cardot et al. (2007) Cardot, H., A. Mas, and P. Sarda (2007).  CLT in functional linear regression models.  Probab. Theory Rel. 138, 325–361.
  * Conway et al. (1984) Conway, J. M., K. H. Norris, and C. Bodwell (1984).  A new approach for the estimation of body composition: infrared interactance.  Am. J. Clin. Nutr. 40(6), 1123–1130.
  * Crambes et al. (2009) Crambes, C., A. Kneip, and P. Sarda (2009).  Smoothing spline estimators for functional linear regression.  Ann. Statist. 37(1), 35–72.
  * Dai et al. (2017) Dai, X., H.-G. Müller, and F. Yao (2017).  Optimal Bayes classifiers for functional data and density ratios.  Biometrika 104(3), 545–560.
  * Debebe et al. (2017) Debebe, A., M. Redi-Abshiro, and B. S. Chandravanshi (2017).  Non-destructive determination of ethanol levels in fermented alcoholic beverages using fourier transform mid-infrared spectroscopy.  Chemistry Central Journal 11(1), 1–8.
  * Delaigle and Hall (2010) Delaigle, A. and P. Hall (2010).  Defining probability density for a distribution of random functions.  Ann. Statist. 38(2), 1171–1193.
  * Delaigle and Hall (2012) Delaigle, A. and P. Hall (2012).  Achieving near perfect classification for functional data.  J. R. Statist. Soc. B 74(2), 267–286.
  * Fu (2012) Fu, W. J. (2012, Feb).  Penalized regressions: The bridge versus the lasso.  J. Comput. Graph. Stat. 7(3), 397–416.
  * Gabushin (1967) Gabushin, V. (1967).  Inequalities for the norms of a function and its derivatives in metric $l_{p}$.  Mathematical Notes of the Academy of Sciences of the USSR 1, 194–198.
  * Galeano et al. (2015) Galeano, P., E. Joseph, and R. E. Lillo (2015).  The Mahalanobis distance for functional data with applications to classification.  Technometrics 57(2), 281–291.
  * Gasser et al. (1984) Gasser, T., H. G. Müller, W. Köhler, L. Molinari, and A. Prader (1984).  Nonparametric regression analysis of growth curves.  Ann. Statist. 12, 210–229.
  * Gaynanova and Wang (2019) Gaynanova, I. and T. Wang (2019).  Sparse quadratic classification rules via linear dimension reduction.  J. Multivar. Anal. 169, 278–299.
  * Glowinski (1984) Glowinski, R. (1984).  Numerical methods for nonlinear variational problems.  Springer Series in Computational Physics. Springer-Verlag.
  * Hall and Hooker (2016) Hall, P. and G. Hooker (2016).  Truncated linear models for functional data.  J. R. Statist. Soc. B 78(3), 637–653.
  * Hall et al. (2001) Hall, P., D. S. Poskitt, and B. Presnell (2001).  A functional data-analytic approach to signal discrimination.  Technometrics 43(1), 1–9.
  * Hastie et al. (1995) Hastie, T., A. Buja, and R. Tibshirani (1995).  Penalized discriminant analysis.  Ann. Statist. 23, 73–102.
  * Hastie et al. (2009) Hastie, T., R. Tibshirani, and J. Friedman (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction (the 2nd ed.).  Springer.
  * Hsing and Eubank (2015) Hsing, T. and R. Eubank (2015).  Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators.  Wiley series in probability and statistics. Wiley.
  * James and Hastie (2001) James, G. M. and T. J. Hastie (2001).  Functional linear discriminant analysis for irregularly sampled curves.  J. R. Statist. Soc. B 63(3), 533–550.
  * James et al. (2009) James, G. M., J. Wang, and J. Zhu (2009).  Functional linear regression that’s interpretable.  Ann. Statist. 37, 2083–2108.
  * Kalivas (1997) Kalivas, J. H. (1997).  Two data sets of near infrared spectra.  Chemometrics and Intelligent Laboratory Systems 37(2), 255–259.
  * Kneip et al. (2016) Kneip, A., D. Poß, and P. Sarda (2016).  Functional linear regression with points of impact.  Ann. Statist. 44(1), 1–30.
  * Kraus and Stefanucci (2018) Kraus, D. and M. Stefanucci (2018).  Classification of functional fragments by regularized linear classifiers with domain selection.  Biometrika 106(1), 161–180.
  * Li and Leoni (2018) Li, T. and G. Leoni (2018).  ${L}^{1}$ regularization for compact support.  SIAM Undergraduate Research Online 11, 101–121.
  * Li and Hsing (2010) Li, Y. and T. Hsing (2010).  Uniform convergence rates for nonparametric regression and principal component analysis in functionallongitudinal data.  Ann. Statist. 38(6), 3321–3351.
  * Lin et al. (2017) Lin, Z., J. Cao, L. Wang, and H. Wang (2017).  Locally sparse estimator for functional linear regression models.  J. Comput. Graph. Stat. 26(2), 306–318.
  * Lin et al. (2016) Lin, Z., L. Wang, and J. Cao (2016).  Interpretable functional principal component analysis.  Biometrics 72, 846–854.
  * Martin-Barragan et al. (2014) Martin-Barragan, B., R. Lillo, and J. Romo (2014).  Interpretable support vetor machines for functional data.  Eur. J. Oper. Res. 232, 146–155.
  * Müller (2005) Müller, H.-G. (2005).  Functional modelling and classification of longitudinal data.  Scand. J. Stat. 32(2), 223–240.
  * Picheny et al. (2019) Picheny, V., R. Servien, and N. Villa-Vialaneix (2019).  Interpretable sparse sir for functional data.  Stat. Comput. 29, 255–267.
  * Reyes (2015) Reyes, J. D. L. (2015).  Numerical PDE-constrained optimization.  Springer briefs in optimization. Springer.
  * Roche (2018) Roche, A. (2018).  Variable selection and estimation in multivariate functional linear regression via the lasso.  Technical report, hal-01725351.
  * Sheehy et al. (1999) Sheehy, A., T. Gasser, L. Molinari, and R. H. Largo (1999).  An analysis of variance of the pubertal and midgrowth spurts for length and width.  Ann. Hum. Biol. 26(4), 309–331.
  * Shin (2008) Shin, H. (2008).  An extension of Fisher’s discriminant analysis for stochastic processes.  J. Multivar. Anal. 99, 1191–1216.
  * Tian and James (2013) Tian, T. S. and G. M. James (2013).  Interpretable dimension reduction for classifying functional data.  Comput. Stat. Data Anal. 57, 282–296.
  * Tu et al. (2020) Tu, C. Y., J. Park, and H. Wang (2020).  Estimation of functional sparsity in nonparametric varying coefficient models for longitudinal data analysis.  Stat. Sinica 29, 439–465.
  * Wang and Kai (2015) Wang, H. and B. Kai (2015).  Functional sparsity: Global versus local.  Stat. Sinica 25, 1337–1354.
  * Yao et al. (2005) Yao, F., H.-G. Müller, and J.-L. Wang (2005).  Functional data analysis for sparse longitudinal data.  J. Am. Statist. Assoc. 100(470), 577–590.
  * Zhang and Wang (2016) Zhang, X. and J.-L. Wang (2016).  From sparse to dense functional data and beyond.  Ann. Statist. 44(5), 2281–2321.
  * Zhou et al. (2013) Zhou, J., N.-Y. Wang, and N. Wang (2013).  Functional linear model with zero-value coefficient function at sub-regions.  Stat. Sinica 23, 25–50.
  * Zhou et al. (2012) Zhou, Y., R. Ogden, and P. Reiss (2012).  Wavelet-based LASSO in functional linear regression.  J. Comput. Graph. Stat. 21(3), 600–617.


# High-resolution land cover change from low-resolution labels: Simple
baselines for the 2021 IEEE GRSS Data Fusion Contest

Kolya Malkin
Yale University
&Caleb Robinson
Microsoft AI for Good
&Nebojsa Jojic
Microsoft Research
<EMAIL_ADDRESS>
###### Abstract

We present simple algorithms for land cover change detection in the 2021 IEEE
GRSS Data Fusion Contest [2]. The task of the contest is to create high-
resolution (1m / pixel) land cover change maps of a study area in Maryland,
USA, given multi-resolution imagery and label data. We study several baseline
models for this task and discuss directions for further research.

See https://dfc2021.blob.core.windows.net/competition-data/dfc2021_index.txt
for the data and https://github.com/calebrob6/dfc2021-msd-baseline for an
implementation of these baselines.

##  1 Introduction

We describe a possible road map towards estimating high-resolution (e.g., 1m)
land cover change, using high-resolution multitemporal input imagery and low-
resolution (e.g., 30m) noisy land cover labels for one or more time points.
The examples here use the data from the multitemporal semantic change
detection track of the 2021 IEEE Geoscience and Remote Sensing Society’s Data
Fusion Contest (DFC-MSD) [2]. The input imagery is from the National
Agriculture Imagery Program (NAIP) for the years 2013 and 2017, captured at 1m
resolution, and the available labels are from the 30m-resolution National Land
Cover Database (NLCD) for the years 2013 and 2016. The NLCD labels are derived
from Landsat imagery, and the contest also provides multiyear Landsat data for
the period between 2013 and 2017. All the data is limited to Maryland, USA.

The task exemplifies the situation commonly found worldwide. New imagery comes
in faster than high-quality high-resolution labels are being created. However,
some older, noisy, and low-resolution labels are often available, e.g., 30m
NLCD in the United States, or 500m MODIS land cover available worldwide. How
could one use machine learning to build models that predict high-resolution
change without a large set of high-resolution change examples? The hope
springs from the success of weakly supervised segmentation and label-super
resolution research, e.g. [8, 7], which demonstrated that it is possible to
train high-resolution label predictors for high-resolution input imagery using
regional supervision, where a large block of land (e.g. 30m$\times$30m) is
labeled with a single class, which often designates an area where land cover
follows particular mixing proportions. Weak supervision was the topic of the
2020 IEEE GRSS Data Fusion Contest [1].

Taking the next step to predict land cover _change_ could be as
straightforward as applying these techniques separately for different time
points and comparing the super-resolved labels to estimate change, and in this
paper we evaluate what such simple baselines could accomplish. We demonstrate
some simple lessons in terms of scope of training and model selection that
could aid in development of techniques specifically built for high-resolution
_change_ detection. Such techniques would presumably not use the multiyear
images in isolation, but analyze them jointly to take advantage of the fact
that most land cover _does not_ change from year to year.

The paper is organized as follows. We first describe the data in the DFC-MSD
task and illustrate the types of changes that are being scored in the contest.
Then we study training neural networks to approximate functions $y=f(x)$ where
$x$ are high-resolution (1m) image patches and $y$ are upsampled $30$m labels,
and discuss how the choices for model complexity and scope of training data
can affect the results. In particular, a highly expressive and well-trained
function $f$ could learn to predict blurry labels, because the targets $y$
cone in $30\times 30$ blocks, while very simple color and texture models,
e.g., [7, 9], tend to super-resolve them at 1m, as they are unable to assign
the same label to $30\times 30$ blocks of pixels with varying color and
texture. Finally, we discuss possible approaches that could use all the data
in concert to build more accurate prediction models.

All data can be downloaded from the addresses listed at
https://dfc2021.blob.core.windows.net/competition-data/dfc2021_index.txt;
alternative fast download instructions can be found on the DFC-MSD webpage
[2]. Example code is available in the accompanying code repository at
https://github.com/calebrob6/dfc2021-msd-baseline.

##  2 The change detection task

###  2.1 Data

The input data in the DFC-MSD comprises 9 layers covering the study area of
the state of Maryland in the United States ($\sim$35,000 km2). All layers are
upsampled from their native resolutions to 1m / pixel and provided as 2250
aligned _tiles_ of dimensions not exceeding $4000\times 4000$. The layers are
the following:

  1. (1)

NAIP (2 layers): 1m-resolution 4-band (red, green, blue, and near-infrared)
aerial imagery from the US Department of Agriculture’s National Agriculture
Imagery Program (NAIP) from two points in time: 2013 and 2017.

  2. (2)

Landsat (5 layers): 30m-resolution 9-band imagery from the Landsat-8 satellite
from five time points: 2013, 2014, 2015, 2016, and 2017. Each of these images
is a median-composite from all cloud and cloud-shadow masked surface-
reflectance scenes intersecting Maryland.

  3. (3)

NLCD (2 layers): 30m-resolution coarse land cover labels from the US
Geological Survey’s National Land Cover Database [5], in 15 classes (see Table
1) from two time points: 2013 and 2016. These labels were created in a semi-
automatic way, with Landsat imagery as the principal input.

An example of the input layers is shown in Figure 1.

Figure 1: Four of the input data layers in the DFC-MSD. Above:  NAIP imagery
from 2013 and 2017. Below:  Coarse NLCD land cover labels from 2013 and 2016.
See Table 1 for the color scheme.

###  2.2 Task

The goal of the DFC-MSD is to classify pixels in the study area as belonging
to classes describing land cover change between the times when NAIP 2013 and
2017 aerial photographs were taken. Specifically, we wish to detect the loss
or gain of four land cover classes in a simplified scheme based on that of the
NLCD: water, tree canopy, low vegetation, and impervious surfaces, as well as
the absence of land cover change. An approximate correspondence between NLCD
classes and the target classes is shown in Table 1.

|  |  |  |  | Approximate class freq.
---|---|---|---|---|---
| NLCD class |  |  | Target class | W% | TC% | LV% | I%
| Open Water |  |  | water | 98 | 2 | 0 | 0
| Developed, Open Space |  |  | – | 0 | 39 | 49 | 12
| Developed, Low Intensity |  |  | – | 0 | 31 | 34 | 35
| Developed, Medium Intensity |  |  | impervious | 1 | 13 | 22 | 64
| Developed High Intensity |  |  | impervious | 0 | 3 | 7 | 90
| Barren Land (Rock/Sand/Clay) |  |  | – | 5 | 13 | 43 | 40
| Deciduous Forest |  |  | tree canopy | 0 | 93 | 5 | 0
| Evergreen Forest |  |  | tree canopy | 0 | 95 | 4 | 0
| Mixed Forest |  |  | tree canopy | 0 | 92 | 7 | 0
| Shrub/Scrub |  |  | tree canopy | 0 | 58 | 38 | 4
| Grassland/Herbaceous |  |  | low vegetation | 1 | 23 | 54 | 22
| Pasture/Hay |  |  | low vegetation | 0 | 12 | 83 | 3
| Cultivated Crops |  |  | low vegetation | 0 | 5 | 92 | 1
| Woody Wetlands |  |  | tree canopy | 0 | 94 | 5 | 0
| Emergent Herbaceous Wetlands |  |  | tree canopy | 8 | 86 | 5 | 0
Table 1: The approximate correspondence between the NLCD classes and the four
target classes. The last four columns show _estimated_ frequencies of target
classes of pixels labeled with each NLCD class, computed from a high-
resolution land cover product, specific to Maryland, in a similar class scheme
[3]. The colors next to NLCD classes are those used in Figure 1. The two
colors next to the target classes are those used to show loss and gain
respectively in Figures 2 and 4.

We require that any pixel classified as a gain of one land cover class must be
simultaneously classified as a loss of another class. Thus, the task is
equivalent to creating two aligned land cover maps in the four target classes
for the two time points and computing their disagreement: pixels assigned
class $c_{1}$ in 2013 and $c_{2}$ in 2017 are considered loss of $c_{1}$ and
gain of $c_{2}$; pixels in which the maps agree are considered no change. See
Figure 2 for examples of the desired outputs. For example, in the second row,
a part of the forest is cut (loss of tree canopy). In its place we see some
gain in impervious surfaces (buildings and roads), as well as some cleared
fields, which are placed in the low vegetation category. In a few places, new
trees were planted, resulting in no change in land cover between the two time
points. In the first row, where forest is also cut, most of it is still
undeveloped and consists of cleared fields. The undisclosed evaluation data
consists of many small areas of interest such as the ones in the figure, as
well as a set of areas where there has been no change. Ambiguous areas are
blocked from evaluation. These include shadows and areas which humans have
difficulty recognizing.

This problem is made difficult by the very different spectral characteristics
of the two NAIP layers (due to differences in sensors, seasons, and time of
day when they were created), as well as the lack of high-resolution label
data: the only available training labels are the NLCD layers, which are
coarse, created for different points in time than NAIP, and have a different
class scheme.

Figure 2: Examples of paired NAIP 2013 and NAIP 2017 imagery and the desired
predictions. See Table 1 for the color scheme. Pixels with no change are shown
in black.

###  2.3 Contest description

The DFC-MSD is organized in two phases: the _validation_ phase, in which
contestants have 100 attempts to submit their proposed land cover change maps
for a subset of the study area to an evaluation server and receive feedback on
their performance, and the _test_ phase, in which contestants have 10 attempts
to submit change maps for a different subset of the study area.

For the validation phase, a set of 50 out of the 2250 tiles were selected for
evaluation. Within these tiles, regions of interest were identified and high-
resolution change maps created by the contest organizers. Some regions of
interest contain pixels with land cover change, while others have only the no
change label. While participants submit predictions for all 50 tiles,
evaluation is performed only over the areas of interest; predictions for
pixels outside areas of interest are ignored. The class statistics for the
validation set are shown in Table 2.

The scoring metric is the average intersection over union (IoU) between the
predicted and ground truth labels for eight classes (loss and gain of each of
the four target classes, excluding no change). (Recall that the IoU for a
class $c$ is defined as

$\frac{\text{\\# pixels labeled $c$ in prediction \emph{and} in ground
truth}}{\text{\\# pixels labeled $c$ in prediction \emph{or} in ground
truth}},$ (1)

or, equivalently, $\frac{f}{2-f}$, where $f$ is the F1 score for class $c$.)
Notice that while IoU for the no change class is not computed, predictions of
change for pixels with no change still contribute to the denominator in (1)
and decrease the score. The evaluation server outputs the IoU for each class
upon submission.

For the test phase, a different set of 50 tiles and regions of interest is
selected. The submission and evaluation formats remain the same.

Class |  water |  tree canopy |  low vegetation |  impervious
---|---|---|---|---
loss |  0.70% |  8.12% |  10.20% |  1.33%
gain |  0.86% |  1.05% |  7.55% |  10.90%
Table 2: Class statistics for the validation set. The total number of pixels
in regions of interest is 13,496,574, or about 1.7% of the total area of the
validation tiles; 59.29% are labeled as no change.

##  3 Methods

###  3.1 Baselines

We study four baseline algorithms, each of which takes as input the 2013 and
2017 NAIP layers and the 2013 and 2016 NLCD layers.

  1. (1)

NLCD diff: We use solely the NLCD layers to assign each pixel a high-
resolution class according to the second column of Table 1; the Barren Land
and Developed, Low Intensity classes are assigned impervious and the
Developed, Open Space class is assigned low vegetation. This results in two
(coarse) land cover maps in the four target classes. The map derived from NLCD
2013 is used as the prediction for 2013; the map derived from NLCD 2016 is
used as the prediction for 2017. Land cover change is computed from these two
maps.

  2. (2)

We train small fully convolutional neural networks (5 layers of 64-filter
$3\times 3$ convolutions with ReLU activation and a logistic regression layer;
151k parameters) to predict the NLCD labels from NAIP imagery. Two FCNs are
trained: one for NAIP 2013, targeting NLCD 2013, and one for NAIP 2017,
targeting NLCD 2016.

We found the following scheme to work best for creating target labels from the
(probabilistic) FCN predictions: the predicted probabilities (at each pixel)
of the classes Developed, Open Space, Developed, Low Intensity, and Barren
Land, which do not have a dominant high-resolution target class, are set to 0,
and the probabilities are renormalized. The remaining classes are then mapped
to the four target labels according to the second column of Table 1: the
output probability of each target label is the sum of the probabilities of all
NLCD classes mapping to it. The target class with highest probability is then
taken as the prediction. The land cover change is computed from the two
resulting prediction layers.

    1. (a)

FCN / tile:  Two separate FCNs are trained for each of the validation tiles.
Overfitting models with a small receptive field to the coarse NLCD labels and
evaluating them results in a super-resolution of these labels: the FCNs cannot
detect the 30m block structure in the NLCD labels they were trained to
predict, so colors and textures are mapped to the classes of the low-
resolution blocks in which they are most likely to appear.

    2. (b)

FCN / all:  The same as (a), but a single FCN is trained on the entire
dataset. This may give an advantage due to the much larger training data
offering more potential for generalization. However, it may also hurt
predictions if pixels with similar appearance tend to belong to different
classes in distant parts of the study area.

  3. (3)

U-Net / all:  The same as (2b), but with a model of the U-Net family [10],
specifically, a U-Net with a ResNet-18 [4] encoder structure where the first 3
blocks in the ResNet are used in the downsampling path, and convolutions with
128, 64, and 64 filters are used in the corresponding upsampling layers (1.2m
total parameters). We use an implementation from the PyTorch Segmentation
Models library111https://github.com/qubvel/segmentation_models.pytorch; see
the accompanying GitHub repo for details.

The IoU scores for each of these models are shown in Table 3. Some examples of
their predictions appear in Figure 4.

Algorithm |  $-$W |  $-$TC |  $-$LV |  $-$I |  $+$W |  $+$TC |  $+$LV |  $+$I | avg.
---|---|---|---|---|---|---|---|---|---
NLCD diff | 0.148 | 0.167 | 0.282 | 0.014 | 0.031 | 0.001 | 0.106 | 0.362 | 0.139
FCN / tile | 0.641 | 0.436 | 0.407 | 0.091 | 0.255 | 0.073 | 0.302 | 0.528 | 0.342
FCN / all | 0.584 | 0.644 | 0.601 | 0.400 | 0.292 | 0.181 | 0.607 | 0.716 | 0.503
U-Net / all | 0.325 | 0.484 | 0.476 | 0.307 | 0.237 | 0.205 | 0.342 | 0.517 | 0.361
Table 3: Baseline IOUs. $-$C and $+$C denote loss and gain of class C,
respectively.

###  3.2 Discussion

Figure 3: Distribution of changes between 2013 and 2016 NLCD layers. Above:
Distribution of 2013 classes (left axis) and distribution of 2016 classes for
pixels of each 2013 class. Below:  The same, restricted to the set of pixels
whose labels differ between 2013 and 2016.

From Table 3, we see that NLCD difference alone is a poor predictor of land
cover change. There may be several reasons for this. First, this algorithm
gives only coarse change predictions, in 30m blocks, that miss fine details;
many changes are in small, thin regions (such as new houses or roads). Second,
the 2013 and 2016 NLCD layers were created in a joint pipeline in which
certain kinds of change involving the Developed classes were explicitly
forbidden [5] (see Figure 3). Third, the NLCD difference simply misses some
changes (e.g., in the lower right corner of Figure 1), either because of
errors made by the algorithm that produced the NLCD labels, or because NLCD
and NAIP represent different time points.

Interestingly, fitting FCNs to predict the coarse NLCD labels is a strong
baseline and an efficient way to super-resolve the NLCD labels to 1m
resolution. In some cases, the FCN trained on a single tile simply ignores
changes in areas that have no change in NLCD: it has overfit to the local
correspondences of (incorrect) NLCD labels and textures that appear. For
example, in Figure 4(d), it predicts the construction area as tree canopy –
the label suggested by both 2013 and 2016 NLCD labels. However, in other
cases, the FCN trained on one tile correctly classifies areas that appear
ambiguous to other models: in Figure 4(e), the silty part of the body of water
is predicted correctly as water in both years, as suggested by the Open Water
label in both NLCD layers, and thus shown as no change. However, to the FCN
trained on all tiles, the local information that this shade of water is indeed
water is trumped by evidence from other tiles that this color and texture
tends to be low vegetation, resulting in an incorrect prediction.

Differences between class appearances across geographic regions motivated
prior work on alternative loss functions for weakly supervised segmentation
(label super-resolution) [8]; it has also been found that tile-by-tile
algorithms sometimes perform better than those trained on a large dataset [7]
– in segmentation problems with wide variation of appearances over space, it
is sometimes the case that “less is more”. We believe that methods that
combine global models with tile-by-tile analysis warrant further study.

We also attempted to train models targeting _high-resolution_ labels from the
Chesapeake Conservancy land cover dataset [3], which use a class scheme
slightly different from ours. While the use of these labels is not allowed in
the DFC-MSD, as they are not available outside of a small geographic region,
we found that they do not improve the predictions: surprisingly, models
trained on these high-resolution labels over the entire study area perform no
better than those trained on coarse NLCD labels.

The last two rows of Table 3 suggest that sometimes “less is more” when it
comes to neural architectures as well. The U-Net has about $80\times$ more
parameters than the FCN. Figure 5 shows the raw predictions of NLCD labels
from the two models, both trained on the entire dataset. The FCN, which has a
small ($5\times 5$) receptive field, is unable to learn that NLCD labels come
in $30\times 30$ blocks, while the U-Net has a large receptive field and
learns to make blurry predictions. (In fact, the winning weak supervision
entry in the 2020 Data Fusion Contest [9] used very simple texture/color
mixture models to super-resolve MODIS labels, followed by smoothing with a
small FCN.)

We also found that targeting the same year of labels for both 2013 and 2017
networks does not change the results much, as the models primarily learn their
textures from the vast area where there is no change anyway. In fact, we
should note that training separate models for different time points is a
starting point, but we expect that the best approaches will analyze imagery
jointly, since most pixels do not change their class from one year to the
next. Such situations are often exploited in computer vision, e.g., in
background subtraction or co-segmentation. The probabilistic index map model
[6], for example, learns to assign a distribution over labels to each pixel in
an aligned image stack assuming that within each image assignments are
consistent, but the palette of assigned colors (or textures) to each index can
freely change.

(a) |
---|---
(b) |
(c) |
(d) |
(e) |
Figure 4: Examples of land cover change predictions by baseline models. Figure
5: Outputs of a FCN and a U-Net trained on the entire dataset to predict NLCD
2013 classes (bottom left) from NAIP 2013 imagery (top left).

##  4 Conclusion

In this paper, we have described simple algorithms for land cover change
detection and raised questions for future research to answer: How can machine
learning models optimally combine noisy, low-resolution labels for detection
and classification of change in high-resolution images? How should they strike
the balance between local analysis and models that work over a geographically
diverse study area? How can low-resolution (Landsat) imagery aid in change
detection? The 2021 IEEE GRSS Data Fusion Contest will help to answer these
questions, fostering innovation in weakly supervised segmentation and change
detection in this important application domain.

## References

  * [1] 2020 IEEE GRSS Data Fusion Contest, 2020.  http://www.grss-ieee.org/community/technical-committees/data-fusion/2020-ieee-grss-data-fusion-contest/.
  * [2] 2021 IEEE GRSS Data Fusion Contest Track MSD, 2021.  http://www.grss-ieee.org/community/technical-committees/data-fusion/2021-ieee-grss-data-fusion-contest-track-msd/.
  * [3] Chesapeake Conservancy.  Land Cover Data Project, 2017.  https://chesapeakeconservancy.org/wp-content/uploads/2017/01/LandCover101Guide.pdf.
  * [4] K. He, X. Zhang, S. Ren, and J. Sun.  Deep Residual Learning for Image Recognition.  Computer Vision and Pattern Recognition (CVPR), 2016.
  * [5] C. Homer, J. Dewitz, S. Jin, G. Xian, C. Costello, P. Danielson, L. Gass, M. Funk, J. Wickham, S. Stehman, R. Auch, and K. Riitters.  Conterminous United States land cover change patterns 2001–2016 from the 2016 National Land Cover Database.  ISPRS Journal of Photogrammetry and Remote Sensing, 162:184–199, 2020.
  * [6] N. Jojic and Y. Caspi.  Capturing image structure with probabilistic index maps.  Computer Vision and Pattern Recognition (CVPR), 2004.
  * [7] N. Malkin, A. Ortiz, and N. Jojic.  Mining self-similarity: Label super-resolution with epitomic representations.  European Conference on Computer Vision (ECCV), 2020.
  * [8] N. Malkin, C. Robinson, L. Hou, R. Soobitsky, J. Czawlytko, D. Samaras, J. Saltz, L. Joppa, and N. Jojic.  Label super-resolution networks.  International Conference on Learning Representations (ICLR), 2019\.
  * [9] C. Robinson, N. Malkin, L. Hu, B. Dilkina, and N. Jojic.  Weakly supervised semantic segmentation in the 2020 IEEE GRSS Data Fusion Contest.  International Geoscience and Remote Sensing Symposium (IGARSS), 2020\.
  * [10] O. Ronneberger, P. Fischer, and T. Brox.  U-Net: Convolutional networks for biomedical image segmentation.  International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI), 2015.