# Feature Selection for Vertex Discriminant Analysis

Alfonso Landeros2, Tong Tong Wu1, Kenneth Lange2,3,4

Department of Biostatistics and Computational Biology1
University of Rochester
Rochester, New York 14642-0630

Departments of Computational Medicine2,
Human Genetics3, and Statistics4
University of California
Los Angeles, CA 90095
Phone: 310-206-8076
E-mail<EMAIL_ADDRESS>

Research supported in part by
USPHS grants GM53275 and HG006139
and NSF grant CCF-1934962

###### Abstract

We revisit vertex discriminant analysis (VDA) from the perspective of proximal
distance algorithms. By specifying sparsity sets as constraints that directly
control the number of active features, VDA is able to fit multiclass
classifiers with no more than $k$ active features. We combine our sparse VDA
approach with repeated cross validation to fit classifiers across the full
range of model sizes on a given dataset. Our numerical examples demonstrate
that grappling with sparsity directly is an attractive approach to model
building in high-dimensional settings. Applications to kernel-based VDA are
also considered.

Support vector machines (SVM) stand out in the universe of classification
algorithms for their simplicity, interpretability, speed, and theoretical
guarantees (Smola and Schölkopf, 1998; Vapnik, 1999). The circle of SVM
algorithms have been bolstered by fundamental advances in statistical theory
and computer hardware but also challenged by the growing size of modern
datasets. The Achilles heal of SVMs is their limitation to binary
classification. Strategies to overcome this limitation, most notably one-
versus-rest and one-versus-one, rely on multiple SVMs to discriminate between
classes. This necessarily increases time to fit a model, obscures statistical
interpretation, and violates the principle of parsimony. Vertex discriminant
analysis (VDA) (Lange and Wu, 2008) addresses this particular weakness of
SVMs. Our subsequent paper (Wu and Lange, 2010) derives a fast cyclic
coordinate descent algorithm and attacks feature selection by a combination of
Euclidean group and lasso penalties. Our second followup paper (Wu and Wu,
2012) generalizes VDA to nonlinear discriminant analysis by the kernel trick.

The current paper revisits the problem of feature selection in VDA from the
perspective of distance to set penalties. This allows us to incorporate a
spectrum of constraints, with particular emphasis on sparsity sets. The
experiments recorded here demonstrate that the direct control of sparsity is
highly effective in building parsimonious models on both simulated and real-
world datasets.

##  1 Proposed Method

As a prelude to introducing a generic method for feature constraints, let us
review VDA with $\epsilon$-insensitive pseudo-distances and derive a quadratic
surrogate that both simplifies the associated empirical loss and converts
minimization to iterative least squares. Next, we present a distance penalized
version of VDA that reduces constrained minimization to unconstrained
minimization. In particular, we focus on distance to sparsity sets as a
mechanism for directly enforcing parsimony. We then present two competitive MM
algorithms designed with high-dimensional optimization in mind. Finally, we
comment on practical issues, including initialization strategies, convergence
criteria, and hyperparameter tuning.

Before delving into specifics, let us clarify our notation. We restrict our
attention to problems formulated over Euclidean space $\mathbb{R}^{n}$.
Vectors are written in lowercase boldface, for instance $\boldsymbol{y}$, and
matrices are written in uppercase boldface, for instance $\boldsymbol{X}$. The
notation $\|\cdot\|$ indicates the Euclidean norm on the ambient space. In
some cases we decorate a norm with a subscript, for example the Frobenius norm
$\|\cdot\|_{F}$, to avoid ambiguity. The superscript $\top$ indicates a vector
or matrix transpose. The lowercase letter $k$ denotes the number of selected
parameters; equivalently, the lowercase letter $s=1-k/p$ is reserved for the
sparsity level. The maximum number of parameters $p$ is understood from
context. Finally, the letters $m$ and $t$ are reserved for inner iterations
and outer iterations, respectively. For example, $\boldsymbol{B}_{m}$ is an
estimate at iteration $m$ to a problem
$\min_{\boldsymbol{B}}~{}f_{\rho(t)}(\boldsymbol{B})$ parameterized by a
penalty coefficient $\rho(t)$, whereas $\boldsymbol{B}_{t}$ denotes one of its
solutions, either exact or approximate, at level $t$.

###  1.1 Majorizing $\epsilon$-Insensitive Empirical Risk

Given a dataset with $n$ samples, $p$ predictors, and $c$ categories, the goal
in VDA is to minimize the empirical risk

$\displaystyle f(\boldsymbol{B})$ $\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}R(\boldsymbol{y}_{i},\boldsymbol{x}_{i},\boldsymbol{B}).$

Here $\boldsymbol{x}_{i}\in\mathbb{R}^{p}$ is the feature vector for sample
$i$, $\boldsymbol{y}_{i}\in\mathbb{R}^{c-1}$ is the sample’s label encoded as
a vertex of a regular simplex in $\mathbb{R}^{c-1}$, and the entries of
$\boldsymbol{B}=(b_{jk})$ are regression coefficients to be estimated. When
$c=2$, $3$, or $4$, the simplex is $[-1,1]$, an equilateral triangle, or a
regular tetrahedron centered around the origin, respectively. The predicted
value for sample $i$’s class vertex $\boldsymbol{y}_{i}$ is
$\boldsymbol{B}^{\top}\boldsymbol{x}_{i}$. We collect the rows
$\boldsymbol{x}_{i}^{\top}$ into a feature matrix $\boldsymbol{X}$ and the
rows $\boldsymbol{y}_{i}^{\top}$ into a class matrix $\boldsymbol{Y}$. To fit
$\boldsymbol{X}\boldsymbol{B}$ to $\boldsymbol{Y}$, we focus attention on the
risk function

$\displaystyle f(\boldsymbol{B})$ $\displaystyle=$
$\displaystyle\frac{1}{2n}\sum_{i=1}^{n}\max\\{0,\|\boldsymbol{y}_{i}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i}\|_{2}-\epsilon\\}^{2},$
(1)

which is the square of $\epsilon$-insensitive pseudo-distance. The residual
$\boldsymbol{y}_{i}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i}$ for sample $i$
incurs no loss when $\boldsymbol{B}^{\top}\boldsymbol{x}_{i}$ lies within a
Euclidean $\epsilon$-ball around its assigned vertex $\boldsymbol{y}_{i}$.
Otherwise, the residual contributes the reduced loss
$(\|\boldsymbol{y}_{i}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i}\|-\epsilon)^{2}$
to the risk. If $\mathcal{Y}$ denotes the collection of $c$ vertices and
$\hat{\boldsymbol{B}}$ the estimated coefficient matrix, then an unassigned
feature vector $\boldsymbol{x}$ is assigned to a class through the map

$\displaystyle\boldsymbol{x}$ $\displaystyle\mapsto$ $\displaystyle\mathop{\rm
argmin}\nolimits_{\boldsymbol{y}\in\mathcal{Y}}\|\boldsymbol{y}-\hat{\boldsymbol{B}}^{\top}\boldsymbol{x}\|,$

singling out the closest vertex. Inclusion of an intercept can be beneficial
when training class numbers are unbalanced, and easily is achieved by adding
an intercept row to $\boldsymbol{B}$ and a corresponding constant column
$\boldsymbol{1}$ to $\boldsymbol{X}$.

Squared $\epsilon$-insensitive pseudo-distances admit a quadratic surrogate
that greatly simplifies minimization of the loss $f(\boldsymbol{B})$. To
identify the surrogate, consider the majorization

$\displaystyle\inf_{\boldsymbol{v}}\\{f_{1}(\boldsymbol{v})+f_{2}(\boldsymbol{u}-\boldsymbol{v})\\}$
$\displaystyle\leq$ $\displaystyle
f_{1}(\boldsymbol{v}_{m})+f_{2}(\boldsymbol{u}-\boldsymbol{v}_{m})$

of the infimal convolution of two function $f_{1}(\cdot)$ and $f_{2}(\cdot)$.
In addition to the indicated inequality for all $\boldsymbol{u}$, majorization
requires the tangency condition of equality when $\boldsymbol{u}$ equals the
anchor point $\boldsymbol{v}_{m}$. This is achieved by choosing
$\boldsymbol{v}_{m}$ to be the minimum point of the infimal convolution
corresponding to $\boldsymbol{u}_{m}$. For the closed set
$C=[-\epsilon,\epsilon]$, the $0/\infty$ indicator
$f_{1}(\boldsymbol{v})=\delta_{C}(\|\boldsymbol{v}\|)$ and the squared
Euclidean distance $f_{2}(\boldsymbol{v})=\|\boldsymbol{v}\|^{2}$ we deduce
the majorization

$\displaystyle\max\\{0,\|\boldsymbol{u}\|-\epsilon\\}^{2}$ $\displaystyle=$
$\displaystyle\min_{\boldsymbol{v}}\delta_{[-\epsilon,\epsilon]}(\|\boldsymbol{v}\|)+\|\boldsymbol{u}-\boldsymbol{v}\|^{2}\;\>\leq\;\>\|\boldsymbol{u}-\boldsymbol{v}_{m}\|^{2};$
$\displaystyle\boldsymbol{v}_{m}$ $\displaystyle=$
$\displaystyle\begin{cases}\boldsymbol{u}_{m},&\|\boldsymbol{u}_{m}\|\leq\epsilon\\\
\frac{\epsilon}{\|\boldsymbol{u}_{m}\|}\boldsymbol{u}_{m},&\|\boldsymbol{u}_{m}\|>\epsilon\end{cases}$

relevant to VDA.

Applying this majorization to model (1) term-by-term, we arrive at the overall
majorization

$\displaystyle f(\boldsymbol{B})$ $\displaystyle\leq$
$\displaystyle\frac{1}{2n}\sum_{i=1}^{n}\|\boldsymbol{r}_{i}-\boldsymbol{r}_{mi}\|^{2},$
$\displaystyle\boldsymbol{r}_{i}$ $\displaystyle=$
$\displaystyle\boldsymbol{y}_{i}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i},$
$\displaystyle\boldsymbol{r}_{mi}$ $\displaystyle=$
$\displaystyle\begin{cases}\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i},&\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|\leq\epsilon\\\
\frac{\epsilon}{\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|}(\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}),&\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|>\epsilon.\end{cases}$

Here the estimate at iteration $m$, denoted $\boldsymbol{B}_{m}$, serves as an
anchor for the surrogate function. It is straightforward to verify that
tangency occurs when $\boldsymbol{B}=\boldsymbol{B}_{m}$. To make further
progress, consider two cases concerning the squared term. When the residual is
modest, we have

$\displaystyle\boldsymbol{r}_{i}-\boldsymbol{r}_{mi}$ $\displaystyle=$
$\displaystyle(\boldsymbol{y}_{i}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i})-(\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i})\;\>=\;\>\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i}.$

In the opposite case, the penalty argument reads

$\displaystyle\boldsymbol{r}_{i}-\boldsymbol{r}_{mi}$ $\displaystyle=$
$\displaystyle\left(1-\frac{\epsilon}{\|\boldsymbol{y}_{i}^{\top}-\boldsymbol{x}_{i}^{\top}\boldsymbol{B}_{m}\|}\right)\boldsymbol{y}_{i}+\frac{\epsilon}{\|\boldsymbol{y}_{i}^{\top}-\boldsymbol{x}_{i}^{\top}\boldsymbol{B}_{m}\|}\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i}.$

In either case, we have successfully isolated the linear predictor
$\boldsymbol{B}^{\top}\boldsymbol{x}_{i}$. This leads to the quadratic
surrogate

$\displaystyle g(\boldsymbol{B}\mid\boldsymbol{B}_{m})$ $\displaystyle=$
$\displaystyle\frac{1}{2n}\sum_{i=1}^{n}\|\boldsymbol{z}_{mi}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i}\|^{2}\;\>=\;\>\frac{1}{2n}\|\boldsymbol{Z}_{m}-\boldsymbol{X}\boldsymbol{B}\|_{F}^{2},$
(2)

where $\boldsymbol{X}\boldsymbol{B}$ predicts $\boldsymbol{Y}$ and
$\boldsymbol{Y}$ is replaced by $\boldsymbol{Z}_{m}$ with rows

$\displaystyle\boldsymbol{z}_{mi}^{\top}$ $\displaystyle=$
$\displaystyle\begin{cases}\boldsymbol{x}_{i}^{\top}\boldsymbol{B}_{m},&\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|\leq\epsilon\\\
w_{mi}\boldsymbol{y}_{i}^{\top}+(1-w_{m,i})\boldsymbol{x}_{i}^{\top}\boldsymbol{B}_{m},&\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|>\epsilon\end{cases}$
$\displaystyle w_{mi}$ $\displaystyle=$
$\displaystyle\frac{\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|-\epsilon}{\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|}.$

In summary, our surrogate can be interpreted as an ordinary sum of squares
with shifted responses. Observe that the tangency requirement
$g(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})=f(\boldsymbol{B}_{m})$ implies

$\displaystyle f(\boldsymbol{B}_{m})$ $\displaystyle=$
$\displaystyle\frac{1}{2n}\sum_{i=1}^{n}\|\boldsymbol{z}_{mi}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{2n}\sum_{\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|>\epsilon}\\!\\!w_{mi}^{2}\|\boldsymbol{y}_{i}-\boldsymbol{B}_{m}^{\top}\boldsymbol{x}_{i}\|^{2}.$

Hence at iteration $m$, only samples with residual errors beyond $\epsilon$
contribute to the loss. These samples tend to drive the subsequent update
$\boldsymbol{B}_{m+1}$.

###  1.2 Feature Selection via Sparsity Sets

Having derived a quadratic surrogate for the empirical risk (1), we now focus
on feature selection. Following our previous work on proximal distance
algorithms (Keys et al., 2019), we penalize $f(\boldsymbol{B})$ by the squared
Euclidean distance of $\boldsymbol{B}$ from the constraint set $S$. This
creates the objective

$\displaystyle f_{\rho}(\boldsymbol{B})$ $\displaystyle=$ $\displaystyle
f(\boldsymbol{B})+\frac{\rho}{2}\mathop{\rm
dist}\nolimits(\boldsymbol{B},S)^{2}$ (3)

for $\rho>0$ large. In binary classification $c=2$, so a natural choice for
$S$ is

$\displaystyle S_{k}$ $\displaystyle=$
$\displaystyle\Big{\\{}\boldsymbol{b}\in\mathbb{R}^{p}:\sum_{j=1}^{p}1_{\\{b_{j}\neq
0\\}}\leq k\Big{\\}}$

for $\boldsymbol{b}=(b_{j})$. In words, the sparsity set $S_{k}$ consists of
those vectors with at most $k$ nonzero components. When an intercept
$\boldsymbol{b}_{0}$ is included in the model, the intercept is usually
omitted in the sparsity set. Minimizing the penalized objective
$f_{\rho}(\boldsymbol{B})$ with $S=S_{k}$ delivers a solution minimizing
empirical risk with at most $k<p$ distinguishing features.

When the number of classes $c>2$, the most plausible way to define sparsity is
through the number of active features. A feature corresponds to a row
$\boldsymbol{b}_{j}$ of $\boldsymbol{B}$, so we encode sparsity through the
constraint set

$\displaystyle S_{k}^{\mathrm{row}}$ $\displaystyle=$
$\displaystyle\Bigg{\\{}\boldsymbol{B}=\begin{pmatrix}\boldsymbol{b}_{1}\\\
\cdots\\\
\boldsymbol{b}_{p}\end{pmatrix}\in\mathbb{R}^{p\times(c-1)}:\sum_{j=1}^{p}1_{\\{\|\boldsymbol{b}_{j}\|=0\\}}\leq
k\Bigg{\\}}.$

Thus, at most $k$ active rows (features) are allowed in the model. Note that
$S_{k}^{\mathrm{row}}\equiv S_{k}$ when $c=2$. Regardless of the constraint
set $S$, distance majorization replaces $\mathop{\rm
dist}\nolimits(\boldsymbol{B},S)^{2}$ with
$\|\boldsymbol{B}-P_{S}(\boldsymbol{B}_{m})\|_{F}^{2}$, where
$P_{S}(\boldsymbol{B})$ denotes the Euclidean projection of $\boldsymbol{B}$
onto $S$. For a sparsity set, projection acts by ranking features according to
their norms $\|\boldsymbol{b}_{j}\|$ and then identifying the top $k$ features
via a partial sort. Omitted rows are then set to $\boldsymbol{0}^{\top}$. The
intercept row is ignored in projection. In summary, projection onto a sparsity
set amounts to computing order statistics on estimates of model coefficients
and thresholding to select the highest ranking features.

In the special case $S=S_{k}^{\mathrm{row}}$, combining majorization (2) with
distance majorization leads to the overall quadratic surrogate

$\displaystyle g_{\rho}(\boldsymbol{B}\mid\boldsymbol{B}_{m})$
$\displaystyle=$
$\displaystyle\frac{1}{2n}\|\boldsymbol{Z}_{m}-\boldsymbol{X}\boldsymbol{B}\|_{F}^{2}+\frac{\rho}{2}\|\boldsymbol{P}_{m}-\boldsymbol{B}\|_{F}^{2}$
(4)

with the shifted response matrix $\boldsymbol{Z}_{m}$ and projection
$P_{S}(\boldsymbol{B}_{m})$. This convex quadratic serves as the main
ingredient in deriving our novel proximal distance algorithms for VDA.

###  1.3 Algorithms

Let us briefly recall the MM principle (Lange, 2016) applied to the penalized
objective (3). Given an anchor point $\boldsymbol{B}_{m}$, the surrogate (4)
majorizes (3) so that $f_{\rho}(\boldsymbol{B})\leq
g_{\rho}(\boldsymbol{B}\mid\boldsymbol{B}_{m})$ holds for every
$\boldsymbol{B}$ in the essential domain of the penalized objective. At each
iteration we select a potentially nonunique stationary point
$\boldsymbol{B}_{m+1}\in\mathop{\rm
argmin}\nolimits_{\boldsymbol{B}}\,g_{\rho}(\boldsymbol{B}\mid\boldsymbol{B}_{m})$
minimizing the surrogate (4) around the anchor point $\boldsymbol{B}_{m}$.
These observations combined with the tangency condition reveal the chain of
inequalities

$\displaystyle
f_{\rho}(\boldsymbol{B}_{m+1})\overset{\text{majorization}}{\leq}g_{\rho}(\boldsymbol{B}_{m+1}\mid\boldsymbol{B}_{m})\overset{\text{definition}}{\leq}g_{\rho}(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})\overset{\text{tangency}}{=}f_{\rho}(\boldsymbol{B}_{m}).$

Thus, minimizing surrogate (4) guarantees the descent property in minimizing
the penalized loss (3). The remainder of this section identifies specific
iteration schemes based on this useful principle. Additional details are
available in Appendix A.

####  1.3.1 Direct Minimization

The directional derivative of the surrogate (4) in the direction
$\boldsymbol{W}$ is

$\displaystyle
d_{\boldsymbol{W}}g_{\rho}(\boldsymbol{B}\mid\boldsymbol{B}_{m})$
$\displaystyle=$ $\displaystyle-n^{-1}\mathop{\rm
tr}\nolimits[(\boldsymbol{Z}_{m}-\boldsymbol{X}\boldsymbol{B})^{\top}\boldsymbol{X}\boldsymbol{W}]+\rho\mathop{\rm
tr}\nolimits[(\boldsymbol{B}-\boldsymbol{P}_{m})^{\top}\boldsymbol{W}].$

At a stationary point
$d_{\boldsymbol{W}}g_{\rho}(\boldsymbol{B}\mid\boldsymbol{B}_{m})$ vanishes
for all $\boldsymbol{W}$. It follows that

$\displaystyle\boldsymbol{0}$ $\displaystyle=$
$\displaystyle-n^{-1}\boldsymbol{X}^{\top}(\boldsymbol{Z}_{m}-\boldsymbol{X}\boldsymbol{B})+\rho(\boldsymbol{B}-\boldsymbol{P}_{m})\;\>=\;\>\nabla
g_{\rho}(\boldsymbol{B}\mid\boldsymbol{B}_{m})$ (5)

and that the $\boldsymbol{B}$ update is

$\displaystyle\boldsymbol{B}_{m+1}$ $\displaystyle=$
$\displaystyle(n^{-1}\boldsymbol{X}^{\top}\boldsymbol{X}+\rho\boldsymbol{I}_{p})^{-1}(n^{-1}\boldsymbol{X}^{\top}\boldsymbol{Z}_{m}+\rho\boldsymbol{P}_{m}).$

To fit $\boldsymbol{B}$ across a spectrum of feature subsets and values of
$\rho$, we extract the thin singular value decomposition (SVD)
$\boldsymbol{X}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{\top}$ and
invoke the Woodbury matrix identity

$\displaystyle(\boldsymbol{A}+\boldsymbol{E}\boldsymbol{C}\boldsymbol{F})^{-1}=\boldsymbol{A}^{-1}-\boldsymbol{A}^{-1}\boldsymbol{E}(\boldsymbol{C}^{-1}+\boldsymbol{F}\boldsymbol{A}^{-1}\boldsymbol{E})^{-1}\boldsymbol{F}\boldsymbol{A}^{-1}$

with $\boldsymbol{A}=\rho\boldsymbol{I}_{p}$, $\boldsymbol{E}=\boldsymbol{V}$,
$\boldsymbol{F}=\boldsymbol{V}^{\top}$, and
$\boldsymbol{C}=n^{-1}\boldsymbol{\Sigma}^{2}$. The substitutions derived in
Appendix A.1 then yield the preferred MM update

$\displaystyle\boldsymbol{B}_{m+1}$ $\displaystyle=$
$\displaystyle\boldsymbol{P}_{m}+\boldsymbol{V}n^{-1}(n^{-1}\boldsymbol{\Sigma}^{2}+\rho\boldsymbol{I}_{r})^{-1}\left[\boldsymbol{\Sigma}\boldsymbol{U}^{\top}\boldsymbol{Z}_{m}-\boldsymbol{\Sigma}^{2}\boldsymbol{V}^{\top}\boldsymbol{P}_{m}\right]$
(6)

involving the inverse of an $r\times r$ diagonal matrix, where $r$ equals the
rank of $\boldsymbol{\Sigma}$. While the required thin SVD is expensive to
extract in high dimensions, it is done only once. This strategy successfully
decouples the annealing parameter $\rho$ and model size $k$ from the multitude
of matrix operations.

####  1.3.2 Steepest Descent

Alternatively, one can take advantage of the quadratic form of (4) to derive a
step size $\gamma_{m}$ in gradient descent update
$\boldsymbol{B}_{m+1}=\boldsymbol{B}_{m}-\gamma_{m}\nabla
g_{\rho}(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})$. The Taylor expansion

$\displaystyle g_{\rho}(\boldsymbol{B}\mid\boldsymbol{B}_{m})$
$\displaystyle=$ $\displaystyle
g_{\rho}(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})+\mathop{\rm
tr}\nolimits\left\\{\nabla
g_{\rho}(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})^{\top}(\boldsymbol{B}-\boldsymbol{B}_{m})\right\\}$
$\displaystyle+$ $\displaystyle\frac{1}{2}\mathop{\rm
tr}\nolimits\left\\{(\boldsymbol{B}-\boldsymbol{B}_{m})^{\top}\nabla^{2}g_{\rho}(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})(\boldsymbol{B}-\boldsymbol{B}_{m})\right\\},$

is exact. In the particular case
$\boldsymbol{B}=\boldsymbol{B}_{m}-\gamma\nabla
g_{\rho}(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})$, elementary calculus leads
to the steepest descent update

$\displaystyle\gamma_{m}$ $\displaystyle=$
$\displaystyle\frac{\|\boldsymbol{G}_{m}\|_{F}^{2}}{n^{-1}\|\boldsymbol{X}\boldsymbol{G}_{m}\|_{F}^{2}+\rho\|\boldsymbol{G}_{m}\|_{F}^{2}}\qquad\boldsymbol{B}_{m+1}=\boldsymbol{B}_{m}+\gamma_{m}\boldsymbol{G}_{m},$
(7)

where $\boldsymbol{G}_{m}=\nabla
g_{\rho}(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})$ is defined by formula (5).

###  1.4 Practical Considerations and Implementation

We now address several issues in formulating a viable proximal distance
algorithm based on our general framework.

####  1.4.1 Initialization

The penalized loss (3) is non-convex when $S$ is a sparsity set. Thus,
solutions may be sensitive to the initial guess $\boldsymbol{B}_{0}$. The
special cases $k=p$ and $k=0$ are natural starting points for selecting a
model via cross validation. In these two cases the penalized loss is convex
but unfortunately not strongly convex. Choosing an initial point therefore
requires some care. We recommend minimizing the regularized loss

$\displaystyle h_{\lambda}(\boldsymbol{B})$ $\displaystyle=$ $\displaystyle
f(\boldsymbol{B})+\frac{\lambda}{2}\|\boldsymbol{B}\|_{F}^{2},\quad\lambda\;\>>\;\>0,$
(8)

via the surrogate
$g_{\lambda}(\boldsymbol{B}\mid\boldsymbol{B}_{m})=\frac{1}{2n}\|\boldsymbol{Z}_{m}-\boldsymbol{X}\boldsymbol{B}\|_{F}^{2}+\frac{\lambda}{2}\|\boldsymbol{B}\|_{F}^{2}$.
Strong convexity now holds, so minimization of $f_{\lambda}(\boldsymbol{B})$
furnishes a neutral starting point. Any random initial value
$\boldsymbol{B}_{0}$ should lead to the same solution. Note that in the
special case with $k=p$ active features, the distance penalty vanishes and
minimization of surrogate (4) reduces to a potentially ill-defined least
squares problem. We therefore recommend substituting the regularized surrogate
(8) for surrogate (4) to stitch the solution path from $k=p$ to $k=p-1$.

####  1.4.2 Convergence Criteria

Minimizing the penalized loss (3) with $\rho>0$ large generally slows down
convergence of iterative methods but is required to achieve optimality
(Beltrami, 1970). We surmount the poor convergence characteristics of penalty
methods by solving a sequence of subproblems $\eqref{eq:penalized}$
parameterized by an annealing schedule $\\{\rho(t)\\}_{t\geq 0}$. The idea is
to choose an annealing schedule with $\rho(0)$ small (for example,
$\rho(0)=1$) and gradually increase the penalty coefficient. Unfortunately,
the choice of annealing schedule sending $\rho\to\infty$ affects the quality
of solutions. It is therefore crucial to assess the quality of solutions in
minimizing (3) for a particular value $\rho(t)$ in the annealing schedule, the
inner iterations, and in minimizing across all subproblems, the outer
iterations. These issues merit a brief discussion of convergence criteria.

Our previous work relied on the condition

$\displaystyle|f_{\rho}(\boldsymbol{B}_{m+1})-f_{\rho}(\boldsymbol{B}_{m})|$
$\displaystyle\leq$ $\displaystyle\delta_{f}[1+f_{\rho}(\boldsymbol{B}_{m})],$

where $\delta_{f}$ is a control parameter that sets a minimum level of
progress per inner iteration on a relative scale (Keys et al., 2019; Landeros
et al., 2022). In retrospect, the gradient condition

$\displaystyle\|\nabla f_{\rho}(\boldsymbol{B}_{m})\|$ $\displaystyle\leq$
$\displaystyle\delta_{g},$

is probably a safer criterion because convergence can be quite slow for large
$\rho$ even under Nesterov acceleration. The tangency condition of the MM
principle conveniently allows one to compute the gradient $\nabla
f_{\rho}(\boldsymbol{B}_{m})$ through the equivalent expression $\nabla
g_{\rho}(\boldsymbol{B}_{m}\mid\boldsymbol{B}_{m})$. Controlling the distance
term across outer iterations $t$ yields an approximate solution
$\boldsymbol{B}_{t}=\mathop{\rm
argmin}\nolimits~{}f_{\rho(t)}(\boldsymbol{B})$. Thus, we impose conditions on
the distance penalty

$\displaystyle
q_{t}\leq\delta_{d}\quad\mathrm{or}\quad|q_{t}-q_{t-1}|\leq\delta_{q},$

where $q_{t}\equiv\mathop{\rm
dist}\nolimits(\boldsymbol{B}_{t},S_{k}^{\mathrm{row}})$ quantifies distance
to a sparsity set. Here the control parameter $\delta_{d}$ defines a
satisfactory level of closeness to $S_{k}^{\mathrm{row}}$, and the dual
criteria involving $\delta_{q}$ avoid excessive computing if progress slows
dramatically. Our choices allow one to gradually guide solution estimates
towards a constraint set while avoiding slow convergence due to large values
of $\rho$.

####  1.4.3 Proximal Distance Iteration

Algorithm 1 summarizes the general minimization procedure in proximal distance
iteration, including an early exit condition on the annealing path and
Nesterov acceleration. Warm starts are implicit in Algorithm 1 as one moves
from a $\rho(t)$-penalized subproblem to the next subproblem. The final
projection step is justified when an approximate solution $\boldsymbol{B}$ is
close to a sparsity set $S_{k}^{\mathrm{row}}$ in Euclidean distance; that is,
when feature selection finally stabilizes.

Algorithm 1 Proximal Distance Iteration and its Convergence

1:Set tolerances $\delta_{d},\delta_{q},\delta_{g}$, and annealing schedule
$\\{\rho(t)\\}_{t\geq 0}$.

2:Set target model size $k$ and sparsity level $s=1-k/p$.

3:Set maximum outer iterations $i_{\mathrm{outer}}$ and maximum inner
iterations $i_{\mathrm{inner}}$.

4:Set threshold for Nesterov acceleration, $i_{\mathrm{Nesterov}}$.

5:Initialize $\boldsymbol{B}$, and $q_{0}\leftarrow\mathop{\rm
dist}\nolimits(\boldsymbol{B},S_{k}^{\mathrm{row}})$.

6:Initialize counter for Nesterov acceleration $i\leftarrow 1$.

7:for outer iterations $t=1,2,\ldots,i_{\mathrm{outer}}$ do

8: Initialize $\boldsymbol{B}_{1}\leftarrow\boldsymbol{B}$ and
$\boldsymbol{N}_{1}\leftarrow\boldsymbol{B}$.

9: for inner iterations $m=1,2,\ldots,i_{\mathrm{inner}}$ do

10: if $\|\nabla f_{\rho(t)}(\boldsymbol{B}_{m})\|<\delta_{g}$ then
$\triangleright$ Convergence check for fixed $\rho\equiv\rho(t)$.

11: Break inner loop.

12: else

13: Solve the subproblem $\boldsymbol{B}_{m+1}\leftarrow\mathop{\rm
argmin}\nolimits g_{\rho(t)}(\boldsymbol{B}\mid\boldsymbol{N}_{m})$.
$\triangleright$ e.g. using algorithm map (6) or (7).

14: if $f_{\rho(t)}(\boldsymbol{B}_{m+1})<f_{\rho(t)}(\boldsymbol{B}_{m})$ AND
$m\geq i_{\mathrm{Nesterov}}$ then $\triangleright$ Stabililty check in
Nesterov acceleration.

15: Accelerate
$\boldsymbol{N}_{m+1}\leftarrow\boldsymbol{B}_{m+1}+\frac{i-1}{i+2}(\boldsymbol{B}_{m+1}-\boldsymbol{B}_{m})$.

16: Increment $i\leftarrow i+1$.

17: else

18: Reset Nesterov acceleration; $i\leftarrow 1$ and
$\boldsymbol{N}_{m+1}\leftarrow\boldsymbol{B}_{m+1}$.

19: end if

20: end if

21: end for

22: Update $\boldsymbol{B}\leftarrow\boldsymbol{B}_{m}$ and set
$q_{t}\leftarrow\mathop{\rm
dist}\nolimits(\boldsymbol{B},S_{k}^{\mathrm{row}})$.

23: if $q_{t}<\delta_{d}$ OR $|q_{t}-q_{t-1}|<\delta_{q}[1+q_{t-1}]$ then
$\triangleright$ Convergence check for penalized problem.

24: Break outer loop.

25: end if

26:end for

27:Project final estimate $\boldsymbol{B}\leftarrow
P_{S_{k}^{\mathrm{row}}}(\boldsymbol{B})$.

Let us now briefly address the theoretical convergence properties of Algorithm
1. The sparsity constraints described in Section 1.2 are defined by closed but
nonconvex sets. In principle, projection operators onto nonconvex can be
multivalued. For example, projection of the point $(5,2,2,3)$ onto the
sparsity set $S_{3}$ is nonunique due to ties. Fortunately, the theoretical
results from Keys et al. (2019) buttress Algorithm 1 under Zangwill’s global
convergence theorem when the penalty coefficient $\rho$ is fixed (Luenberger
and Luenberger, 1984). First, $S_{k}^{\mathrm{row}}$ is closed for any choice
of $k$ and nonempty. Second, the $\epsilon$-insensitive loss
$f(\boldsymbol{B})$ is coercive whenever the predictor matrix $\boldsymbol{X}$
has full rank. This implies the penalized loss $f_{\rho}(\boldsymbol{B})$ is
also coercive. Finally, observe that the surrogate generated by distance
majorization is related to the possibly multivalued proximal mapping (Parikh,
2014) by

$\displaystyle\bigcup_{\boldsymbol{P}_{m}\in
P_{S_{k}^{\mathrm{row}}}(\boldsymbol{B}_{m})}~{}\mathop{\rm
argmin}\nolimits_{\boldsymbol{B}}\left[f(\boldsymbol{B})+\frac{\rho}{2}\|\boldsymbol{B}-\boldsymbol{P}_{m}\|^{2}\right]=\mathop{\rm
prox}\nolimits_{\rho^{-1}f}\left[P_{S_{k}^{\mathrm{row}}}(\boldsymbol{B}_{m})\right].$

Because set projections preserve closedness, it suffices to show that the
projection operator $P_{S_{k}^{\mathrm{row}}}(\boldsymbol{B})$ is single-
valued almost everywhere to satisfy the hypothesis of Proposition 8 from Keys
et al. (2019). For fixed $\rho$ and $k$ this guarantees that all limit points
of the proximal distance iterates are stationary points,
$\\{\boldsymbol{B}:\nabla f_{\rho}(\boldsymbol{B})=\boldsymbol{0}\\}$. In
fact, convexity of the loss $f(\boldsymbol{B})$ guarantees that every
surrogate $g_{\rho}(\boldsymbol{B}\mid\boldsymbol{B}_{m})$ is $\rho$-strongly
convex. The stronger conclusion of Proposition 11 from Keys et al. (2019)
applies and allays concerns regarding selection of projected points when
$P_{S_{k}^{\mathrm{row}}}(\boldsymbol{B})$ is set-valued.

####  1.4.4 Tuning the Sparsity Level

Because a causal or optimal feature set is rarely known in advance, the
effectiveness of VDA hinges on its ability to generate candidate feature sets
and select an optimal model. Cross validation is an obvious way to evaluate
the classification error of candidate models. Before discussing implementation
specifics, we note that the sparsity set $S_{k}^{\mathrm{row}}$ functions as a
discrete hyperparameter determining model selection. Changing the size of the
target feature set $k$ fundamentally changes the nature of the optimization
problem. In contrast, perturbing the penalty constant in a lasso penalized
model merely shifts or perturbs parameter estimates along a solution path
(Hastie et al., 2001). Model size is only indirectly determined by the penalty
constant. The main benefits of the sparsity set approach are that it (a)
allows us to formulate feature selection problems directly in terms of model
size, (b) relies on simple projections rather than complex combinatorial
searches to guide solutions towards an optimal or nearly optimal sparsity
pattern, and (c) ultimately delivers models with excellent statistical
properties.

For a given dataset $\\{(\boldsymbol{y}_{i},\boldsymbol{x}_{i})\\}_{i=1}^{n}$,
we implement cross validation by partitioning the samples into two subsets,
one for $F$-fold cross validation and one for evaluating out-of-sample
prediction. The former subset is further split into a training set, used to
estimate a solution path spanning several model sizes $k$, and a validation
set, used to identify an optimal model size. Samples in the latter testing set
never appear in the cross validation procedure. Given any set $\mathcal{T}$ of
size $|\mathcal{T}|$ and model size $k$, let

$\displaystyle\mathrm{CV}_{\mathcal{T}}(\boldsymbol{B},k)$ $\displaystyle=$
$\displaystyle\frac{1}{F|\mathcal{T}|}\sum_{f=1}^{F}\sum_{(\boldsymbol{y},\boldsymbol{x})\in\mathcal{T}}1\\{\boldsymbol{y}\neq\hat{\boldsymbol{y}}_{f}\\}$

denote the classification error, on a $[0,1]$ scale, averaged over $F$ folds.
Here the $\hat{\boldsymbol{y}}_{f}$ are the predicted values under the
parameter matrix $\boldsymbol{B}$. This metric allows us to choose an optimal
$k$. The same random folds apply to all values of $k$. For each $k$ we
standardize features, find the best fitting $\boldsymbol{B}$ on the training
data, and evaluate the cross validation error on the training, validation, and
testing sets. The optimal $k$ over a grid is chosen to minimize validation
error over a predetermined grid $\mathcal{G}$ of plausible model sizes.
Because estimates of cross validation error and therefore the selected model
are sensitive to the choice of training and validation sets, we repeat the
cross validation procedure several times by shuffling samples in the cross
validation subset. These replicates allow us to investigate the stability of
the solution paths and the optimal choice of $k$. Algorithm 2 summarizes the
entire process.

Algorithm 2 Cross Validation for Sparse VDA Algorithms

1:Split dataset into testing and CV subsets.

2:for each replicate $r\leftarrow 1,R$ do

3: Shuffle CV subset.

4: for each $\mathrm{fold}\leftarrow 1,F$ do

5: Split CV subset into training and validation subsets.

6: Estimate means $\boldsymbol{\mu}$ and variances $\boldsymbol{\sigma}$ for
features using the training subset.

7: Standardize all three subsets using
$(\boldsymbol{\mu},\boldsymbol{\sigma})$.

8: Fit solution path $\\{\hat{\boldsymbol{B}}_{k}\\}_{k\in\mathcal{G}}$ to
training data over a grid $\mathcal{G}$ of model sizes.

9: Evaluate classification error on training, validation, and testing subsets.

10: end for

11: Estimate $\mathrm{CV}_{\mathrm{train}}$, $\mathrm{CV}_{\mathrm{val.}}$,
and $\mathrm{CV}_{\mathrm{test}}$ by averaging over folds.

12: Rank models by the scores
$\mathrm{CV}_{\mathrm{val.}}(\hat{\boldsymbol{B}}_{k},k)$.

13: Identify and record optimal model size $k_{\mathrm{opt}}$.

14:end for

15:for each metric $\leftarrow$ ($\mathrm{CV}_{\mathrm{train}}$,
$\mathrm{CV}_{\mathrm{val.}}$, $\mathrm{CV}_{\mathrm{test}}$,
$k_{\mathrm{opt}}$) do

16: Compute its median value over all $R$ replicates.

17: Construct a 95% credible equal-tailed interval.

18:end for

##  2 Extension to Nonlinear VDA

Linear VDA for feature selection easily generalizes to nonlinear VDA based on
reproducing kernels. Construction of decision boundaries in sparse nonlinear
VDA operates by selecting representative samples called support vectors in the
SVM literature (Wu and Wu, 2012). Here we focus on defining and solving the
optimization problem generated by nonlinear VDA. The goal is to estimate a
classification vector
$\boldsymbol{\psi}(\boldsymbol{x})=[\psi_{1}(\boldsymbol{x}),\psi_{2}(\boldsymbol{x}),\ldots,\psi_{c-1}(\boldsymbol{x})]$
by minimizing the objective
$(2n)^{-1}\sum_{i=1}^{n}\max\\{\|\boldsymbol{y}_{i}-\boldsymbol{\psi}(\boldsymbol{x}_{i})\|-\epsilon,0\\}^{2}$.
Given a reproducing kernel $W(\cdot,\cdot)$,
$\boldsymbol{\psi}(\boldsymbol{x})$ assumes the functional form

$\displaystyle\psi_{j}(\boldsymbol{x})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}b_{ij}W(\boldsymbol{x},\boldsymbol{x}_{i})+b_{0j},\qquad
j=1,2,\ldots,(c-1)$

in predicting the vertex components from a feature vector $\boldsymbol{x}$.
Model parameters $\boldsymbol{B}=(b_{ij})$, including the intercepts $b_{0j}$,
can be estimated by minimizing by criterion (3) with $\boldsymbol{W}$
substituted for $\boldsymbol{X}$. Nonlinear prediction is achieved through
through the introduction of the reproducing kernel. The underlying model is
still linear in its parameters.

##  3 Numerical Experiments

In this section we demonstrate the feature selection capabilities of our VDA
algorithms leveraging proximal distance iteration. Our first simple example
demonstrates how VDA reliably recovers sparse models without inflating (a)
misclassification errors, (b) false positives in selected features, or (c)
false negatives in discarded features. Next we report benchmark results for
both linear and nonlinear sparse VDA for datasets from the UCI Machine
Learning Repository (Dua and Graff, 2019) and Elements of Statistical Learning
(Hastie et al., 2001). We conclude with a brief application to classical
microarray expression data (Dettling and Bühlmann, 2002). Appendix B provides
detailed descriptions of simulations, datasets, data preprocessing, and cross
validation settings. Appendix C documents our choices for additional settings
in VDA and describes our computing environment.

###  3.1 Sparse Signal Recovery

Model selection is likely to succeed only when (a) classes are sufficiently
distinct, with between-class variation large compared to within-class
variation, and (b) sample size is sufficiently large to overcome noise and
masking. We now demonstrate linear VDA’s viability in these ideal
circumstances. Based on the simulation model described in Section B.5), we
generate synthetic datasets with $p=50$ features, $c=10$ classes, and $n=100$
or $n=1000$ samples equally distributed among classes. We also vary the
within-class variance $\sigma^{2}$ and indirectly modulate between-class
variance through an intermediate parameter $d$. The number of causal features
$p_{\mathrm{causal}}=10$ is kept fixed throughout each scenario. Because the
simulation procedure can generate highly correlated causal features, the
number $p_{\mathrm{causal}}$ should be viewed as an upper bound on estimated
model size.

For a given dataset, we estimate an initial solution by minimizing the
regularized model (8) with the choice $\lambda=10^{-3}$ to avoid excessive
shrinkage. This seeds VDA with a model containing all $p$ active features. In
varying $k$ from $p-1$ to $0$, linear VDA generates several candidate sparse
solutions $\boldsymbol{B}_{k}$ by minimizing the distance penalized objective
(4). Thus, feature selection proceeds from a fully dense model
$\boldsymbol{B}_{p}$ to a completely sparse $\boldsymbol{B}_{0}$. In each
case, estimating $\boldsymbol{B}_{k-1}$ uses the previous solution
$\boldsymbol{B}_{k}$ as a warm-start. This procedure is carried out once
without cross validation to illustrate how our method generates candidate
models.

Figure 1:  Sparse recovery under various conditions. The number of features,
causal features, and classes are fixed to $50$, $10$, and $10$, respectively.
Columns report classification errors, mean squared error, and ROC curves
grouped by separation level $d$ (shapes; a proxy for between-class variation)
and within-class variation $\sigma$ (colors). The true sparsity level is
highlighted by a vertical line ($p_{\mathrm{causal}}=10$, or 80% sparsity).

Our simulation study emphasizes on the following metrics:

  * •

classification accuracy on the (entire) dataset,

  * •

mean squared error in model coefficients (MSE), defined on a scale relative to
the ground truth $\boldsymbol{B}_{0}$ as
$\mathrm{MSE}(\boldsymbol{B},\boldsymbol{B}_{0})=p^{-1}\|\boldsymbol{B}_{0}\|_{F}^{-2}\|\boldsymbol{B}-\boldsymbol{B}_{0}\|_{F}^{2}$,
and

  * •

the number of positive and negative calls relative to $\boldsymbol{B}_{0}$.

Figure 1 reports our VDA findings as receiver-operator characteristic (ROC)
curves summarizing false positive and false negative rates. It is clear that a
low model size $k$ preserves a high classification accuracy (over 70% across
scenarios). At some level the target model becomes too small for accurate
classification. The trade-off between parsimony and accuracy is reflected as
an increased risk of selecting the wrong subset of features. This in turn
increases classification errors. The green curves in Figure 1 show this
tendency with $n=100$ and $\sigma=10$. When the ratio of between-class
variance to within-class variance is sufficiently high (for example, $d=3$ and
$\sigma=1$), the MSE relative to $\boldsymbol{B}_{0}$ is minimized near the
true sparsity level, 80%. This suggests that one can expect good estimates of
model coefficients provided classes are adequately captured by a linear
classifier and the dataset is sufficiently large. The situation is only
improved as sample size $n$ increases.

###  3.2 Linear Classifier Benchmarks

We now turn our attention to fitting linear classifiers to datasets from
various applications as spelled out in Appendix B. On each dataset, we first
fit an initial linear VDA without sparsity. We then apply $5$-fold cross
validation to estimate classification errors on training, validation, and test
subsets over the full range of model sizes $k\in\\{1,2,\ldots,p\\}$, moving
from dense to sparse. The entire procedure is replicated 50 times, with each
replicate shuffling the samples between training and validation subsets. In
each replicate and each fold, we select an optimal model on the basis of
minimum error on the validation subset and minimal model size, as described in
Section 1.4.4. Finally, we summarize variability in classification errors and
selected model sizes by constructing equal-tailed 95% credible intervals from
all 50 replicates. The fitted classifier includes an intercept term. The
findings in Table 1 show that our latest versionn of VDA consistently
identifies a model more parsimonious than the full model.

Table 1:  Linear classifier benchmarks. 50 replicates are summarized as 95% equal-tailed credible intervals (minimum,median, and maximum). Classification errors correspond to the optimal sparsity level, or, equivalently, model size, on the basis of minimal validation error.  | Sparsity | Train | Validation | Test
---|---|---|---|---
| Med. | Min. | Max. | Med. | Min. | Max. | Med. | Min. | Max. | Med. | Min. | Max.
iris | 50.0 | 0.0 | 69.37 | 2.92 | 2.5 | 3.33 | 3.33 | 2.5 | 5.83 | 6.67 | 5.56 | 7.78
breast | 22.22 | 0.0 | 44.44 | 2.38 | 2.24 | 2.52 | 2.93 | 2.42 | 3.44 | 3.8 | 3.65 | 4.09
vowel | 20.0 | 0.0 | 40.0 | 45.15 | 44.08 | 45.98 | 54.93 | 53.06 | 58.1 | 52.64 | 51.96 | 53.44
letter | 0.0 | 0.0 | 0.0 | 34.3 | 34.25 | 34.39 | 34.79 | 34.64 | 35.02 | 34.5 | 34.39 | 34.61
zoo | 37.5 | 6.25 | 62.5 | 0.55 | 0.0 | 1.1 | 6.58 | 4.41 | 11.1 | 6.67 | 3.33 | 10.0
lymphography | 38.89 | 5.56 | 65.42 | 11.43 | 8.2 | 14.76 | 25.71 | 20.0 | 32.38 | 24.81 | 19.55 | 31.26
waveform | 23.81 | 0.0 | 46.55 | 10.97 | 10.22 | 11.57 | 14.93 | 12.92 | 16.53 | 16.06 | 15.39 | 16.66
optdigits | 29.69 | 10.08 | 37.5 | 2.44 | 2.36 | 2.52 | 4.29 | 4.08 | 4.55 | 4.16 | 4.0 | 4.27
splice | 92.22 | 83.58 | 93.76 | 5.57 | 5.22 | 5.88 | 6.0 | 5.46 | 6.5 | 6.15 | 5.82 | 6.49
HAR | 45.54 | 18.67 | 61.06 | 1.22 | 1.16 | 1.29 | 2.19 | 2.0 | 2.35 | 5.47 | 5.26 | 5.67

Figure 2:  Cross validation errors on validation subsets for the splice and
optdigits datasets. Results are based on solution path traversal from the full
feature set (0% sparsity) to the smallest possible model (100% sparsity).
Annotations correspond to results for the optimal model (black circle). The
horizontal black bars highlight a credible interval for the optimal model
itself to illustrate stability across replicates. Headings summarize the
number of samples in the validation set, number of features, and number of
classes.

A few of the selected examples are noteworthy. Our VDA classifier clearly
fails to identify a more parsimonious model on the letter recognition example
(letters) based on the tight intervals for classification errors.
Interestingly, VDA drops at least one feature in the classic iris dataset, in
which sepal length and sepal width fail to discriminate versicolor and
virginica. It is conceivable that features from the dataset optdigits, which
summarize white pixel counts of $32\times 32$ images along $4\times 4$ blocks,
are in fact redundant. Results on the splice example are also encouraging
because using features closest to splice junctions is known to improve
classification. In this case, the 92.22% sparsity level corresponds to a model
with 15 active features compared to 180 in the dataset. Note that each feature
in this example is a binary value indicating the presence or absence of a
particular nucleotide (A,G,C,T) across 60 sites. Thus, our fitted VDA
classifiers reliably discriminate intron-exon, exon-intron, and non-splice
sequences using approximately $15/3\approx 5$ sites.

While promising, our feature selection results with VDA have an important
caveat, namely stability. Figure 2 illustrates how classification errors vary
as VDA moves towards an increasingly parsimonious model. There is a clear dip
near the selected model in the splice example. Yet it is possible that a
fitted classifier is selected because it represents the smallest model that
does not increase errors compared to models earlier in the solution path.

###  3.3 Nonlinear Classifier Benchmarks

Table 2:  Nonlinear classifier benchmarks. 50 replicates are summarized by 95% equal-tailed credible intervals.  | Sparsity | Train | Validation | Test
---|---|---|---|---
| Med. | Min. | Max. | Med. | Min. | Max. | Med. | Min. | Max. | Med. | Min. | Max.
circles | 83.0 | 5.24 | 96.05 | 23.2 | 22.15 | 24.3 | 30.0 | 26.18 | 33.02 | 33.09 | 32.2 | 34.49
clouds | 90.25 | 47.9 | 97.0 | 1.1 | 0.8 | 1.48 | 3.2 | 2.09 | 5.2 | 6.47 | 6.12 | 6.8
vowel | 16.82 | 3.08 | 40.36 | 0.0 | 0.0 | 0.0 | 6.63 | 5.11 | 9.25 | 6.84 | 6.33 | 7.32
zoo | 4.92 | 1.64 | 49.96 | 0.0 | 0.0 | 0.0 | 20.86 | 16.79 | 30.5 | 20.0 | 16.67 | 26.67
waveform | 0.0 | 0.0 | 62.26 | 0.0 | 0.0 | 0.0 | 47.87 | 45.72 | 54.39 | 48.24 | 44.7 | 50.9

We next turn to cross validation for selecting the avatars in nonlinear,
kernel-based VDA. Sparse models have fewer avatars. Table 2 summarizes our
results. Nonlinear VDA is most successful in estimating a flexible decision
boundary with few avatars in the simulated datasets clouds and circles. The
reported testing error interval for circles compares favorably to error rates
previously reported by Wu and Wu (2012), even though it is larger than the
Bayes error (30%) used in generating the data. The testing subset for dataset
vowel is by design more difficult than the training data. Our improvement in
prediction (43.07%) over the linear VDA classifier with model selection in
Table 1 (65.26%) is consistent with our previous findings (Wu and Wu, 2012).
Unfortunately, our nonlinear classifiers diverge from established results on
the waveform example, where testing errors range from 17% to 30% using a
variety of linear and nonlinear methods.

###  3.4 Cancer Benchmarks

Table 3:  Cancer benchmarks summarized by 95% equal-tailed credible intervals based on 50 replicates.  | Sparsity | Train | Validation | Test
---|---|---|---|---
| Med. | Min. | Max. | Med. | Min. | Max. | Med. | Min. | Max. | Med. | Min. | Max.
colon | 96.95 | 65.97 | 99.88 | 0.0 | 0.0 | 0.0 | 19.98 | 12.01 | 28.05 | 22.22 | 11.74 | 30.56
srbctA | 98.77 | 91.02 | 99.61 | 0.0 | 0.0 | 0.0 | 4.04 | 0.0 | 10.14 | 5.13 | 0.0 | 14.81
leukemiaA | 99.69 | 84.41 | 99.97 | 0.0 | 0.0 | 0.0 | 7.02 | 3.42 | 13.77 | 4.76 | 0.0 | 14.29
lymphomaA | 99.28 | 50.84 | 99.92 | 1.49 | 0.0 | 3.03 | 9.8 | 2.5 | 17.44 | 16.67 | 6.18 | 21.6
brain | 95.93 | 57.23 | 99.38 | 0.0 | 0.0 | 0.0 | 23.74 | 14.76 | 34.59 | 12.5 | 4.17 | 25.0
prostate | 87.91 | 34.5 | 99.42 | 0.0 | 0.0 | 0.0 | 11.02 | 6.37 | 15.78 | 5.0 | 1.67 | 7.96

We conclude our numerical experiments with applications to gene expression in
tumors, where it is common for the number of features to dominate the number
of samples, $p\gg n$. Here Table 3 reports results for linear VDA classifiers
across 6 datasets for different cancers as prepared by Dettling and Bühlmann
(2002). The leukemia, prostate, and colon datasets involve binary
classification while the rest involve multiple classes. Our errors compare
quite favorably to previous results for VDA (Wu and Lange, 2010) and
hierarchical clustering (Dettling and Bühlmann, 2002). Although our cross
validation errors are not directly comparable to the best leave-one-out errors
noted in Table 7 of Dettling and Bühlmann (2002), we succeed in selecting
models with fewer active genes then reported in Table 9. The left-skewed
distributions for optimal sparsity (model size) in Figure 3 support this
contention in spite of the wide ranges reported for brain and prostate.

Figure 3:  Distribution of optimal models identified by our cross validation
procedure for VDA, based on 50 replicates.

##  4 Discussion and Concluding Remarks

The novelty of the current version of VDA stems from the application of
$\epsilon$-insensitive loss, a new majorization, distance penalties, and
projection onto sparsity sets. Once we majorize the loss, our algorithm
reduces to the classical proximal update

$\displaystyle\mathop{\rm prox}\nolimits_{\rho^{-1}f}(\boldsymbol{P}_{m})$
$\displaystyle=$ $\displaystyle\mathop{\rm
argmin}\nolimits_{\boldsymbol{B}}\left[f(\boldsymbol{B})+\frac{\rho}{2}\|\boldsymbol{B}-\boldsymbol{P}_{m}\|^{2}\right];\quad\boldsymbol{P}_{m}\in
P_{S_{k}^{\mathrm{row}}}(\boldsymbol{B}_{m}).$

Computing this proximal update involves solving the nonlinear equations

$\displaystyle(\boldsymbol{X}^{\top}\boldsymbol{W}_{\boldsymbol{B}}\boldsymbol{X}+\rho\boldsymbol{I})\boldsymbol{B}=\boldsymbol{X}^{\top}\boldsymbol{W}_{\boldsymbol{B}}\boldsymbol{Y}+\rho\boldsymbol{P}_{m},$

where $\boldsymbol{W}_{\boldsymbol{B}}$ is a diagonal matrix depending on the
residuals $\boldsymbol{y}_{i}-\boldsymbol{B}^{\top}\boldsymbol{x}_{i}$. This
structure suggests the standard iteratively reweighted least squares (IRLS)
strategy noticed by other authors in robust classification (Groenen et al.,
2008; van den Burg and Groenen, 2016; Nguyen and McLachlan, 2017). Sequential
minimization via our surrogate (4) is effectively an IRLS algorithm for
fitting VDA classifiers.

Our work has its limitations. While our empirical results are satisfying, we
have not characterized the statistical properties of our estimators. Notably,
we still lack conditions implying the consistency of VDA (Lange and Wu, 2008;
Wu and Lange, 2010). The simulation study in Section 3.1 makes the case that
our model estimates successfully reduce mean squared error (Figure 1).
Statistical theory would be especially helpful in understanding the stability
of model selection under cross validation. Figure 3 suggests that VDA is
strongly biased towards selecting parsimonious models on datasets believed to
have a small number of informative features. Unfortunately, cross validation
replicates often yield candidate models with different numbers of active
features. Dense models can be recovered, as illustrated by the prostate
example summarized in Figure 3. Another limitation is that only one or a few
driving features may be selected out of a cluster of highly correlated
features. This tendency may be beneficial in finding causal features in
biological models. Finally, our code implementation only relies on BLAS calls
to exploit parallelism in linear algebra operations. Our algorithms may have
additional parallel structure gone unnoticed.

Our VDA benchmarks do not report computational times required for solving a
single penalized problem or running cross validation. Without supporting
theory, it seems premature to pit VDA against the wide array of readily
available software for classification. Anecdotally, cumulative VDA fitting
times (across every fold, model size, and cross validation replicate) range
from about 5 seconds for the iris, zoo, lymphography data to several minutes
for the splice, letters, colon, and lymphoma data. The most time consuming
dataset to fit was the HAR example, which required $>7$ hours to generate 50
cross validation replicates with 5 folds each, or about $8.4$ minutes per
replicate on average. These observations are encouraging, particularly in the
more interesting case of sparse linear VDA. While our model selection
technique permits flexibility in nonlinear decision boundaries, we find that
fitting nonlinear VDA classifiers is noticeably slower than linear VDA. The
additional computational cost and loss of interpretability in nonlinear VDA is
unjustified in cases where a linear classifier is strictly better (see for
example the zoo benchmarks in Tables 1 and 2). The SVD-based variant of our
software slows noticeably with a large number of samples. Additional
computational tricks, beyond parallel linear algebra operations, are needed to
scale our VDA algorithms to datasets with sample sizes $n>10^{5}$. Mature
classification packages such as LIBSVM (Chang and Lin, 2011) may provide
insight in overcoming computing limitations.

Repeated cross validation is not the only technique for evaluating candidate
models. For instance, the method of stability selection described by
Meinshausen and Bühlmann (2010) is attractive for its theoretical finite-
sample guarantees in controlling family-wise errors. Alternatively, selecting
features in groups can help overcome multicollinearity (Alexander and Lange,
2011). Ideally group selection should also control candidate model sizes.
Without going into details, it is possible to project onto sets that limit the
number of groups selected and the number of features selected per group. The
downside of doubly sparse cross validation is that involves searching over a
2-dimensional tuning-constant grid. The percentile lasso (Roberts and Nowak,
2014) is another attractive device for stabilizing model selection. Our
repeated cross validation procedure, similar to the percentile lasso, attempts
to control variability by using the median of optimal model sizes across
replicates. The distributions reported in Figure 3 suggest that adapting
feature selection to account for the variety of models generated by replicate
cross validations may, in fact, improve credible intervals. Incorporating a
Bayesian perspective into sparse VDA might be based on the the peaks in Figure
3, but again model selection is apt to become even more time consuming.

In conclusion, we have demonstrated both the prediction accuracy and the
feature selection capabilities of our modified VDA algorithms across a range
of simulated and real datasets. Our numerical examples confirm the promise of
incorporating feature selection into VDA with distance to constraint-set
penalties. Combining our estimation procedures with cross validation overcomes
the need to specify a particular model size. Replicating cross validation
further improves the reliability of feature selection with sparse VDA and
provides a straightforward mechanism for honest reporting. In the case of
kernel VDA, our proposed methods allow for the simultaneous selection of the
support vectors defining the decision boundaries separating multiple classes.
Our future work on VDA will seek to address gaps in theory, the caveats of
repeated cross validation, and the extension of VDA to clustering problems in
analogy with SVM-based clustering.

## References

  * Alexander and Lange (2011) Alexander, D. H. and K. Lange (2011).  Stability selection for genome-wide association.  Genetic Epidemiology 35(7), 722–728.
  * Alizadeh et al. (2000) Alizadeh, A. A., M. B. Eisen, R. E. Davis, C. Ma, I. S. Lossos, A. Rosenwald, J. C. Boldrick, H. Sabet, T. Tran, X. Yu, J. I. Powell, L. Yang, G. E. Marti, T. Moore, J. Hudson, L. Lu, D. B. Lewis, R. Tibshirani, G. Sherlock, W. C. Chan, T. C. Greiner, D. D. Weisenburger, J. O. Armitage, R. Warnke, R. Levy, W. Wilson, M. R. Grever, J. C. Byrd, D. Botstein, P. O. Brown, and L. M. Staudt (2000).  Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling.  Nature 403(6769), 503–511.
  * Alon et al. (1999) Alon, U., N. Barkai, D. A. Notterman, K. Gish, S. Ybarra, D. Mack, and A. J. Levine (1999).  Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays.  Proceedings of the National Academy of Sciences 96(12), 6745–6750.
  * Beltrami (1970) Beltrami, E. J. (1970).  An Algorithmic Approach to Nonlinear Analysis and Optimization.  Academic Press.
  * Breiman et al. (1984) Breiman, L., J. H. Friedman, R. A. Olshen, and C. J. Stone (1984).  Classification And Regression Trees.  Routledge.
  * Chang and Lin (2011) Chang, C.-C. and C.-J. Lin (2011).  LIBSVM: A library for support vector machines.  ACM Transactions on Intelligent Systems and Technology 2(3), 27:1–27:27.
  * Dettling and Bühlmann (2002) Dettling, M. and P. Bühlmann (2002).  Supervised clustering of genes.  Genome Biology 3(12), research0069.1.
  * Dua and Graff (2019) Dua, D. and C. Graff (2019).  UCI Machine Learning Repository.
  * Golub et al. (1999) Golub, T. R., D. K. Slonim, P. Tamayo, C. Huard, M. Gaasenbeek, J. P. Mesirov, H. Coller, M. L. Loh, J. R. Downing, M. A. Caligiuri, C. D. Bloomfield, and E. S. Lander (1999).  Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring.  Science 286(5439), 531–537.
  * Groenen et al. (2008) Groenen, P. J. F., G. Nalbantov, and J. C. Bioch (2008).  SVM-Maj: A majorization approach to linear support vector machines with different hinge errors.  Advances in Data Analysis and Classification 2(1), 17–43.
  * Hastie et al. (2001) Hastie, T., J. Friedman, and R. Tibshirani (2001).  The Elements of Statistical Learning.  Springer Series in Statistics. Springer New York.
  * Keys et al. (2019) Keys, K. L., H. Zhou, and K. Lange (2019).  Proximal Distance Algorithms: Theory and Practice.  Journal of Machine Learning Research 20(66), 1–38.
  * Khan et al. (2001) Khan, J., J. S. Wei, M. Ringnér, L. H. Saal, M. Ladanyi, F. Westermann, F. Berthold, M. Schwab, C. R. Antonescu, C. Peterson, and P. S. Meltzer (2001).  Classification and diagnostic prediction of cancers using gene expression profiling and artificial neural networks.  Nature Medicine 7(6), 673–679.
  * Landeros et al. (2022) Landeros, A., O. H. M. Padilla, H. Zhou, and K. Lange (2022).  Extensions to the Proximal Distance Method of Constrained Optimization.  arXiv:2009.00801 [cs, math].
  * Lange (2016) Lange, K. (2016).  MM Optimization Algorithms.  SIAM-Society for Industrial and Applied Mathematics.
  * Lange and Wu (2008) Lange, K. and T. T. Wu (2008).  An MM Algorithm for Multicategory Vertex Discriminant Analysis.  Journal of Computational and Graphical Statistics 17(3), 527–544.
  * Luenberger and Luenberger (1984) Luenberger, P. i. E. E. S. a. O. R. D. G. and D. G. Luenberger (1984).  Linear and Nonlinear Programming.  Addison-Wesley.
  * Meinshausen and Bühlmann (2010) Meinshausen, N. and P. Bühlmann (2010).  Stability selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72(4), 417–473.
  * Nguyen and McLachlan (2017) Nguyen, H. D. and G. J. McLachlan (2017).  Iteratively-Reweighted Least-Squares Fitting of Support Vector Machines: A Majorization–Minimization Algorithm Approach.  In Proceedings of the 2017 Future Technologies Conference, pp. 439–446. The Science and Information Organization.
  * Parikh (2014) Parikh, N. (2014).  Proximal Algorithms.  Foundations and Trends® in Optimization 1(3), 127–239.
  * Pomeroy et al. (2002) Pomeroy, S. L., P. Tamayo, M. Gaasenbeek, L. M. Sturla, M. Angelo, M. E. McLaughlin, J. Y. H. Kim, L. C. Goumnerova, P. M. Black, C. Lau, J. C. Allen, D. Zagzag, J. M. Olson, T. Curran, C. Wetmore, J. A. Biegel, T. Poggio, S. Mukherjee, R. Rifkin, A. Califano, G. Stolovitzky, D. N. Louis, J. P. Mesirov, E. S. Lander, and T. R. Golub (2002).  Prediction of central nervous system embryonal tumour outcome based on gene expression.  Nature 415(6870), 436–442.
  * Roberts and Nowak (2014) Roberts, S. and G. Nowak (2014).  Stabilizing the lasso against cross-validation variability.  Computational Statistics & Data Analysis 70, 198–211.
  * Singh et al. (2002) Singh, D., P. G. Febbo, K. Ross, D. G. Jackson, J. Manola, C. Ladd, P. Tamayo, A. A. Renshaw, A. V. D’Amico, J. P. Richie, E. S. Lander, M. Loda, P. W. Kantoff, T. R. Golub, and W. R. Sellers (2002).  Gene expression correlates of clinical prostate cancer behavior.  Cancer Cell 1(2), 203–209.
  * Smola and Schölkopf (1998) Smola, A. J. and B. Schölkopf (1998).  Learning with Kernels, Volume 4.  Citeseer.
  * van den Burg and Groenen (2016) van den Burg, G. J. J. and P. J. F. Groenen (2016).  GenSVM: A Generalized Multiclass Support Vector Machine.  Journal of Machine Learning Research 17(224), 1–42.
  * Vapnik (1999) Vapnik, V. (1999).  The Nature of Statistical Learning Theory.  Springer science & business media.
  * Wang and Shen (2007) Wang, L. and X. Shen (2007).  On L1-Norm Multiclass Support Vector Machines.  Journal of the American Statistical Association 102(478), 583–594.
  * Wu and Lange (2010) Wu, T. T. and K. Lange (2010).  Multicategory vertex discriminant analysis for high-dimensional data.  The Annals of Applied Statistics 4(4), 1698–1721.
  * Wu and Wu (2012) Wu, T. T. and Y. Wu (2012).  Nonlinear vertex discriminant analysis with reproducing kernels.  Statistical Analysis and Data Mining: The ASA Data Science Journal 5(2), 167–176.

##  Appendix A Derivations

###  A.1 Directly minimizing surrogate (4)

According to the orthogonality relation
$\boldsymbol{V}^{\top}\boldsymbol{V}=\boldsymbol{I}_{r}$ and Woodbury’s
formula, the required inverse is

$\displaystyle[n^{-1}\boldsymbol{V}\boldsymbol{\Sigma}^{2}\boldsymbol{V}^{\top}+\rho\boldsymbol{I}_{p}]^{-1}$
$\displaystyle=$
$\displaystyle\rho^{-1}\left[\boldsymbol{I}_{p}-\rho^{-1}\boldsymbol{V}(n\boldsymbol{\Sigma}^{-2}+\rho^{-1}\boldsymbol{V}^{\top}\boldsymbol{V})^{-1}\boldsymbol{V}^{\top}\right]$
$\displaystyle=$
$\displaystyle\rho^{-1}\left[\boldsymbol{I}_{p}-\rho^{-1}\boldsymbol{V}(n\boldsymbol{\Sigma}^{-2}+\rho^{-1}\boldsymbol{I}_{r})^{-1}\boldsymbol{V}^{\top}\right]$
$\displaystyle=$
$\displaystyle\rho^{-1}\left[\boldsymbol{I}_{p}-n^{-1}\boldsymbol{V}(\rho\boldsymbol{I}+n^{-1}\boldsymbol{\Sigma}^{2})^{-1}\boldsymbol{\Sigma}^{2}\boldsymbol{V}^{\top}\right].$

Multiplying the above inverse by $\rho\boldsymbol{P}_{m}$ yields

$\displaystyle\boldsymbol{P}_{m}-n^{-1}\boldsymbol{V}(n^{-1}\boldsymbol{\Sigma}^{2}+\rho\boldsymbol{I}_{r})^{-1}\boldsymbol{\Sigma}^{2}\boldsymbol{V}^{\top}\boldsymbol{P}_{m}.$

Likewise multiplying the above inverse by
$n^{-1}\boldsymbol{X}^{\top}\boldsymbol{Z}_{m}$ gives

$\displaystyle
n^{-1}\rho^{-1}\left[\boldsymbol{V}\boldsymbol{\Sigma}\boldsymbol{U}^{\top}-n^{-1}\boldsymbol{V}(n^{-1}\boldsymbol{\Sigma}^{2}+\rho\boldsymbol{I}_{r})^{-1}\boldsymbol{\Sigma}^{2}\boldsymbol{V}^{\top}\boldsymbol{V}\boldsymbol{\Sigma}\boldsymbol{U}^{\top}\right]\boldsymbol{Z}_{m}$
$\displaystyle=$ $\displaystyle
n^{-1}\rho^{-1}\boldsymbol{V}\left[\boldsymbol{I}_{r}-n^{-1}(n^{-1}\boldsymbol{\Sigma}^{2}+\rho\boldsymbol{I}_{r})^{-1}\boldsymbol{\Sigma}^{2}\right]\boldsymbol{\Sigma}\boldsymbol{U}^{\top}\boldsymbol{Z}_{m}$
$\displaystyle=$ $\displaystyle
n^{-1}\rho^{-1}\boldsymbol{V}(n^{-1}\boldsymbol{\Sigma}^{2}+\rho\boldsymbol{I}_{r})^{-1}\left[(n^{-1}\boldsymbol{\Sigma}^{2}+\rho\boldsymbol{I}_{r})-n^{-1}\boldsymbol{\Sigma}^{2}\right]\boldsymbol{\Sigma}\boldsymbol{U}^{\top}\boldsymbol{Z}_{m}$
$\displaystyle=$ $\displaystyle
n^{-1}\boldsymbol{V}(n^{-1}\boldsymbol{\Sigma}^{2}+\rho\boldsymbol{I}_{r})^{-1}\boldsymbol{\Sigma}\boldsymbol{U}^{\top}\boldsymbol{Z}_{m}.$

Summing these two results leads to the desired update

$\displaystyle\boldsymbol{B}_{m+1}$ $\displaystyle=$
$\displaystyle\boldsymbol{P}_{m}+n^{-1}\boldsymbol{V}(n^{-1}\boldsymbol{\Sigma}^{2}+\rho\boldsymbol{I}_{r})^{-1}\left[\boldsymbol{\Sigma}\boldsymbol{U}^{\top}\boldsymbol{Z}_{m}-\boldsymbol{\Sigma}^{2}\boldsymbol{V}^{\top}\boldsymbol{P}_{m}\right].$

##  Appendix B Dataset Descriptions and Provenance

###  B.1 Preprocessing

Throughout each example we standardize each feature vector
$\boldsymbol{x}_{i}$ to have mean $\boldsymbol{0}$ and unit variance by
estimating sample means
$\boldsymbol{\mu}=\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{x}_{i}$ and variances
$\boldsymbol{\sigma}^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(\boldsymbol{x}_{i}-\boldsymbol{\mu})^{2}$,
where operations are applied element-wise. We then apply a $Z$-score
transformation
$\boldsymbol{x}\mapsto(\boldsymbol{x}-\boldsymbol{\mu})/\boldsymbol{\sigma}$.
Location and scale estimates are based on a training subset only, and the
$Z$-score transformation is applied to all subsets in cross validation. In the
nonlinear case, we estimate the variance in distances between dissimilar
samples (that is, when true class labels are distinct) using the median and
set a scale parameter $\gamma=13/10\times\mathrm{median}$ for radial basis
functions (RBF) kernels.

For real-world datasets we take care to drop missing data. In addition, we
make the choice to drop any non-numeric features (e.g. categorical data) that
are not already coded. We note this in each description where applicable.

###  B.2 Gaussian Clouds (clouds)

As described by Wu and Wu (2012), we simulate six Gaussian clouds with centers
$(x,y)$ equally spaced along the unit circle in $\mathbb{R}^{2}$.

$\theta=\begin{cases}0\,\text{or}\,\pi&\text{for class 1},\\\
\frac{\pi}{3}\,\text{or}\,\frac{4\pi}{3}&\text{for class 2},\\\
\frac{2\pi}{3}\,\text{or}\,\frac{5\pi}{3}&\text{for class
3},\end{cases}\quad\text{with}\quad(x,y)=(\cos(\theta),\sin(\theta)).$

Each class has equal probability and centers are sampled uniformly conditional
on the chosen class. Points about a given cloud are simulated with variance-
covariance matrix $\sigma^{2}\boldsymbol{I}_{2\times 2}$. Thus, higher
variability increases class overlap and hence classification difficulty. We
set $\sigma=0.25$ and split the data into training, validation, and testing
subsets of sizes 200, 50, and 1000, respectively, in 5-fold cross validation.

###  B.3 Nested Circles (circles)

The recipe described by Wu and Wu (2012) samples two features from the circle
$\\{(x,y):x^{2}+y^{2}<c\\}$ and assigns them to one of $c$ classes. The
intermediate values $r=\lceil x^{2}+y^{2}\rceil$ and $r^{\prime}$ are used to
simulate class assignment via the rule

$\text{class}=\begin{cases}r&\text{with probability}\,p\\\
r^{\prime}&\text{with probability}\,1-p.\end{cases}$

Here $\lceil\cdot\rceil$ is the ceiling function, $r^{\prime}$ is any class in
$\\{1,2,\ldots,c\\}$ other than $r$ which is sampled uniformly, and $1-p$ is
the Bayes error controlling the difficulty of classification. Our benchmark
takes $c=3$ and $p=0.8$ to simulate a dataset with $1250$ samples. We split
the data into training, validation, and testing subsets of sizes 200, 50, and
1000, respectively, in 5-fold cross validation.

###  B.4 Waveform (waveform)

This dataset is based on convex combinations of triangular waveforms that was
originally presented by Breiman et al. (1984) and featured in Example 12.7.1
of The Elements of Statistical Learning (Hastie et al., 2001). We simulate
$p=21$ features across $c=3$ classes using the recipe

$\displaystyle X_{j}$ $\displaystyle=Uh_{1}(j)+(1-U)h_{2}(j)+\epsilon_{j}$
$\displaystyle\text{for class 1},$ $\displaystyle X_{j}$
$\displaystyle=Uh_{1}(j)+(1-U)h_{3}(j)+\epsilon_{j}$ $\displaystyle\text{for
class 2},\qquad j=1,2,\ldots,21$ $\displaystyle X_{j}$
$\displaystyle=Uh_{2}(j)+(1-U)h_{3}(j)+\epsilon_{j}$ $\displaystyle\text{for
class 3},$

where $U$ is a uniform deviate on $[0,1]$ and the $\epsilon_{j}$ are
independent standard normal deviates. The triangular waveforms
$h_{1}(j)=\max\\{6-|j-11|,0\\}$, $h_{2}(j)=h_{1}(j-4)$, and
$h_{3}(j)=h_{1}(j+4)$ are centered at $j=11$, $j=15$, and $j=7$, respectively.
A total of 1375 samples are simulated with training, validation, and testing
subsets of sizes 300, 75, and 1000, respectively, in $5$-fold cross
validation.

###  B.5 Sparse Recovery Model (10clouds)

This recipe is adapted from the clouds example and a recipe from Wang and Shen
(2007) (see Section 3.1). Given a target number of features $p$ and classes
$c$, we simulate features by sampling from
$\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}_{p\times p})$, where

$\displaystyle\theta_{\ell}$
$\displaystyle=\frac{\pi}{3}+\frac{2(\ell-1)\pi}{c},$ $\displaystyle a_{\ell}$
$\displaystyle=\begin{cases}\cos\theta_{\ell},&\ell~{}\text{is even}\\\
\sin\theta_{\ell},&\text{otherwise}\end{cases},$ $\displaystyle\mu_{j}$
$\displaystyle=\begin{cases}d\times\mathop{\rm
sgn}\nolimits(a_{\ell}),&\text{for}~{}j=1,2,\ldots,c\\\
0,&\text{otherwise}\end{cases},~{}\text{and}$ $\displaystyle\Sigma_{ij}$
$\displaystyle=\begin{cases}\sigma,&i=j~{}\text{and}~{}j=1,2,\ldots,c\\\
1/\sigma,&i=j~{}\text{and}~{}j>c\\\ 0,&\text{otherwise}\end{cases}.$

For each sample $\boldsymbol{x}_{i}$ we simulate the response
$\boldsymbol{y}_{i}=\boldsymbol{B}^{\top}\boldsymbol{x}_{i}$ where each row
$\boldsymbol{b}_{j}$ of $\boldsymbol{B}$ is assigned

$\boldsymbol{b}_{j}=\begin{cases}c^{-1}\mathop{\rm
sgn}\nolimits(a_{\ell})\boldsymbol{v}_{j},&j=1,2,\ldots,c\\\
\boldsymbol{0},&\text{otherwise}\end{cases}$

and each $\boldsymbol{v}_{j}$ is one of $c$ vertices assigned to classes.
These choices increase the likelihood that the simulated dataset contains
close to $c$ informative features. Large values of parameter $d$ induce
separation between classes whereas $\sigma$ controls variability within each
class.

###  B.6 UCI Datasets

We select 8 datasets from the UCI Machine Learning Repository (Dua and Graff,
2019), all of which are in the overdetermined regime. Table 4 summarizes the
characteristics of each selected example.

Table 4:  Summary of datasets from UCI Machine Learning Repository and corresponding cross validation settings. Values in parentheses indicate the number of samples, features, or classes dropped from the original dataset.  | # classes | # samples | # features | # folds | Train / Test
---|---|---|---|---|---
iris | 3 | 150 | 4 | 3 | 120 / 30
lymphography | 4 | 148 | 18 | 3 | 105 / 43
zoo | 7 | 101 | 16 (1) | 3 | 91 / 10
bcw | 2 | 699 | 9 | 5 | 562 / 137
splice | 3 | 3186 | 180 | 5 | 2549 / 637
letters | 26 | 20000 | 16 | 5 | 16000 / 4000
optdigits | 10 | 5620 | 64 | 5 | 3823 / 1797
HAR | 6 | 10299 | 561 | 5 | 7352 / 2947

###  B.7 Cancer Datasets

Our experiments on cancer microarray expression data is based on 6
preprocessed datasets as described by Dettling and Bühlmann (2002) and
previously available from
https://stat.ethz.ch/$\sim$dettling/supercluster.html. The detailed study of
Dettling and Bühlmann (2002) reports leave-one-out cross validation error
rates on leukemia (1.39%), prostate cancer (4.90%), colon cancer (16.13%),
Small-Blue-Round-Cell-Tumor (SBRCT) cancer (0.00%), lymphoma (0.00%), and
brain cancer (11.90%) using a supervised clustering method that compared
favorably against existing literature at the time. Notably, all 6 examples are
in the underdetermined regime.

Table 5:  Summary of cancer microarray expression data and corresponding cross validation settings.  | # classes | # samples | # features | # folds | Train / Test | References
---|---|---|---|---|---|---
leukemia | 2 | 72 | 3571 | 3 | 58 / 14 | Golub et al. (1999)
prostate | 2 | 102 | 6033 | 3 | 82 / 20 | Singh et al. (2002)
colon | 2 | 62 | 2000 | 3 | 50 / 12 | Alon et al. (1999)
SRBCT | 4 | 63 | 2308 | 3 | 50 / 13 | Khan et al. (2001)
lymphoma | 3 | 62 | 4026 | 3 | 50 / 12 | Pomeroy et al. (2002)
brain | 5 | 42 | 5597 | 3 | 34 / 8 | Alizadeh et al. (2000)

###  B.8 Vowel Dataset

This dataset is taken from Elements of Statistical Learning. It consists of
$n=990$ samples with $p=10$ features and $c=11$ classes. We use the original
train / test split of 528 / 462 in 5-fold cross-validation. A full original
description of the dataset is available at
https://hastie.su.domains/ElemStatLearn/.

###  B.9 Visualizations of Simulated Examples

Figure 4 illustrates the shape of our simulated examples.

Figure 4: Representative subsets used in 5-fold cross validation for each
simulated example. Here $n$, $p$, and $c$ correspond to numbers of samples,
features, and classes, respectively.

##  Appendix C Additional Details

This section addresses additional details not covered in the preceding
appendices.

###  C.1 Class Vertices and Maximal Dead Zone $\epsilon$

We prefer encoding $c$ class labels as vertices of a unit simplex with
vertices

$\displaystyle\boldsymbol{v}_{j}$
$\displaystyle=\begin{cases}(c-1)^{-1/2}\boldsymbol{1},&\text{if}~{}j=1\\\
a\boldsymbol{1}+b\boldsymbol{e}_{j-1},&\text{if}~{}2\leq j\leq c,\end{cases}$
$\displaystyle a$ $\displaystyle=-\frac{1+\sqrt{c}}{(c-1)^{3/2}}\qquad
b=\sqrt{\frac{c}{c-1}}.$

We use the maximal value of $\epsilon$ that avoids overlapping dead zones,

$\epsilon=\frac{1}{2}\sqrt{\frac{2c}{c-1}}$

###  C.2 Computing Environment

All numerical examples are run on a Linux desktop (Manjaro, 5.10.89-1)
equipped with 32 GB RAM and an Intel 10900KF processor locked at 4.9 GHz. Only
8 cores are used. Our algorithms are implemented in the Julia language,
version 1.7.1.