# Multiway sparse distance weighted discrimination

Bin Guo1, Lynn E. Eberly1,2, Pierre-Gilles Henry2,
Christophe Lenglet2, Eric F. Lock1

1 Division of Biostatistics, School of Public Health
2 Center for Magnetic Resonance Research
University of Minnesota

###### Abstract

Modern data often take the form of a multiway array. However, most
classification methods are designed for vectors, i.e., 1-way arrays. Distance
weighted discrimination (DWD) is a popular high-dimensional classification
method that has been extended to the multiway context, with dramatic
improvements in performance when data have multiway structure. However, the
previous implementation of multiway DWD was restricted to classification of
matrices, and did not account for sparsity. In this paper, we develop a
general framework for multiway classification which is applicable to any
number of dimensions and any degree of sparsity. We conducted extensive
simulation studies, showing that our model is robust to the degree of sparsity
and improves classification accuracy when the data have multiway structure.
For our motivating application, magnetic resonance spectroscopy (MRS) was used
to measure the abundance of several metabolites across multiple neurological
regions and across multiple time points in a mouse model of Friedreich’s
ataxia, yielding a four-way data array. Our method reveals a robust and
interpretable multi-region metabolomic signal that discriminates the groups of
interest. We also successfully apply our method to gene expression time course
data for multiple sclerosis treatment. An R implementation is available in the
package MultiwayClassification at
http://github.com/lockEF/MultiwayClassification.

Keywords: Distance weighted discrimination; Multiway Classification; Sparsity;
Tensors

##  1 Introduction

Development and wide deployment of advanced technologies have produced tools
that generate massive amounts of data with complex structures. These data are
often represented as a multiway or multidimensional array, which extends the
two-way data matrix to higher dimensions. This paper concerns the task of
classification from multiway data. As our motivating data application we
consider magnetic resonance spectroscopy (MRS) data for a study of
Friedreich’s Ataxia in mice. The MRS data measures the concentration of
several metabolites across multiple regions of the brain, and across multiple
time points after a treatment, yielding a four-way data array: mice $\times$
metabolites $\times$ regions $\times$ time. We are interested in identifying
signal that distinguishes the treatment groups using the totality of the MRS
data.

A naive approach is to transform the multi-dimensional array to a vector, and
then apply high-dimensional classifiers designed for vector-valued data to the
transformed vector. However, the effects of the same metabolite in different
brain regions and at different time points on the classification into the two
classes are very likely correlated. Ignoring the dependence across different
dimensions may result in inaccurate classifications and complicate
interpretation. Thus, the multiway structure should be considered in the
model. In addition, only some metabolites may be useful for distinguishing the
classes, and these distinctions may only be present for some time points.
Exclusion of uninformative features can improve classification performance and
interpretation, and this motivates an approach that also accommodates a sparse
structure. In what follows we briefly review existing high-dimensional
classifiers (Section 1.1), sparse high-dimensional classifiers (Section 1.2),
and classifiers for multiway data (Section 1.3); our methodological
contributions are summarized in Section 1.4.

###  1.1 High-dimensional classification of vectors

Traditional approaches to classification like logistic regression and Fisher’s
linear discriminant analysis (LDA) are prone to overfitting with a large
number of features, and this has motivated several classification methods for
high-dimensional vector-valued data. We can roughly divide these methods into
two categories: non-linear classifiers such as k-nearest neighbor
classification (Cover & Hart, 1967) and random forests (Breiman, 2001), and
linear classifiers such as penalized LDA (Witten & Tibshirani, 2011), support
vector machines (SVM) (Cortes & Vapnik, 1995), and distance weighted
discrimination (DWD) (Marron et al., 2007). Non-linear classifiers are
flexible to use as they require minimal assumptions, but may not have
straightforward interpretations. In contrast, linear classifiers use a
weighted sum of the measured features. SVM and DWD are both commonly used in
biomedical research. However, DWD tends to outperform SVM in scenarios with
high-dimensional data and a lower sample size; here SVM suffers from the data
piling problem, which means many cases will pile up at the discriminating
margins as a symptom of overfitting (Marron et al., 2007).

###  1.2 High-dimensional classification of vectors with sparsity

In their original formulations, SVM, DWD, and other linear classification
approaches use all available variables. However, in practice there may be only
a few important variables affecting the outcome, especially in biomedical
applications to disease classification with imaging or genomics data (Zou,
2019). Thus, methods that use all variables for classification may include too
many noise features and deteriorate the classification performance due to
error accumulation (Fan & Fan, 2008). Sparse methods, in which only a subset
of variables are used for classification, can improve performance in this
respect and also improve interpretation by identifying a small number of
informative features. Zou (2019) gave a detailed review of high-dimensional
classification methods that can account for sparsity. Many sparse classifiers
were derived based on high dimensional extensions of the classical
discriminant analysis, such as linear programming discriminant (LPD) (Cai &
Liu, 2011), penalized LDA (Witten & Tibshirani, 2011), regularized optimal
affine discriminant (ROAD) (Fan et al., 2012) and direct sparse discriminant
analysis (DSDA) (Mai et al., 2012). Other methods have extended margin-based
approaches like SVM and DWD. Hastie et al. (2009) gives a nice complete
introduction of SVM from its geometric view to its statistical loss $+$
penalty formulation. A class of methods extended SVM to enforce sparsity by
using different loss functions or penalties, such as lasso penalized SVM
(Bradley & Mangasarian, 1998), elastic-net penalized SVM (Wang et al., 2008)
and SCAD penalized SVM (Zhang et al., 2006). For DWD, the objective function
can also be represented as an analogous loss $+$ penalty formulation (Liu et
al., 2011). Wang & Zou (2016) proposed sparse DWD (SDWD) by adding lasso and
elastic-net penalties on the model coefficients, which can improve performance
and efficiency of DWD in high-dimensional classification.

###  1.3 High-dimensional classification of multiway arrays

There is a growing literature on multiway classification by extending
classifiers of vectors to multiway arrays using factorization and dimension
reduction techniques. Ye et al. (2004) and Bauckhage (2007) extended LDA and
related approaches to multiway data. Tao et al. (2005) proposed a supervised
tensor learning framework by performing a rank-1 decomposition on the
coefficients to reduce dimension, in which the coefficients are factorized
into a single set of weights for each dimension. Wimalawarne et al. (2016)
investigated tensor-based classification where a logistic loss function and a
penalty term with different continuous tensor norms for the coefficients are
considered. Pan et al. (2018) developed a classification approach that allows
for multiway data and covariate adjustment; their proposal, termed Covariate-
Adjusted Tensor Classification in High Dimensions (CATCH), assumes a multiway
structure in the residual covariance but not in the signals discriminating the
classes. Lyu et al. (2017) proposed multiway versions of DWD and SVM under the
assumption that the coefficient array is low-rank. Their implementation of
multiway DWD was shown to dramatically improve performance over two-way
classifiers when the data have multiway structure, and also tended to
outperform analogous extensions of SVM. However, their method is restricted to
use for three-way data and does not account for sparsity.

###  1.4 High-dimensional classification of multiway arrays with sparsity:
present contributions

In this paper, we developed a general framework for multiway classification
that is applicable to any number of dimensions and any degree of sparsity. Our
proposed approach builds a connection between the two approaches of multiway
DWD and sparse DWD. The central assumption is that the signal discriminating
the groups can be efficiently represented by meaningful patterns in each
dimension, which we identify by imposing a low-rank structure on the
coefficient array. Adding penalty terms to multiway DWD can enforce sparsity
and give better performance even when the data are not sparse. The rest of the
paper is organized as follows. In Section 2 we introduce our formal
mathematical framework and notation. In Section 3 we briefly review the
standard DWD, sparse DWD and multiway DWD methods. In Section 4 we describe
our proposed multiway sparse DWD method and algorithms for estimation. We show
the multiway approach has improved performance and interpretation compared
with existing methods through simulation studies (Section 5). We considered
two data applications in Section 6, and these results showed that the multiway
sparsity model is robust to data with any degree of sparsity and has
competitive classification performance over other classification techniques.
The article concludes with discussion and future directions in Section 7.

##  2 Notation and framework

Throughout this article bold lowercase characters ($\mathbf{a}$) denote
vectors, bold uppercase characters ($\mathbf{A}$) denote matrices, and
blackboard bold uppercase characters ($\mathbb{A}$) denote multiway arrays of
the specified dimension (e.g., $\mathbb{A}:P_{1}\times P_{2}\times\cdots\times
P_{K}$). Square brackets index entries within an array, e.g.,
$\mathbb{A}[p_{1},p_{2},\cdots,p_{K}]$. The operators $||\cdot||_{1}$ and
$||\cdot||_{2}$ define the generalization of the $L_{1}$ and $L_{2}$ norms,
respectively:

$||\mathbb{A}||_{1}=\sum_{p_{1}=1}^{P_{1}}\cdots\sum_{p_{K}=1}^{P_{K}}|\mathbb{A}[p_{1},p_{2},\cdots,p_{K}]|\;\;,\;\;||\mathbb{A}||_{2}^{2}=\sum_{p_{1}=1}^{P_{1}}\cdots\sum_{p_{K}=1}^{P_{K}}\mathbb{A}[p_{1},p_{2},\cdots,p_{K}]^{2}.$

The generalized inner product for two arrays $\mathbb{A}$ and $\mathbb{B}$ of
the same dimension is

$\mathbb{A}\cdot\mathbb{B}=\sum_{p_{1}=1}^{P_{1}}\cdots\sum_{p_{K}=1}^{P_{K}}\mathbb{A}[p_{1},\ldots,p_{K}]\mathbb{B}[p_{1},\ldots,p_{K}],$

the generalized outer product for $\mathbb{A}:P_{1}\times\cdots\times P_{K}$
and $\mathbb{B}:Q_{1}\times\cdots\times Q_{L}$ is
$\mathbb{A}\circ\mathbb{B}:P_{1}\times\cdots\times P_{K}\times Q_{1}\dots
Q_{L}$ where

$\mathbb{A}\circ\mathbb{B}[p_{1},\ldots,p_{K},q_{1},\ldots,q_{L}]=\mathbb{A}[p_{1},\ldots,p_{K}]\mathbb{B}[q_{1},\ldots,q_{L}].$

For our context, $\mathbb{X}$: $N\times P_{1}\times\cdots\times P_{K}$ gives
data in the form of a $K$-way array for $N$ subjects, where $P_{k}$ is the
size of the $k^{th}$ dimension for $k=1,\ldots,K$. Each subject belongs to one
of two classes denoted by $-1$ and $+1$; let $y_{i}\in\\{-1,1\\}$ give the
class label for each subject and $\mathbf{y}=[y_{1},\ldots,y_{N}]$. Our goal
is to predict the class labels $\mathbf{y}$ based on the multiway covariates
$\mathbb{X}$.

##  3 DWD and its extensions

###  3.1 Distance Weighted Discrimination (DWD)

Here we briefly describe the standard DWD for high-dimensional vector-valued
data $\mathbf{x}_{i}$ for each subject, given by the rows of
$\mathbf{X}:N\times P$. The goal of DWD and related methods is to find the
hyperplane $\mathbf{b}=[b_{1},\cdots,b_{P}]$ which best separates the two
classes via the subject _scores_ $\mathbf{X}\mathbf{b}$. To solve this
standard binary classification problem, SVM (Cortes & Vapnik, 1995) identifies
the hyperplane that maximizes the margin separating the two classes. However,
SVM can suffer from the data piling problem as shown in Marron et al. (2007),
which means many data points may pile up on the margin when the SVM is applied
to high-dimensional data. To tackle this issue, they proposed distance
weighted discrimination (DWD) which finds the separating hyperplane that
minimizes the sum of the inverse distance from the data points to the
hyperplane. The standard DWD is formulated as the following optimization
problem:

$\underset{\mathbf{r},\mathbf{b},b_{0},\bm{\xi}}{\operatorname{argmin}}\sum_{i=1}^{N}\frac{1}{r_{i}}+C\xi_{i}$

subject to $r_{i}=y_{i}(\mathbf{x}_{i}^{T}\mathbf{b}+b_{0})+\xi_{i}\geq 0$ and
$\xi_{i}\geq 0$ for $i=1,\ldots,N$, and $||\mathbf{b}||\leq 1$. Here, $C$ is
the penalty parameter and $b_{0}$ is an intercept term. The classification
rule is given by the sign of $\mathbf{x}_{i}\mathbf{b}+b_{0}$, and thus
$\xi_{i}$ can be considered a penalty for misclassification.

###  3.2 Sparse DWD

The standard DWD may not be suitable for high dimensional classification when
the underlying signal is sparse, as it does not conduct variable selection. To
further improve performance and efficiency of DWD in high-dimensional
classification, Wang & Zou (2016) proposed sparse DWD by adding penalties on
the model coefficients to enforce sparsity. The generalization of standard DWD
to sparse DWD is based on the fact that the objective function for standard
DWD can be decomposed into two components: loss function and penalty (Liu et
al., 2011):

$\underset{b_{0},\mathbf{b}}{\arg\min}\frac{1}{N}\sum_{i=1}^{N}V\left(y_{i}\left(b_{0}+\mathbf{x}_{i}^{T}\mathbf{b}\right)\right)+\frac{\lambda_{2}}{2}\|\mathbf{b}\|_{2}^{2}$
(1)

where the loss function is given by

$V(u)=\left\\{\begin{array}[]{ll}{1-u,}&{\text{ if }u\leq 1/2}\\\
{1/(4u),}&{\text{ if }u>1/2.}\end{array}\right\\}$ (2)

To account for sparsity, many variations of penalties can be added in the
model, such as lasso and elastic-net (Wang & Zou, 2016). The elastic-net
penalty often outperforms the lasso in prediction, and thus the elastic-net
penalized DWD is attractive, with the objective function

$\underset{b_{0},\mathbf{b}}{\arg\min}\frac{1}{N}\sum_{i=1}^{N}V\left(y_{i}\left(b_{0}+\mathbf{x}_{i}^{T}\mathbf{b}\right)\right)+P_{\lambda_{1},\lambda_{2}}(\mathbf{b})$
(3)

with $V(\cdot)$ defined as in (2) and

$P_{\lambda_{1},\lambda_{2}}(\mathbf{b})=\lambda_{1}||\mathbf{b}||_{1}+\frac{\lambda_{2}}{2}||\mathbf{b}||_{2}^{2}$

in which $\lambda_{1}$ and $\lambda_{2}$ are tuning parameters for
regularization. Both parameters control the shrinkage of the coefficients
toward $0$. However, the $L_{1}$ penalty controlled by $\lambda_{1}$ may
result in some model coefficients being shrunk exactly to 0, and in this way,
sparsity (retention of relatively few of the covariates for the
classification) is imposed. In practice, $\lambda_{1}$ and $\lambda_{2}$ can
be determined by cross-validation. Comparing the objectives (1) and (3), note
that the penalized elastic net DWD is equivalent to standard DWD when
$\lambda_{1}=0$.

###  3.3 Multiway DWD

The standard DWD (Section 3.1) and sparse DWD (Section 3.2) are designed for
high-dimensional classification on vector-valued data. Lyu et al. (2017)
proposed a multiway DWD model which extends the standard DWD from single
vector to multiway features. The assumption of this model is that the multiway
coefficient matrix can be decomposed into patterns that are particular to each
dimension, giving a low-rank-representation.

Consider classification of subjects with matrix valued covariates
$\mathbf{X}_{i}:P_{1}\times P_{2}$, concatenated to form the three-way array
$\mathbb{X}:N\times P_{1}\times P_{2}$ (for example, mice by metabolites by
brain regions). In this multiway context, the coefficients take the form of a
matrix $\mathbf{B}:P_{1}\times P_{2}$, and the separating hyperplane is given
by $f\left(\mathbf{X}_{i}\right)=\mathbf{X}_{i}\cdot\mathbf{B}$ where
$\mathbf{B}$ is the coefficient matrix. The rank-1 multiway DWD model assumes
that the coefficient matrix $\mathbf{B}$ has the rank-1 decomposition:
$\mathbf{B}=\mathbf{u}_{1}\mathbf{u}_{2}^{\intercal}$, where
$\mathbf{u}_{1}:P_{1}\times 1$ and $\mathbf{u}_{2}:P_{2}\times 1$ denote
vectors of weights for each dimension. Under this assumption, the hyperplane
to separate the two classes is:

$\displaystyle
f\left(\mathbf{X}_{i}\right)=\mathbf{X}_{i}\cdot\mathbf{B}=\sum_{i_{1}=1}^{P_{1}}\sum_{i_{2}=1}^{P_{2}}\mathbf{X}_{i}[i_{1},i_{2}]\mathbf{u}_{1}[i_{1}]\mathbf{u}_{2}[i_{2}].$

For example, for MRS data of the form
$\mathbf{X}_{i}:\text{Metabolites}\times\text{Regions}$, $\mathbf{u}_{1}$
gives a discriminating profile across the metabolites and $\mathbf{u}_{2}$
weights that profile across the different regions. The rank $R$ model allows
for additional patterns in each dimension via
$\mathbf{B}=\mathbf{U}_{1}\mathbf{U}_{2}^{\intercal}$ where
$\mathbf{U}_{1}:P_{1}\times R$ and $\mathbf{U}_{2}:P_{2}\times R$. The
coefficients are estimated by iteratively updating the weights in each
dimension to optimize an objective function. It has been shown that this
multiway DWD model can improve classification accuracy when the underlying
true model has a multiway structure and can provide a simple and
straightforward interpretation (Lyu et al., 2017).

##  4 Proposed methods

We propose a general framework for classifying high-dimensional multiway data
that combines aspects of sparse DWD (Section 3.2) and multiway DWD (M-DWD)
(Section 3.3). The proposed method can be considered as a multiway version of
sparse DWD that allows for any number of dimensions
$\mathbb{X}_{i}:P_{1}\times\cdots\times P_{K}$, combined across subjects to
form $\mathbb{X}:N\times P_{1}\times\cdots\times P_{K}$. In the following we
first describe our generalization of multiway DWD to an arbitrary number of
dimensions (Section 4.1), then we describe our sparsity inducing objective
function classification method that assumes the coefficient array has a low-
rank decomposition (Section 4.2).

###  4.1 Generalized multiway DWD

The generalized rank-1 multiway model assumes that the coefficient array
$\mathbb{B}:P_{1}\times\cdots\times P_{K}$ has a rank-1 decomposition
$\mathbb{B}=\mathbf{u}_{1}\circ\mathbf{u}_{2}\circ\cdots\circ\mathbf{u}_{K}$,
where $\\{\mathbf{u}_{k}\\}_{k=1}^{K}$ denote the vector of weights for each
dimension. The hyperplane for separating the two classes is
$f(\mathbb{X}_{i})=\mathbb{X}_{i}\cdot\mathbb{B}$, and thus for the rank-$1$
model

$f(\mathbb{X}_{i})=\sum_{i_{1}=1}^{P_{1}}\cdots\sum_{i_{K}=1}^{P_{K}}\mathbb{X}[i,i_{1},\cdots,i_{K}]\mathbf{u}_{1}[i_{1}]\mathbf{u}_{2}[i_{2}]\cdots\mathbf{u}_{K}[i_{K}].$

For our motivating application, $\mathbf{u}_{1}$ identifies a metabolite
profile that is weighted over the regions ($\mathbf{u}_{2}$) and time points
$(\mathbf{u}_{3})$ (or equivalently, a time profile that is weighted across
the metabolites and regions) to discriminate the classes.

The rank-1 model assumes the classification is given by combining a single
pattern in each dimension. However, in practice multiple patterns may
contribute, e.g., a different metabolite profile may affect the classification
for different regions of the brain. Thus, we propose a rank-R model for the
coefficient array $\mathbb{B}$ that is assumed to have a rank-R
Candecomp/Parafac (CP) factorization (Harshman, 1970):

$\mathbb{B}=\llbracket\mathbf{U}_{1},\dots,\mathbf{U}_{K}\rrbracket=\sum_{r=1}^{R}\mathbf{u}_{1r}\circ\cdots\circ\mathbf{u}_{Kr}$

where $\mathbf{U}_{k}:P_{k}\times R$ for $k=1,\dots,K$ with columns
$\mathbf{u}_{kr}$ as the weight for the $k^{th}$ dimension and $r^{th}$ rank
component. The coefficient array $\mathbb{B}$ in the rank-1 multiway model is
a special case of the rank-R multiway model when $R=1$. Moreover, as the CP
factorization extends the matrix rank, this model is equivalent to that for
the multiway DWD approach in Section 3.3 when $K=2$. Note that with no
constraints on the coefficient array $\mathbb{B}$ the number of free
parameters is $\prod_{k=1}^{K}P_{k}$, and with the rank-$R$ constraint
$\mathbb{B}$ has $R(P_{1}+\cdots+P_{K})$ free parameters. Thus, in addition to
facilitating interpretation of relevant patterns in each dimension, a low-rank
approach can reduce over-fitting and improve performance when there is
multiway structure.

###  4.2 Objective functions

To allow for sparsity in the generalized multiway DWD model, we consider an
extension of the sparse DWD objective (3),

$h(\mathbf{y},\mathbb{X};\mathbb{B},b_{0})=\frac{1}{N}\sum_{i=1}^{N}V\left(\mathbf{y}_{i}(b_{0}+\mathbb{X}_{i}\cdot\mathbb{B})\right)+P_{\lambda_{1},\lambda_{2}}(\mathbb{B}),$
(4)

which we minimize under the restriction that rank$(\mathbb{B})$=R. For the
penalty $P_{\lambda_{1},\lambda_{2}}(\mathbb{B})$ we consider different
extensions of the elastic net to low-rank multiway coefficients, which do or
do not distribute the $L_{1}$ and $L_{2}$ penalties across the factorization
components:

$\displaystyle P^{\mathbb{B}}_{\lambda_{1},\lambda_{2}}(\mathbb{B})$
$\displaystyle=\lambda_{1}\|\mathbb{B}\|_{1}+\frac{\lambda_{2}}{2}\|\mathbb{B}\|_{2}^{2}$
$\displaystyle P^{\mathbf{U}}_{\lambda_{1},\lambda_{2}}(\mathbb{B})$
$\displaystyle=\lambda_{1}\sum_{r=1}^{R}\prod_{k=1}^{K}\|\mathbf{u}_{kr}\|_{1}+\frac{\lambda_{2}}{2}\sum_{r=1}^{R}\prod_{k=1}^{K}\|\mathbf{u}_{kr}\|_{2}^{2},\text{
or}$ $\displaystyle P_{\lambda_{1},\lambda_{2}}(\mathbb{B})$
$\displaystyle=\lambda_{1}\sum_{r=1}^{R}\prod_{k=1}^{K}\|\mathbf{u}_{kr}\|_{1}+\frac{\lambda_{2}}{2}\|\mathbb{B}\|_{2}^{2}.$

For $R=1$, the three penalties are equivalent. For $R>1$, we find that the
solution of the objective for the penalty
$P^{\mathbf{U}}_{\lambda_{1},\lambda_{2}}(\mathbb{B})$ is often shrunk to a
lower rank than the specified rank $R$ (see Appendix B). This is due to the
distributed $L_{2}$ penalty term, and similar behavior has been observed in
other contexts; for example, when $K=2$ the distributed $L_{2}$ penalty is
equivalent to the nuclear norm penalty for matrices (Lock, 2018), which favors
a smaller rank. However, the distributed version of the $L_{1}$ penalty in
$P^{\mathbf{U}}_{\lambda_{1},\lambda_{2}}(\cdot)$ and
$P_{\lambda_{1},\lambda_{2}}(\cdot)$ is intuitive because it enforces sparsity
in the weights for each dimension (e.g., separately for the metabolites,
regions, and time points) and is amenable to iterative approaches to
optimization (see Section 4.3). Thus, in what follows, we use
$P_{\lambda_{1},\lambda_{2}}(\cdot)$ as our penalty.

This objective function for generalized sparse multiway DWD subsumes standard
DWD, multiway DWD, and sparse DWD. The objective function is equivalent to
that for standard DWD when $K=1$ and $\lambda_{1}=0$ (Marron et al., 2007; Liu
et al., 2011), it is equivalent to sparse DWD (Wang & Zou, 2016) when $K=1$
and $\lambda_{1}>0$, and it is equivalent to multiway DWD (Lyu et al., 2017)
if $\lambda_{1}=0$ and $K=2$.

###  4.3 Optimization

Here we describe the estimation algorithm for fixed penalty parameters
$\lambda_{1}$ and $\lambda_{2}$; selection of these parameters is discussed in
Section 4.4. To obtain estimates of the coefficient array $\mathbb{B}$ and
intercept $\mathbf{b}_{0}$, we iteratively optimize the objective function
$h(\mathbf{y},\mathbb{X};\mathbb{B},b_{0})$ (4) for each dimension to obtain
the estimated weights for that dimension with other dimension’s weights fixed.
This general iterative estimation approach is described in Algorithm 1.

  1. 1.

_Initialization._ Generate $K$ random matrices
$\tilde{\mathbf{U}}_{k}:P_{k}\times R$ for $k=1,...,K$. In our implementation,
entries are generated independently from a Uniform$[0,1]$ distribution.
Compute the coefficient array
$\tilde{\mathbb{B}}=\llbracket\mathbf{U}_{1},\dots,\mathbf{U}_{K}\rrbracket$,
and initialize $\tilde{b}_{0}=0$.

  2. 2.

_Iteration._ Update $\tilde{\mathbf{U}}_{1}$ and $\tilde{b}_{0}$ by optimizing
the conditional objective,

$\displaystyle\\{\tilde{b}_{0},\tilde{\mathbf{U}}_{1}\\}=\underset{\\{b_{0},\mathbf{U}_{1}\\}}{\arg\min}\;\;h(\mathbf{y},\mathbb{X};\llbracket\mathbf{U}_{1},\tilde{\mathbf{U}}_{2},\dots,\tilde{\mathbf{U}}_{K}\rrbracket,b_{0}).$
(5)

Similarly update $\tilde{\mathbf{U}}_{2},\ldots,\tilde{\mathbf{U}}_{K}$, with
the value for $b_{0}$ updated at each step. The details of the optimization
sub-step (5) are described below.

  3. 3.

_Convergence_. Set
$\mathbb{B}^{\text{New}}=\llbracket\tilde{\mathbf{U}}_{1},\dots,\tilde{\mathbf{U}}_{K}\rrbracket$.
If $||\mathbb{B}^{\text{New}}-\tilde{\mathbb{B}}||_{2}^{2}>\epsilon$ for some
pre-specified threshold $\epsilon$, update
$\tilde{\mathbb{B}}=\mathbb{B}^{\text{New}}$ and repeat Step 2. Otherwise, our
final estimates are $\hat{\mathbb{B}}=\tilde{\mathbb{B}}$ and
$\hat{b}_{0}=\tilde{b}_{0}$.

Algorithm 1 General estimation steps.

The procedure for solving the optimization problem in equation (5) depends on
whether the model is rank-$1$ $(R=1)$ or higher rank ($R\geq 2$). If $R=1$,
then $\mathbb{B}=\mathbf{u}_{1}\circ\cdots\circ\mathbf{u}_{2}$ and the
objective in (5) can be expressed as

$\frac{1}{N}\sum_{i=1}^{N}V\left(y_{i}\left(b_{0}+\tilde{\mathbf{x}}_{i}^{(1)T}\mathbf{u}_{1}\right)\right)+\left(\lambda_{1}\prod_{k=2}^{K}\|\tilde{\mathbf{u}}_{k}\|_{1}\right)\|\mathbf{u}_{1}\|_{1}+\left(\frac{\lambda_{2}}{2}\prod_{k=2}^{K}\|\tilde{\mathbf{u}}_{k}\|^{2}_{2}\right)\|\mathbf{u}_{1}\|_{2}^{2},$

where

$\displaystyle\tilde{\mathbf{x}}_{i}^{(1)}[p_{1}]=\mathbb{X}_{i}[p_{1},\dots]\cdot(\tilde{\mathbf{u}}_{2}\circ\cdots\circ\tilde{\mathbf{u}}_{K})\;\text{
for }p_{1}=1,\ldots,P_{1}.$ (6)

This is equivalent to the vector sparse DWD objective (3) as a function of
$b_{0}$ and $\mathbf{u}_{1}$, and thus can be solved using existing software
such as the sdwd function in the R package SDWD (Wang & Zou, 2015). For rank
$R\geq 2$, we solve (5) using a coordinate descent algorithm based on the
majorization-minimization (MM) principle (Hunter & Lange, 2004), thus
extending the algorithm described in Wang & Zou (2016). The steps for this
procedure are given in Algorithm 2, and analogous updates are used for
$\mathbf{U}_{2},\ldots,\mathbf{U}_{K}$.

To derive Step b of Algorithm 2, consider replacing
$\tilde{\mathbf{U}}_{1}[j,r]$ with $\mathbf{U}_{1}[j,r]$ in
$\tilde{\mathbb{B}}$ to obtain $\mathbb{B}^{*}$. Then,

$\displaystyle
h(\mathbf{y},\mathbb{X};\mathbb{B}^{*},b_{0})=\frac{1}{N}\sum_{i=1}^{N}V\left(\tilde{\mu}_{i}+y_{i}\tilde{\mathbf{X}}_{i}^{(1)}[j,r](\mathbf{U}_{1}[j,r]-\tilde{\mathbf{U}}_{1}[j,r])\right)+P_{\lambda_{1},\lambda_{2}}(\mathbb{B}^{*}),$

and we approximate $h(\mathbf{y},\mathbb{X};\mathbb{B}^{*},b_{0})$ with a
quadratic form analogous to that in Wang & Zou (2016),

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}V(\mu_{i})+\frac{1}{N}\sum_{i=1}^{N}V^{\prime}\left(\tilde{\mu}_{i}\right)(\mathbf{U}_{1}[j,r]-\tilde{\mathbf{U}}_{1}[j,r])+2(\mathbf{U}_{1}[j,r]-\tilde{\mathbf{U}}_{1}[j,r])^{2}+P_{\lambda_{1},\lambda_{2}}(\mathbb{B}^{*}).$
(7)

Note that

$P_{\lambda_{1},\lambda_{2}}(\mathbb{B}^{*})=\lambda_{1}q_{r}|\mathbf{U}_{1}[j,r]|+\frac{\lambda_{2}}{2}\left(\mathbf{W}[r,r]\mathbf{U}_{1}[j,r]^{2}+\sum_{r^{\prime}\neq
r}^{R}\mathbf{W}[r,r^{\prime}]\tilde{\mathbf{U}}_{1}[j,r^{\prime}]\right)+C,$

for a value $C$ that is constant with respect to $\mathbf{U}_{1}[j,r]$. Thus,
the minimizer of (7) over $\mathbf{U}_{1}[j,r]$ is given by
$\tilde{\mathbf{U}}_{1}[j,r]^{\text{new}}$ in Step b. A similar justification
is used for the update of the intercept in Step c.

Step a: Compute $\tilde{\mathbf{X}}_{i}^{(1)}:P_{1}\times R$ with columns
defined as in (6) for each rank component $r=1,\ldots,R$:
$\tilde{\mathbf{X}}_{i}^{(1)}[p_{1},r]=\mathbb{X}_{i}[p_{1},\dots]\cdot(\tilde{\mathbf{u}}_{2r}\circ\cdots\circ\tilde{\mathbf{u}}_{Kr})$
and let
$\tilde{\mu}_{i}=y_{i}(\tilde{b}_{0}+\tilde{\mathbf{U}}_{1}\cdot\tilde{\mathbf{X}}_{i}^{(1)})$
for $i=1,\ldots,N$ (note:
$\tilde{\mathbf{U}}_{1}\cdot\tilde{\mathbf{X}}_{i}^{(1)}=\mathbb{X}_{i}\cdot\tilde{\mathbb{B}}$).
Define $\mathbf{W}:R\times R$ by
$\mathbf{W}=(\mathbf{U}_{2}^{T}\mathbf{U}_{2})\cdots(\mathbf{U}_{K}^{T}\mathbf{U}_{K})$.

Step b: Update $\mathbf{U}_{1}$ via the MM principle and cyclic coordinate
decent over $j=1,\dots,P_{1}$ and $r=1,\dots,R$.

  * •

Compute
$z=4\tilde{\mathbf{U}}_{1}[j,r]-\frac{1}{N}\sum_{i=1}^{N}V^{\prime}(\tilde{\mu}_{i})\tilde{\mathbf{X}}_{i}^{(1)}[j,r]y_{i}-\lambda_{2}\sum_{r^{\prime}\neq
r}^{R}\mathbf{W}[r,r^{\prime}]\tilde{\mathbf{U}}_{1}[j,r^{\prime}]$, where
$V^{\prime}$ is the derivative of $V$ in (2).

  * •

Compute
$\tilde{\mathbf{U}}_{1}[j,r]^{\text{new}}=\frac{S(z,\lambda_{1}q_{r})}{4+\lambda_{2}\mathbf{W}[r,r]}$,
where $q_{r}=\prod_{k=2}^{K}\|\tilde{U}_{kr}\|_{1}$ is the weight for the
$L_{1}$ penalty, and $S(z,g)=sign(z)(|z|-g)_{+}$ is the soft-thresholding
operator in which $\omega_{+}=max(\omega,0)$.

  * •

Update
$\tilde{\mu}_{i}=\tilde{\mu_{i}}+y_{i}(\tilde{\mathbf{U}}_{1}[j.r]^{\text{new}}-\tilde{\mathbf{U}}_{1}[j,r])\tilde{\mathbf{X}}_{i}[j,r]$
for $i=1,\ldots,N$.

  * •

Set $\tilde{\mathbf{U}}_{1}[j,r]=\tilde{\mathbf{U}}_{1}[j,r]^{\text{new}}$.

Step c: Update the intercept.

  * •

Compute
$\tilde{b}_{0}^{\text{new}}=\tilde{b}_{0}-\sum_{i=1}^{N}V^{\prime}(\tilde{\mu}_{i})y_{i}/(4N)$

  * •

Update
$\tilde{\mu}_{i}=\tilde{\mu}_{i}+y_{i}(\tilde{b}_{0}^{\text{new}}-\tilde{b}_{0})$

  * •

Set $\tilde{b}_{0}=\tilde{b}_{0}^{\text{new}}$.

Step d: Repeat steps b-c until the difference between new estimates and
previous estimates is smaller than a prespecified threshold, for example,
$\sum_{j=1}^{P_{1}}\sum_{r=1}^{R}(\tilde{\mathbf{U}}_{1}[j,r]^{\text{new}}-\tilde{\mathbf{U}}_{1}[j,r])^{2}+(\tilde{b}_{0}^{\text{new}}-\tilde{b}_{0})^{2}<\epsilon$.
At the end of iterations, we have the updated weight matrix $\mathbf{U}_{1}$.

Algorithm 2 Update for $\mathbf{U}_{1}$ when $R\geq 2$.

In general, Algorithm 1 does not guarantee convergence to a global optimum,
and it can converge to local optima. We conducted simulation studies to assess
the convergence of the proposed method and other alternative methods in
Section 5.2, and in Appendix D. In practice, we find that starting with
multiple initial values and pruning the paths with inferior objective values
can improve the algorithm and obtain more robust results.

###  4.4 Selection of tuning parameters

To determine the best pair of $\lambda_{1}$ and $\lambda_{2}$, K-fold cross-
validation can be used across a grid of $\lambda_{1}$ and $\lambda_{2}$
values. For example, by default we set the following candidates for the tuning
parameters:
$\lambda_{1}=(10^{-4},0.001,0.005,0.01,0.025,0.05,0.1,0.25,0.5,0.75,1)$ and
$\lambda_{2}=(0.25,0.50,0.75,1.00,3,5)$. The misclassification rate under
cross-validation is one potential criterion; however, different parameter
values may yield the same misclassification rate, especially if the classes
are perfectly separated or nearly perfectly separated (see Section S7). We
thus recommend using t-test statistics for the test-data DWD scores between
the two classes as a more general measure of separation for selecting tuning
parameters. For each pair of $\lambda_{1}$ and $\lambda_{2}$ we compute the
predicted DWD scores for all subjects via K-fold cross-validation, and then
report t-test statistics for the difference of the predicted DWD scores
between the two classes using the test-data. We selected as optimal the pair
of $\lambda_{1}$ and $\lambda_{2}$ with the maximum t-test statistic. For a
fixed $\lambda_{2}$, the warm-start trick is used in the algorithm to select
the optimal $\lambda_{1}$. That is, we use the solution at a smaller
$\lambda_{1}$ as the initial value (the warm-start) to compute the solution at
the next $\lambda_{1}$, improving computational efficiency. The performance of
cross-validation for selecting parameters is evaluated through simulation
studies presented in Appendix A.

##  5 Simulations

In our simulation studies, we first compare the proposed rank-1 multiway
sparse model with the existing methods in Section 5.1. In Section 5.2 we
explore convergence and its impact on performance. In Section 5.3, we show the
rank-R sparse model performs well when the true data generating model has a
higher rank. Additional simulation studies are presented in the appendices,
including simulations to assess cross-validation for parameter selection, rank
mispecification, correlated predictors, and different penalization approaches.

###  5.1 Rank-1 model simulation design

To evaluate the performance of the rank-1 multiway sparse DWD model (M-SDWD),
we compared it with a rank-1 multiway sparse DWD model with $\lambda_{1}=0$
(M-SDWD $\lambda_{1}=0$), non-sparse multiway DWD (M-DWD), and the full model
sparse DWD (Full SDWD). For M-DWD we use an extended version of multiway DWD
(Lyu et al., 2017) that allows for data of any dimension $K>2$. The full model
sparse DWD is a naive way to analyze multiway data by applying the standard
sparse DWD (Wang & Zou, 2016) to the vectorized multiway data, without
considering multiway structure (i.e., no rank constraint). In this simulation
study, data were generated under several conditions, including different
multiway array dimensions, sample sizes, and sparsity levels. For all
scenarios, training datasets with two classes of equal size
$(N_{0}=N_{1}=N/2,N=40\;\text{or}\;100)$ were generated. The predictors have
the form of a three-way array of dimensions $\mathbb{X}_{i}:P_{1}\times
P_{2}\times P_{3}$. We consider two settings of different dimensionality:
higher dimensional ($30\times 15\times 15$) and lower dimensional ($15\times
4\times 5$). In each training dataset, for the $N_{0}$ samples corresponding
to class -$1$, the entries of $\mathbb{X}_{i}$ were generated independently
from a $N(0,1)$ distribution. For the other $N_{1}$ samples corresponding to
class $1$, the entries of $\mathbb{X}_{i}$ were generated independently from a
normal distribution with variance $1$ and the mean for each entry given by the
array $\sqrt{\alpha}\mu_{1}$ where
$\mu_{1}=\mathbf{u}_{1}\circ\mathbf{u}_{2}\circ\mathbf{u}_{3}$; the non-zero
values of $\mathbf{u}_{k},k=1,2,3,$ were generated independently from a
$N(0,1)$ distribution. Here, $\alpha$ gives the signal to noise ratio, i.e.,
the variance of the mean difference between the two classes over the residual
variance. We set $\alpha=0.2$ for this simulation study. The appendices
present results with different signal to noise ratios (Appendix C) and with
correlated predictors (Appendix F).

We consider three scenarios with varying degrees of sparsity: more sparsity,
less sparsity and no sparsity. Under the more sparsity scenario, about $1/3$
of the mean weights in each dimension
($\\{\mathbf{u}_{1},\mathbf{u}_{2},\mathbf{u}_{3}\\}$) are set to zero.
Specifically, under the high dimensional case ($30\times 15\times 15$), 25
entries of $\mathbf{u}_{1}$, and 10 entries each of $\mathbf{u}_{2}$ and
$\mathbf{u}_{3}$ are set to zero, that is, $5\times 5\times 5$ of the
variables have signal discriminating the classes; under the low dimensional
case ($15\times 4\times 5$), only 5 entries of $\mathbf{u}_{1}$, and 2 entries
each of $\mathbf{u}_{2}$ and $\mathbf{u}_{3}$ are non-zero. Under the less
sparsity scenario, sparsity is considered for only one dimension. In
particular, 10 entries of $\mathbf{u}_{1}$ under both the high and low
dimensional cases are set to zero, and the rest of the $\mathbf{u}_{k}$ are
non-zero. Under the no sparsity case, all variables for each $\mathbf{u}_{k}$
are nonzero. For the high dimensional case we also considered an additional
case with even more sparsity in the model: only 3 variables of each dimension
have signals and the rest of variables for each dimension are zero. Each
scenario was replicated 200 times.

Under all these scenarios, we computed the correlation of the estimated
hyperplane and the true hyperplane. Here the “true” hyperplane is the mean
difference between the classes $\mu_{1}$, which corresponds to the Bayes
linear classifier. We assess predictive performance by considering
misclassification rates for test data that were generated from the same
distributions as the training data with the same sample sizes ($N=40$ or
$100$). To assess recovery of the sparsity structure, we computed the true
positive rates and true negative rates for the proportions of non-zero/zero
weights that were correctly estimated. All these statistics were computed
based on the vectorized weights for multi-way based methods. The margin of
error across the replicates for each statistic was also computed.

Table 1 and Table 2 summarize the simulation results for lower and higher
dimensional data, respectively. Both tables show the proposed rank-1 multiway
sparse DWD model has the best performance, with higher correlation and lower
misclassification rates than other methods, when the true model is more
sparse. When the true model is less sparse or not sparse, the proposed method
has comparable performances to multiway DWD (M-DWD). All multiway based
methods have higher correlations and lower misclassification rates than the
full sparse DWD model under different scenarios. Overall the proposed rank-1
multiway sparse DWD model performs well, but some of the advantages are not
obvious. In particular, the mean correlations are relatively low for the very
sparse models, since the correlations with the truth are very small for some
simulation replications. This is because the algorithm may converge to local
minima, which is explored further in Section 5.2.

Table 1: Simulation results under the low dimensional scenario ($15\times
4\times 5$). In the Sparsity column, the numbers in parentheses indicate the
number of non-zero variables in each dimension. “Cor” is the correlation
between the estimated linear hyperplane and the true hyperplane. “Mis” is the
average misclassification rate. “TP” is the true positive rate, i.e., the
proportion of non-zero coefficients that are correctly estimated to be non-
zero.“TN” is the true negative rate, i.e., the proportion of zero coefficients
that are correctly estimated to be zero. The margins of error (2* standard
deviations across 200 replicates) for each statistic are also listed following
the $\pm$ symbol.

N | Sparsity | Methods | Cor | Mis | TP | TN
---|---|---|---|---|---|---
40 | More ($5\times 2\times 2$) | M-SDWD | 0.359$\pm$0.054 | 0.390$\pm$0.025 | 0.492$\pm$0.057 | 0.644$\pm$0.056
|  | M-SDWD ($\lambda_{1}=0$) | 0.318$\pm$0.049 | 0.386$\pm$0.025 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.313$\pm$0.049 | 0.388$\pm$0.026 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.284$\pm$0.039 | 0.402$\pm$0.023 | 0.337$\pm$0.041 | 0.748$\pm$0.043
| Less ($5\times 4\times 5$) | M-SDWD | 0.773$\pm$0.034 | 0.156$\pm$0.024 | 0.646$\pm$0.045 | 0.527$\pm$0.053
|  | M-SDWD ($\lambda_{1}=0$) | 0.768$\pm$0.035 | 0.146$\pm$0.022 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.757$\pm$0.042 | 0.154$\pm$0.024 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.532$\pm$0.031 | 0.214$\pm$0.024 | 0.320$\pm$0.038 | 0.776$\pm$0.039
| No ($15\times 4\times 5$) | M-SDWD | 0.894$\pm$0.025 | 0.052$\pm$0.015 | 0.792$\pm$0.038 | -
|  | M-SDWD ($\lambda_{1}=0$) | 0.901$\pm$0.026 | 0.048$\pm$0.014 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.896$\pm$0.030 | 0.054$\pm$0.017 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.639$\pm$0.026 | 0.096$\pm$0.020 | 0.444$\pm$0.045 | -
100 | More ($5\times 2\times 2$) | M-SDWD | 0.574$\pm$0.056 | 0.332$\pm$0.025 | 0.536$\pm$0.055 | 0.713$\pm$0.051
|  | M-SDWD ($\lambda_{1}=0$) | 0.489$\pm$0.055 | 0.329$\pm$0.025 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.489$\pm$0.056 | 0.331$\pm$0.025 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.448$\pm$0.049 | 0.337$\pm$0.024 | 0.318$\pm$0.037 | 0.841$\pm$0.034
| Less ($5\times 4\times 5$) | M-SDWD | 0.902$\pm$0.022 | 0.112$\pm$0.018 | 0.720$\pm$0.038 | 0.551$\pm$0.050
|  | M-SDWD ($\lambda_{1}=0$) | 0.890$\pm$0.025 | 0.113$\pm$0.018 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.898$\pm$0.027 | 0.111$\pm$0.018 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.695$\pm$0.028 | 0.156$\pm$0.021 | 0.357$\pm$0.034 | 0.794$\pm$0.033
| No ($15\times 4\times 5$) | M-SDWD | 0.954$\pm$0.014 | 0.035$\pm$0.012 | 0.870$\pm$0.031 | -
|  | M-SDWD ($\lambda_{1}=0$) | 0.963$\pm$0.012 | 0.035$\pm$0.011 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.965$\pm$0.014 | 0.035$\pm$0.012 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.782$\pm$0.021 | 0.061$\pm$0.015 | 0.433$\pm$0.036 | -

Table 2: Simulation results under the high dimensional scenario ($30\times
15\times 15$). In the Sparsity column, the numbers in parentheses indicate the
number of non-zero variables in each dimension. “Cor” is the correlation
between the estimated linear hyperplane and the true hyperplane. “Mis” is the
average misclassification rate. “TP” is the true positive rate, i.e., the
proportion of non-zero coefficients that are correctly estimated to be non-
zero.“TN” is the true negative rate, i.e., the proportion of zero coefficients
that are correctly estimated to be zero. The margins of error (2* standard
deviations across 200 replicates) for each statistic are also listed following
the $\pm$ symbol.

N | Sparsity | Methods | Cor | Mis | TP | TN
---|---|---|---|---|---|---
40 | Even more ($3\times 3\times 3$) | M-SDWD | 0.214$\pm$0.053 | 0.426$\pm$0.025 | 0.077$\pm$0.012 | 0.738$\pm$0.054
|  | M-SDWD ($\lambda_{1}=0$) | 0.154$\pm$0.042 | 0.415$\pm$0.025 | 0.216$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.171$\pm$0.043 | 0.404$\pm$0.025 | 0.216$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.233$\pm$0.041 | 0.388$\pm$0.024 | 0.039$\pm$0.008 | 0.900$\pm$0.035
| More ($5\times 5\times 5$) | M-SDWD | 0.644$\pm$0.056 | 0.218$\pm$0.033 | 0.489$\pm$0.052 | 0.778$\pm$0.047
|  | M-SDWD ($\lambda_{1}=0$) | 0.516$\pm$0.056 | 0.216$\pm$0.031 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.505$\pm$0.058 | 0.236$\pm$0.033 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.389$\pm$0.039 | 0.238$\pm$0.027 | 0.129$\pm$0.026 | 0.932$\pm$0.026
| Less ($10\times 15\times 15$) | M-SDWD | 0.983$\pm$0.002 | 0.000$\pm$0.000 | 0.872$\pm$0.026 | 0.394$\pm$0.060
|  | M-SDWD ($\lambda_{1}=0$) | 0.988$\pm$0.001 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.990$\pm$0.001 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.654$\pm$0.016 | 0.002$\pm$0.002 | 0.143$\pm$0.019 | 0.943$\pm$0.017
| No ($30\times 15\times 15$) | M-SDWD | 0.987$\pm$0.003 | 0.000$\pm$0.000 | 0.869$\pm$0.028 | -
|  | M-SDWD ($\lambda_{1}=0$) | 0.994$\pm$0.001 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.997$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.724$\pm$0.010 | 0.000$\pm$0.000 | 0.262$\pm$0.032 | -
100 | Even more ($3\times 3\times 3$) | M-SDWD | 0.140$\pm$0.045 | 0.467$\pm$0.016 | 0.018$\pm$0.004 | 0.784$\pm$0.049
|  | M-SDWD ($\lambda_{1}=0$) | 0.084$\pm$0.032 | 0.464$\pm$0.016 | 0.064$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.101$\pm$0.035 | 0.455$\pm$0.017 | 0.064$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.218$\pm$0.047 | 0.431$\pm$0.020 | 0.018$\pm$0.003 | 0.868$\pm$0.038
| More ($5\times 5\times 5$) | M-SDWD | 0.849$\pm$0.038 | 0.089$\pm$0.020 | 0.668$\pm$0.040 | 0.806$\pm$0.041
|  | M-SDWD ($\lambda_{1}=0$) | 0.796$\pm$0.040 | 0.101$\pm$0.021 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.766$\pm$0.049 | 0.121$\pm$0.025 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.636$\pm$0.035 | 0.136$\pm$0.022 | 0.172$\pm$0.023 | 0.957$\pm$0.022
| Less ($10\times 15\times 15$) | M-SDWD | 0.992$\pm$0.001 | 0.000$\pm$0.000 | 0.890$\pm$0.021 | 0.461$\pm$0.060
|  | M-SDWD ($\lambda_{1}=0$) | 0.995$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.996$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.816$\pm$0.009 | 0.000$\pm$0.001 | 0.199$\pm$0.014 | 0.954$\pm$0.008
| No ($30\times 15\times 15$) | M-SDWD | 0.995$\pm$0.001 | 0.000$\pm$0.000 | 0.934$\pm$0.016 | -
|  | M-SDWD ($\lambda_{1}=0$) | 0.997$\pm$0.001 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.999$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.848$\pm$0.005 | 0.000$\pm$0.000 | 0.342$\pm$0.027 | -

###  5.2 Assessment of convergence

To assess the convergence of the four methods, we show the distributions of
correlations between the true and estimated hyperplane for the “more sparsity”
scenario for high dimensional data in Figure 1. The distribution is bimodal
for the three multiway methods, where for some replications the signal is
estimated very well (correlation$\approx$1) and others do not recover the true
signal at all (correlation$\approx$0). We find that the poor performing
replications are due to the algorithm converging to a local optimum, which
does not occur for the full (vectoriozed) SDWD method. The proposed multiway
sparse DWD (M-SDWD) performs better than other methods with more reasonable
correlations. In addition, the bimodality becomes less severe when the signal
to noise ratio increases (see Appendix D). Considering prediction performance
without separating out the replications that do not converge appropriately may
complicate interpretation of the results, hense we present another two tables
that show those simulations with correlation greater than 0.5 for each method,
and the statistics among those simulations with correlation greater than $0.5$
for that method. From Table 3 and Table 4 we can see the proposed method
performs the best under the “even more” and “more” sparsity cases. To reduce
convergence to local optima, we recommend running the algorithm with multiple
initial values for the first several iterations and selecting the optimal path
to proceed.

Figure 1: Histogram of correlations between true hyperplane and estimates with
four classification methods under the more sparsity case based on high
dimensional data $(30\times 15\times 15)$ and sample size $N=100$. Table 3:
Simulation results among simulations with correlation greater than 0.5 under
the high dimensional scenario ($30\times 15\times 15$). In the Sparsity
column, the numbers in parentheses indicate the number of non-zero variables
in each dimension. “Cor” is the correlation between the estimated linear
hyperplane and the true hyperplane. “Mis” is the average misclassification
rate. “TP” is the true positive rate, i.e., the proportion of non-zero
coefficients that are correctly estimated to be non-zero.“TN” is the true
negative rate, i.e., the proportion of zero variables that are correctly
estimated to be zero. The margins of error (2* standard deviations across 200
replicates) for each statistic are also listed following the $\pm$ symbol.

N | Sparsity | Methods | Cor | Mis | TP | TN | Prop. cor$>$0.5
---|---|---|---|---|---|---|---
40 | Even more ($3\times 3\times 3$) | M-SDWD | 0.905$\pm$0.031 | 0.045$\pm$0.018 | 0.145$\pm$0.020 | 0.853$\pm$0.109 | 0.16
|  | M-SDWD ($\lambda_{1}=0$) | 0.783$\pm$0.039 | 0.074$\pm$0.029 | 0.216$\pm$0.000 | 0.000$\pm$0.000 | 0.175
|  | M-DWD | 0.816$\pm$0.035 | 0.068$\pm$0.023 | 0.216$\pm$0.000 | 0.000$\pm$0.000 | 0.175
|  | Full SDWD | 0.739$\pm$0.031 | 0.112$\pm$0.031 | 0.033$\pm$0.006 | 1.014$\pm$0.001 | 0.22
| More ($5\times 5\times 5$) | M-SDWD | 0.888$\pm$0.017 | 0.034$\pm$0.011 | 0.638$\pm$0.052 | 0.802$\pm$0.055 | 0.61
|  | M-SDWD ($\lambda_{1}=0$) | 0.840$\pm$0.022 | 0.039$\pm$0.012 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.585
|  | M-DWD | 0.863$\pm$0.021 | 0.039$\pm$0.013 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.56
|  | Full SDWD | 0.682$\pm$0.021 | 0.044$\pm$0.014 | 0.091$\pm$0.015 | 0.997$\pm$0.001 | 0.41
| Less ($10\times 15\times 15$) | M-SDWD | 0.983$\pm$0.002 | 0.000$\pm$0.000 | 0.872$\pm$0.026 | 0.394$\pm$0.060 | 1
|  | M-SDWD ($\lambda_{1}=0$) | 0.988$\pm$0.001 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 1
|  | M-DWD | 0.990$\pm$0.001 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 1
|  | Full SDWD | 0.682$\pm$0.013 | 0.000$\pm$0.000 | 0.135$\pm$0.015 | 0.956$\pm$0.011 | 0.885
| No ($30\times 15\times 15$) | M-SDWD | 0.987$\pm$0.003 | 0.000$\pm$0.000 | 0.869$\pm$0.028 | - | 1
|  | M-SDWD ($\lambda_{1}=0$) | 0.994$\pm$0.001 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | - | 1
|  | M-DWD | 0.997$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | - | 1
|  | Full SDWD | 0.727$\pm$0.009 | 0.000$\pm$0.000 | 0.262$\pm$0.032 | - | 0.99
100 | Even more ($3\times 3\times 3$) | M-SDWD | 0.922$\pm$0.051 | 0.136$\pm$0.056 | 0.045$\pm$0.009 | 1.000$\pm$0.018 | 0.09
|  | M-SDWD ($\lambda_{1}=0$) | 0.794$\pm$0.064 | 0.146$\pm$0.056 | 0.064$\pm$0.000 | 0.000$\pm$0.000 | 0.085
|  | M-DWD | 0.784$\pm$0.060 | 0.153$\pm$0.042 | 0.064$\pm$0.000 | 0.000$\pm$0.000 | 0.115
|  | Full SDWD | 0.840$\pm$0.037 | 0.194$\pm$0.033 | 0.022$\pm$0.004 | 1.015$\pm$0.001 | 0.215
| More ($5\times 5\times 5$) | M-SDWD | 0.935$\pm$0.010 | 0.048$\pm$0.010 | 0.697$\pm$0.038 | 0.831$\pm$0.042 | 0.905
|  | M-SDWD ($\lambda_{1}=0$) | 0.899$\pm$0.012 | 0.051$\pm$0.010 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.88
|  | M-DWD | 0.926$\pm$0.010 | 0.042$\pm$0.009 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.825
|  | Full SDWD | 0.760$\pm$0.017 | 0.060$\pm$0.011 | 0.156$\pm$0.016 | 0.994$\pm$0.001 | 0.765
| Less ($10\times 15\times 15$) | M-SDWD | 0.992$\pm$0.001 | 0.000$\pm$0.000 | 0.890$\pm$0.021 | 0.461$\pm$0.060 | 1
|  | M-SDWD ($\lambda_{1}=0$) | 0.995$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 1
|  | M-DWD | 0.996$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 1
|  | Full SDWD | 0.818$\pm$0.008 | 0.000$\pm$0.000 | 0.198$\pm$0.013 | 0.955$\pm$0.007 | 0.995
| No ($30\times 15\times 15$) | M-SDWD | 0.995$\pm$0.001 | 0.000$\pm$0.000 | 0.934$\pm$0.016 | - | 1
|  | M-SDWD ($\lambda_{1}=0$) | 0.997$\pm$0.001 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | - | 1
|  | M-DWD | 0.999$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | - | 1
|  | Full SDWD | 0.848$\pm$0.005 | 0.000$\pm$0.000 | 0.342$\pm$0.027 | - | 1

Table 4: Simulation results among simulations with correlation greater than
0.5 under the low dimensional scenario ($15\times 4\times 5$). In the Sparsity
column, the numbers in parentheses indicate the number of non-zero variables
in each dimension. “Cor” is the correlation between the estimated linear
hyperplane and the true hyperplane. “Mis” is the average misclassification
rate. “TP” is the true positive rate, i.e., the proportion of non-zero
coefficients that are correctly estimated to be non-zero.“TN” is the true
negative rate, i.e., the proportion of zero coefficients that are correctly
estimated to be zero. The margins of error (2* standard deviations across 200
replicates) for each statistic are also listed following the $\pm$ symbol.

N | Sparsity | Methods | Cor | Mis | TP | TN | Prop. of Cor$>$0.5
---|---|---|---|---|---|---|---
40 | More (5 $\times$ 2 $\times$ 2) | M-SDWD | 0.829$\pm$0.028 | 0.159$\pm$0.030 | 0.616$\pm$0.077 | 0.768$\pm$0.079 | 0.31
|  | M-SDWD ($\lambda_{1}=0$) | 0.778$\pm$0.030 | 0.179$\pm$0.030 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.34
|  | M-DWD | 0.815$\pm$0.029 | 0.149$\pm$0.029 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.3
|  | Full SDWD | 0.721$\pm$0.038 | 0.159$\pm$0.033 | 0.277$\pm$0.052 | 0.955$\pm$0.019 | 0.22
| Less (5 $\times$ 4 $\times$ 5) | M-SDWD | 0.851$\pm$0.019 | 0.101$\pm$0.017 | 0.684$\pm$0.044 | 0.518$\pm$0.057 | 0.845
|  | M-SDWD ($\lambda_{1}=0$) | 0.853$\pm$0.018 | 0.102$\pm$0.017 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.865
|  | M-DWD | 0.877$\pm$0.017 | 0.096$\pm$0.016 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.84
|  | Full SDWD | 0.696$\pm$0.021 | 0.089$\pm$0.018 | 0.280$\pm$0.033 | 0.869$\pm$0.027 | 0.54
| No (15 $\times$ 4 $\times$ 5) | M-SDWD | 0.930$\pm$0.010 | 0.030$\pm$0.008 | 0.816$\pm$0.036 | - | 0.95
|  | M-SDWD ($\lambda_{1}=0$) | 0.935$\pm$0.011 | 0.031$\pm$0.009 | 1.000$\pm$0.000 | - | 0.96
|  | M-DWD | 0.947$\pm$0.010 | 0.027$\pm$0.009 | 1.000$\pm$0.000 | - | 0.94
|  | Full SDWD | 0.724$\pm$0.016 | 0.029$\pm$0.008 | 0.434$\pm$0.049 | - | 0.765
100 | More (5 $\times$ 2 $\times$ 2) | M-SDWD | 0.872$\pm$0.024 | 0.191$\pm$0.025 | 0.673$\pm$0.058 | 0.772$\pm$0.058 | 0.535
|  | M-SDWD ($\lambda_{1}=0$) | 0.850$\pm$0.023 | 0.186$\pm$0.024 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.515
|  | M-DWD | 0.853$\pm$0.026 | 0.190$\pm$0.025 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.525
|  | Full SDWD | 0.800$\pm$0.027 | 0.177$\pm$0.025 | 0.340$\pm$0.042 | 0.953$\pm$0.016 | 0.45
| Less (5 $\times$ 4 $\times$ 5) | M-SDWD | 0.929$\pm$0.011 | 0.090$\pm$0.014 | 0.749$\pm$0.033 | 0.536$\pm$0.051 | 0.945
|  | M-SDWD ($\lambda_{1}=0$) | 0.927$\pm$0.010 | 0.094$\pm$0.014 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.95
|  | M-DWD | 0.939$\pm$0.010 | 0.093$\pm$0.014 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.95
|  | Full SDWD | 0.768$\pm$0.018 | 0.106$\pm$0.016 | 0.336$\pm$0.032 | 0.838$\pm$0.028 | 0.83
| No (15 $\times$ 4 $\times$ 5) | M-SDWD | 0.962$\pm$0.009 | 0.031$\pm$0.010 | 0.879$\pm$0.030 | - | 0.99
|  | M-SDWD ($\lambda_{1}=0$) | 0.968$\pm$0.008 | 0.033$\pm$0.010 | 1.000$\pm$0.000 | - | 0.995
|  | M-DWD | 0.975$\pm$0.006 | 0.030$\pm$0.010 | 1.000$\pm$0.000 | - | 0.99
|  | Full SDWD | 0.816$\pm$0.013 | 0.036$\pm$0.010 | 0.437$\pm$0.037 | - | 0.925

###  5.3 Rank-R model simulation and results

In this simulation study, we aim to validate the higher rank multiway sparse
model with a scenario in which $R=2$. The procedure for data generation is
similar to that of the rank-1 model. Instead of generating vectors
$\mathbf{u}_{k}$ for each dimension, we generated three matrices
$\mathbf{U}_{k},k=1,2,3$ with dimensions $P_{k}\times 2$. The entries of
$\mathbf{U}_{k}$ are generated independently from a normal distribution
$N(0,1)$, then we compute
$\mu_{1}=\llbracket\mathbf{U}_{1},\dots,\mathbf{U}_{K}\rrbracket$. The $N_{1}$
samples corresponding to class $1$ were generated from a multivariate normal
distribution $N(\mu_{1},\mathbf{I}_{P_{1}P_{2}P_{3}\times P_{1}P_{2}P_{3}})$.
The $N_{0}$ samples in class -$1$ were generated independently from a $N(0,1)$
distribution. Different scenarios are considered in terms of dimensions,
sample size and sparsity levels. The results in Table 5 show that the rank-2
model has the best performance when the data were generated from a true rank-2
model under different scenarios. The performances of these methods are similar
for a higher rank (R=5) and low dimensional case (see Appendix E).

Table 5: Simulation results under the high dimensional scenario ($30\times
15\times 15$) when the true model is rank-2. In the Sparsity column, the
numbers in parentheses indicate the number of non-zero variables in each
dimension. “Cor” is the correlation between the estimated linear hyperplane
and the true hyperplane. “Mis” is the average misclassification rate. “TP” is
the true positive rate, i.e., the proportion of non-zero coefficients that are
correctly estimated to be non-zero.“TN” is the true negative rate, i.e., the
proportion of zero coefficients that are correctly estimated to be zero. The
margins of error (2* standard deviations across 200 replicates) for each
statistic are also listed following the $\pm$ symbol.

N | Sparsity | Methods | Cor | Mis | TP | TN
---|---|---|---|---|---|---
40 | More ($5\times 5\times 5$) | M-SDWD (R=2) | 0.833$\pm$0.020 | 0.027$\pm$0.011 | 0.714$\pm$0.038 | 0.761$\pm$0.043
|  | M-SDWD ($\lambda_{1}=0$, R=2) | 0.747$\pm$0.025 | 0.040$\pm$0.013 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.779$\pm$0.025 | 0.035$\pm$0.014 | 0.678$\pm$0.039 | 0.716$\pm$0.050
|  | M-SDWD ($\lambda_{1}=0$, R=1) | 0.745$\pm$0.026 | 0.043$\pm$0.015 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.736$\pm$0.033 | 0.055$\pm$0.019 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.591$\pm$0.024 | 0.071$\pm$0.015 | 0.150$\pm$0.013 | 0.984$\pm$0.007
| Less ($10\times 15\times 15$) | M-SDWD (R=2) | 0.984$\pm$0.003 | 0.000$\pm$0.000 | 0.942$\pm$0.016 | 0.402$\pm$0.060
|  | M-SDWD ($\lambda_{1}=0$, R=2) | 0.990$\pm$0.002 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.813$\pm$0.011 | 0.000$\pm$0.000 | 0.840$\pm$0.029 | 0.488$\pm$0.062
|  | M-SDWD ($\lambda_{1}=0$, R=1) | 0.821$\pm$0.010 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.786$\pm$0.016 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.750$\pm$0.009 | 0.000$\pm$0.000 | 0.239$\pm$0.018 | 0.937$\pm$0.012
| No ($30\times 15\times 15$) | M-SDWD (R=2) | 0.996$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-SDWD ($\lambda_{1}=0$, R=2) | 0.996$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-SDWD (R=1) | 0.790$\pm$0.012 | 0.000$\pm$0.000 | 0.836$\pm$0.033 | -
|  | M-SDWD ($\lambda_{1}=0$, R=1) | 0.799$\pm$0.010 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.763$\pm$0.016 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.784$\pm$0.006 | 0.000$\pm$0.000 | 0.399$\pm$0.036 | -
100 | More ($5\times 5\times 5$) | M-SDWD (R=2) | 0.925$\pm$0.010 | 0.012$\pm$0.004 | 0.841$\pm$0.026 | 0.754$\pm$0.041
|  | M-SDWD ($\lambda_{1}=0$, R=2) | 0.889$\pm$0.011 | 0.016$\pm$0.005 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.853$\pm$0.012 | 0.015$\pm$0.004 | 0.738$\pm$0.031 | 0.818$\pm$0.039
|  | M-SDWD ($\lambda_{1}=0$, R=1) | 0.829$\pm$0.014 | 0.019$\pm$0.005 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.825$\pm$0.020 | 0.027$\pm$0.010 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.797$\pm$0.014 | 0.026$\pm$0.007 | 0.236$\pm$0.016 | 0.992$\pm$0.002
| Less ($10\times 15\times 15$) | M-SDWD (R=2) | 0.992$\pm$0.001 | 0.000$\pm$0.000 | 0.958$\pm$0.012 | 0.456$\pm$0.061
|  | M-SDWD ($\lambda_{1}=0$, R=2) | 0.996$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.807$\pm$0.010 | 0.000$\pm$0.000 | 0.871$\pm$0.024 | 0.526$\pm$0.062
|  | M-SDWD ($\lambda_{1}=0$, R=1) | 0.807$\pm$0.010 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.797$\pm$0.012 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.856$\pm$0.005 | 0.000$\pm$0.000 | 0.346$\pm$0.019 | 0.935$\pm$0.012
| No ($30\times 15\times 15$) | M-SDWD (R=2) | 0.998$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-SDWD ($\lambda_{1}=0$, R=2) | 0.998$\pm$0.000 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-SDWD (R=1) | 0.786$\pm$0.009 | 0.000$\pm$0.000 | 0.908$\pm$0.023 | -
|  | M-SDWD ($\lambda_{1}=0$, R=1) | 0.790$\pm$0.009 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.758$\pm$0.014 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.871$\pm$0.004 | 0.000$\pm$0.000 | 0.466$\pm$0.029 | -

##  6 Applications

###  6.1 MRS data

We apply the proposed method to our motivating application with MRS data to
illustrate its utility. Friedreich’s ataxia (FRDA) is an early-onset
neurodegenerative disease caused by abnormalities in the frataxin gene
(Pandolfo, 2008) resulting in knockdown of frataxin protein. The motivation of
this study is to assess brain metabolic changes in a transgenic mouse model in
which frataxin knockdown can be turned on by giving the antibiotic
doxycycline. Wild type (WT) mice and transgenic (TG) mice were randomly
assigned to two different treatment groups: doxycycline treated group
($N_{1}=11$) and controls ($N_{0}=10$). Thus, the study included 4 groups of
mice: WT treated with dox (WT-dox), WT without dox (WT-nodox), TG treated with
dox (TG-dox), TG without dox (TG-nodox). Treatment was applied during weeks
1-12, then stopped during weeks 12-24 (recovery). Frataxin knockdown during
weeks 1-12 was expected only in the TG-dox group, but not in WT animals nor in
animals not receiving dox. The concentrations of 13 metabolites were measured
in three different regions (cerebellum, cortex, and cervical spine). These
animals were scanned at three time points: 0 weeks, 12 weeks and 24 weeks. The
data have a multiway structure with four dimensions: mice $\times$ metabolites
$\times$ regions $\times$ time, $\mathbb{X}:21\times 13\times 3\times 3$. Our
goal is to summarize the signal distinguishing the two treatment groups (dox
vs. no dox) across the different dimensions, and assess how the effect of
treatment depends on WT/TG status.

The penalty parameters for multiway sparse DWD were selected through 10-fold
cross validation (see Appendix G). Using the selected tuning parameters,
Figure 2 shows the DWD scores computed under leave-one-out cross-validation
based on the rank-1 multiway sparse DWD model. These scores show robust
separation of the no-dox and dox groups (which correspond to the two classes
$y_{i}=-1$ and $y_{i}=1$) but only for the transgenic (TG) mice. The three
scores for the TG-nodox, WT-nodox, and WT-dox groups are very similar. Because
we train the classifier on both genetic cohorts but only observe a difference
in the transgenic mice, this confirms the hypothesis that wild type mice are
not affected by treatment; further, the results illustrate the robustness of
the approach to within-class heterogeneity.

In order to capture the uncertainty of the estimated weights and construct
$95\%$ confidence intervals, 500 bootstrap samples were generated. For each
bootstrap sample, 21 mice were resampled with replacement, and then the model
was fit to the bootstrap sample to estimate weights for each dimension. The
$2.5\%$ and $97.5\%$ quantiles for each estimated weight across the 500
bootstrap samples were computed to construct the $95\%$ confidence interval.
Figure 3 shows the estimated weights and their $95\%$ bootstrap confidence
intervals. The metabolites with large absolute weights, such as PCho+GPC,
Cr+PCr, and Ins, are considered important for Friedreich’s ataxia research
(França et al., 2009; Iltis et al., 2010; Gramegna et al., 2017), so were
expected to have high weights for distinguishing the two groups. Two
metabolites (Lac and Glu) and one time point (0 weeks) did not inform the
classification as their estimated weights are exactly 0 in the full data fit.
Having no effect at baseline (0 weeks) makes sense as this is prior to the dox
treatment. The observed changes at the later time points are similar to those
observed in the R6/2 mouse, a severe mouse model of Huntington’s disease
(Zacharoff et al., 2012). This is consistent with the fact that both
Friedreich’s ataxia and Huntington’s disease are characterized by impairment
in energy metabolism.

Figure 2: FRDA mouse study: Rank-1 multiway sparse DWD scores under leave-one-
out cross-validation to classify mice that did or did not receive dox
treatment. The transgenic (TG) mice that received dox are clearly
distinguished from the TG mice that did not and from the wild-type (WT) mice.
A kernel density estimate is shown for the scores in each subgroup. Figure 3:
FRDA mouse study: Rank-1 multiway sparse DWD weights for metabolites, regions
and time with $95\%$ bootstrap confidence intervals.

We compared misclassification rates for the M-SDWD, M-SDWD with
$\lambda_{1}=0$, M-DWD, CATCH, Full SDWD and Random Forests by 10-fold cross
validation. The 21 mice were randomly partitioned into 10 test subgroups of
approximately equal size. This procedure was repeated $100$ times with
different random partitions and the average test misclassification rates for
dox vs. no dox were as follows: $22.1\%$ for M-SDWD, $20.9\%$ for M-SDWD with
$\lambda_{1}=0$, $23.6\%$ for M-DWD, $22.3\%$ for CATCH, $42.9\%$ for Full
SDWD, and $22.3\%$ for Random Forest. To better reflect the reality of the
treatment effect, we further considered the misclassification rate for TG-dox
vs. the other three groups under this approach: $3.1\%$ for M-SDWD, $1.9\%$
for M-SDWD with $\lambda_{1}=0$, $4.6\%$ for M-DWD, $4.0\%$ for CATCH,
$24.6\%$ for Full SDWD, and $10.5\%$ for Random Forest. Overall, the multiway
methods perform better than the non-multiway methods (Full SDWD and Random
Forest), illustrating the advantages of accounting for multiway structure.
Moreover, the factorization of weights across dimensions in Figure 3 is useful
for interpretation, and is not provided by the CATCH model. We conclude that
the performance of M-SDWD model is competitive with alternatives. The rank-2
multiway sparse model was also applied, resulting in a higher
misclassification rate ($29.9\%$ for dox vs. no dox and $10.8\%$ for TG-dox
vs. the other three groups), indicating that the rank-1 model is sufficient.

###  6.2 Gene course time data

We further applied multiway sparse DWD to the gene expression time course data
described in Baranzini et al. (2004). The purpose of this study is to classify
clinical response to treatment for multiple sclerosis (MS) patients.
Recombinant human interferon beta (rIFN$\beta$) was given to 53 patients for
controlling the symptoms of MS. For each patient, gene expression was measured
for 76 genes at 7 time points: baseline (i.e. before treatment) and 6 follow-
up time points (3 months, 6 months, 9 months, 12 months, 18 months, and 24
months). The data are a multi-way array with 3 dimensions: patients $\times$
genes $\times$ time points, $\mathbb{X}:53\times 76\times 7$. Based on
clinical characteristics, each patient was designated as a good responder or
poor responder to rIFN$\beta$. The proposed rank-1 multiway sparse DWD model
was used to differentiate good responders from poor responders. Figure 4 shows
the DWD scores under leave-one-out cross validation for the rank-1 multiway
model. The groups of these good and poor responders are nearly perfectly
separated. The coefficient estimates and 95$\%$ bootstrap confidence intervals
for each gene and each time point are shown in Figure 5. The coefficients
across time had little variability, which suggests that the distinction
between good and poor responders is not driven by effects that vary over the
time course. This agrees with the results reported in Baranzini et al. (2004)
where they found there was no group*time interaction effects. Note that the
data are not very sparse as almost all estimated coefficients for gene and
time are not zero. We also compared the performance of the rank-1 model with
other competitors. The mean misclassification rates were $13.0\%$ for M-SDWD,
$12.6\%$ for M-SDWD with $\lambda_{1}=0$, $17.4\%$ for M-DWD, $20.8\%$ for
CATCH, $20.5\%$ for Full SDWD and $31.5\%$ for Random Forest. The rank-1
multiway sparse models have the best classification over other methods. It is
notable that the multiway sparse DWD model outperforms multiway DWD even when
the degree of sparsity in the data is low. We also considered the rank-2
multiway sparse model to classify good and poor responders. A rank-2 model
performs worse than the rank-1 model with a higher misclassification rate
$(18.7\%)$, which agrees with the results in Lyu et al. (2017).

Figure 4: MS treatment study: Rank-1 multiway sparse DWD scores under leave-
one-out cross-validation for good and poor treatment responders, with a kernel
density estimate for each group.  Figure 5: MS treatment study: Rank-1
multiway sparse DWD weights for genes and time with $95\%$ bootstrap
confidence intervals.

##  7 Discussion

We have proposed a general framework for high dimensional classification on
multiway data with any number of dimensions, which can account for sparsity in
the model. Both the simulation and data analysis results have shown that the
proposed multiway sparse DWD model can improve classification accuracy when
the underlying signal is sparse or has multiway structure. The proposed method
is robust to any degree of sparsity in the model, which has been demonstrated
in both simulation and applications. Moreover, the method is robust to the
complexity of the shared signal across the different dimensions, by allowing
for higher rank models. Lastly, the use of multiway structure and sparsity can
facilitate and simplify interpretation. For our motivating application in
Section 6.1, a single application of multiway sparse DWD provided the
following insights: (1) TG mice show a distinct brain metabolomic profile
after receiving doxycycline but WT mice do not, (2) these changes to the
metabolite profile are similar across three neurological regions but most
pronounced in the cerebellum, and (3) these changes are most prominent at the
end of 12 weeks of treatment, but subside post-treatment.

Despite the flexibility of multiway sparse DWD, the results are sensitive to
the choice of rank and tuning parameters. The simulation in Section 5.1 shows
that when the model is not sparse, the rank-1 model with $\lambda_{1}=0$
performs better than the rank-1 model with $\lambda_{1}$ selected by cross
validation; thus, the sparsity penalty only improves results when the signal
is truly sparse. Also, the rank-1 model performs poorly when the true model
has higher rank, as expected. Thus, in practice we suggest applying different
versions of the method and comparing their performance for any given data
application. Also, as discussed in Sections 4.3, 5.2, and Appendix D,
algorithms using cyclic coordinate descent, including ours, may not converge
to a global optimum. In practice, we suggest running the algorithm from
multiple starting values to assess the stability of the results and (if
results are sensitive to starting values) select the solution that yields the
smallest value for the objective function. In this article we only consider
binary classification, and extensions to multi-category classification (e.g.,
as in Huang et al. (2013)) is a worthwhile future direction.

## Acknowledgements

This work was supported in part by the National Institutes of Health (NIH)
grant R01GM130622. The study on Friedreich’s Ataxia mice was supported by the
Friedreich’s Ataxia Research Alliance (FARA).

## Appendices:

##  Appendix A Evaluation of cross-validation method for selecting tuning
parameters

We conducted a simulation study to assess how the choices of tuning parameters
$\lambda_{1}$ and $\lambda_{2}$ affect the performance of the proposed method.
A dataset with sample size $N=100$ was generated, with two classes of equal
size $(N_{0}=N_{1}=50)$. The predictors have the form of a three-way array of
dimensions $\mathbb{X}:P_{1}\times P_{2}\times P_{3}$ where $P_{1}=4,P_{2}=5,$
and $P_{3}=15$. The $N_{0}$ samples corresponding to class -1 were generated
from a multivariate normal distribution
N($\mathbf{0}$,$\mathbf{I}_{P_{1}P_{2}P_{3}\times P_{1}P_{2}P_{3}}$). The
other $N_{1}$ samples corresponding to class 1 were generated from a
multivariate normal distribution
N($\bm{\mu}_{1}$,$\mathbf{I}_{P_{1}P_{2}P_{3}\times P_{1}P_{2}P_{3}}$) where
$\bm{\mu}_{1}=\mathbf{u}_{1}\circ\mathbf{u}_{2}\circ\mathbf{u}_{3}$. Here
$\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ were generated from
$N(\mathbf{0},\mathbf{I}_{P_{k}\times P_{k}})$, $k=1,2$. For $\mathbf{u}_{3}$
5 values are set to zero and 10 nonzero values are generated from
$N(\mathbf{0},\mathbf{I}_{10\times 10})$. We consider 4 candidates for
$\lambda_{1}$ $(10^{-4},0.001,0.005,0.01)$ and 3 candidates for $\lambda_{2}$
$(1,0.5,0.1)$. We applied the multiway sparse DWD method to the simulated data
with fixed parameters and selected parameters by cross-validation. The
simulation is repeated 100 times. The results are shown in Table 6 and Table 7
with average correlations between true values and estimates for
$\mathbf{u}_{k}$ and average percentages of zero or non-zero coefficients in
$\mathbf{u}_{3}$ that are correctly estimated. The accuracy of classification
decreases as $\lambda_{2}$ decreases. For fixed $\lambda_{2}$, as
$\lambda_{1}$ becomes larger, more zero coefficients for $\mathbf{u}_{3}$ are
correctly shrunk to 0. Table 7 shows simulations where the correlations
between true values and estimates are very high, and the classification is
accurate using the cross-validation method to select parameters, although the
selected parameters might be slightly different for each replicate.

Table 6: Simulation results based on multiway sparse DWD with prespecified
$\lambda_{1}$ and $\lambda_{2}$ candidates. “Cor($\mathbf{u}_{k}$)” is the
correlation between the estimated linear hyperplane and the true hyperplane
for $k^{th}$ dimension. “TP($\mathbf{u_{3}}$)” is the true positive rate that
is the proportion of non-zero coefficients in $\mathbf{u}_{3}$ are correctly
estimated to be non-zero.“TN($\mathbf{u}_{3}$)” is the true negative rate that
is the proportion of zero coefficients in $\mathbf{u}_{3}$ are correctly
estimated to be zero. The results are the mean over 100 replicates.

$\lambda_{2}$ | 1 | 0.5 | 0.1
---|---|---|---
$\lambda_{1}$ | 1e-4 | 0.001 | 0.005 | 0.01 | 1e-4 | 0.001 | 0.005 | 0.01 | 1e-4 | 0.001 | 0.005 | 0.01
Cor($\mathbf{u}_{1}$) | 0.96 | 0.96 | 0.97 | 0.97 | 0.93 | 0.93 | 0.93 | 0.93 | 0.89 | 0.91 | 0.93 | 0.92
Cor($\mathbf{u}_{2}$) | 1.00 | 1.00 | 1.00 | 0.99 | 1.00 | 0.99 | 0.99 | 0.98 | 0.95 | 0.95 | 0.95 | 0.93
Cor($\mathbf{u}_{3}$) | 0.98 | 0.99 | 0.98 | 0.98 | 0.97 | 0.97 | 0.97 | 0.96 | 0.94 | 0.93 | 0.93 | 0.92
TP ($\mathbf{u}_{3}$) | 1.00 | 0.98 | 0.81 | 0.85 | 1.00 | 0.98 | 0.90 | 0.83 | 0.99 | 0.94 | 0.80 | 0.72
TN ($\mathbf{u}_{3}$) | 0.03 | 0.37 | 0.81 | 0.91 | 0.03 | 0.41 | 0.81 | 0.92 | 0.09 | 0.59 | 0.90 | 0.95

Table 7: Simulation results based on multiway sparse DWD with selected penalty
parameters by cross-validation. $\lambda_{1}^{G}$ and $\lambda_{2}^{G}$ denote
the geometric means of selected parameters over 100 replicates.

| $\lambda_{1}^{G}=0.0008,\lambda_{2}^{G}=1.56$
---|---
Cor($\mathbf{u}_{1}$) | 0.98
Cor($\mathbf{u}_{2}$) | 0.99
Cor($\mathbf{u}_{3}$) | 0.99
TP($\mathbf{u}_{3}$) | 0.39
TN($\mathbf{u}_{3}$) | 0.93

##  Appendix B Simulation results based on the objective function with
separable $L_{2}$ penalty

In this section, we consider a simulation study to evaluate the performance of
the higher rank model based on the objective function with separable $L_{2}$
penalty (i.e. $P^{U}_{\lambda_{1},\lambda_{2}}(\mathbb{B})$ defined in the
main article). We observed that the separable $L_{2}$ penalty tends to reduce
the rank of the estimated coefficient array, motivating us to use the non-
separable $L_{2}$ penalty (i.e. $P_{\lambda_{1},\lambda_{2}}(\mathbb{B})$) for
the higher rank model. We simulated 200 datasets with high dimensions
($30\times 15\times 15$) under the case of no sparsity. The process of
generating multiway data is identical to that for the higher rank model
($R=2$) described in Section 5.3 of the main article. Table 8 shows average
correlations, misclassification rates and true positive rates over 200
simulation replicates. Multiway sparse DWD (M-SDWD) with $R=2$ has similar
performance to M-SDWD with R=1, which implies M-SDWD (R=2) may not recover the
true rank specified in the model. To measure the rank of the estimated model,
we conduct SVD to obtain singular values for the estimated coefficient matrix.
In most simulations, the second singular values of the estimated coefficient
matrix are zero or nearly zero, as shown in Figure 6, which indicates the
solutions of higher rank (R=2) model are often shrunk to a lower rank. The
proportion of simulations that give estimates with true rank $R=2$ is only
$14.5\%$. Table 9 shows the results of simulations that can truly detect rank
2 components, and we can see the model with R=2 performs much better than
rank-1 multiway models and full SDWD.

Figure 6: Histograms of the two singular values computed based on the estimated coefficient matrix. Table 8: Simulation results based on all simulations under the high dimensional scenario: “Cor” is the correlation between the estimated linear hyperplane and the true hyperplane. “Mis” is the average misclassification rate. “TP” is the true positive rate, i.e., the proportion of non-zero coefficients that are correctly estimated to be non-zero.The margins of error (2* standard deviations across 200 replicates) for each statistic are also listed following the $\pm$. Methods | Cor | Mis | TP
---|---|---|---
M-SDWD (R=2) | 0.833$\pm$0.012 | 0.000$\pm$0.000 | 0.883$\pm$0.029
M-SDWD (R=1) | 0.792$\pm$0.011 | 0.000$\pm$0.000 | 0.872$\pm$0.029
M-SDWD ($\lambda_{1}=0$, R=1) | 0.800$\pm$0.010 | 0.000$\pm$0.000 | 1.000$\pm$0.000
M-DWD | 0.771$\pm$0.014 | 0.000$\pm$0.000 | 1.000$\pm$0.000
Full SDWD | 0.781$\pm$0.006 | 0.000$\pm$0.000 | 0.419$\pm$0.036
Table 9: Simulation results among simulations that give estimates with true rank (R=2) under the high-dimensional scenario: “Cor” is the correlation between the estimated linear hyperplane and the true hyperplane. “Mis” is the average misclassification rate. “TP” is the true positive rate, i.e., the proportion of non-zero coefficients that are correctly estimated to be non-zero.“TN” is the true negative rate, i.e., the proportion of zero variables that are correctly estimated to be zero. The margins of error (2* standard deviations across 200 replicates) for each statistic are also listed following the $\pm$ symbol. Methods | Cor | Mis | TP |  $\%$ rank 2
---|---|---|---|---
M-SDWD (R=2) | 0.984$\pm$0.008 | 0.000$\pm$0.000 | 0.955$\pm$0.040 | 0.145
M-SDWD (R=1) | 0.709$\pm$0.033 | 0.000$\pm$0.000 | 0.844$\pm$0.088 | -
M-SDWD ($\lambda_{1}=0$, R=1) | 0.726$\pm$0.023 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
M-DWD | 0.676$\pm$0.024 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
Full SDWD | 0.782$\pm$0.016 | 0.000$\pm$0.000 | 0.492$\pm$0.088 | -

##  Appendix C Simulation results under different signal to noise ratios

We conducted more simulation studies to compare the proposed method with other
methods under different signal to noise ratios (SNR= (0.1,0.2,0.3,0.5)). We
considered the high-dimensional, more sparsity case of Section 5.1 in the main
article for this simulation. Tables 10 and 11 show the results based on all
simulations and those simulations with correlations larger than 0.5,
respectively. The multiway sparse DWD model performs better than other methods
in terms of its correlation with the true hyperplane, and has competive
misclassification rates across the different SNR levels. The last column in
Table 11 show the proportions of simulations with correlations larger than 0.5
for different methods under different SNRs. As the SNR increases, the
proportion increases as well, as the algorithm tends to converge to the true
solution with more signal.

Table 10: Simulation results for different signal to noise ratios: “Cor” is the correlation between the estimated linear hyperplane and the true hyperplane. “Mis” is the average misclassification rate. “TP” is the true positive rate, i.e., the proportion of non-zero coefficients that are correctly estimated to be non-zero.“TN” is the true negative rate, i.e., the proportion of zero variables that are correctly estimated to be zero. The margins of error (2* standard deviations across 200 replicates) for each statistic are also listed following the $\pm$ symbol. SNR | Methods | Cor | Mis | TP | TN
---|---|---|---|---|---
0.1 | M-SDWD | 0.717$\pm$0.053 | 0.225$\pm$0.029 | 0.484$\pm$0.050 | 0.837$\pm$0.038
| M-SDWD ($\lambda_{1}=0$) | 0.594$\pm$0.054 | 0.220$\pm$0.028 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| M-DWD | 0.560$\pm$0.058 | 0.235$\pm$0.030 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| Full SDWD | 0.469$\pm$0.040 | 0.251$\pm$0.025 | 0.168$\pm$0.028 | 0.917$\pm$0.029
0.2 | M-SDWD | 0.849$\pm$0.038 | 0.089$\pm$0.020 | 0.668$\pm$0.040 | 0.806$\pm$0.041
| M-SDWD ($\lambda_{1}=0$) | 0.796$\pm$0.040 | 0.101$\pm$0.021 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| M-DWD | 0.766$\pm$0.049 | 0.121$\pm$0.025 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| Full SDWD | 0.636$\pm$0.035 | 0.136$\pm$0.022 | 0.172$\pm$0.023 | 0.957$\pm$0.022
0.3 | M-SDWD | 0.879$\pm$0.035 | 0.066$\pm$0.018 | 0.720$\pm$0.036 | 0.776$\pm$0.045
| M-SDWD ($\lambda_{1}=0$) | 0.838$\pm$0.037 | 0.077$\pm$0.020 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| M-DWD | 0.831$\pm$0.042 | 0.084$\pm$0.021 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| Full SDWD | 0.690$\pm$0.031 | 0.104$\pm$0.020 | 0.165$\pm$0.017 | 0.982$\pm$0.011
0.5 | M-SDWD | 0.927$\pm$0.026 | 0.039$\pm$0.014 | 0.777$\pm$0.033 | 0.762$\pm$0.048
| M-SDWD ($\lambda_{1}=0$) | 0.890$\pm$0.030 | 0.048$\pm$0.016 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| M-DWD | 0.884$\pm$0.036 | 0.053$\pm$0.018 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| Full SDWD | 0.760$\pm$0.027 | 0.066$\pm$0.017 | 0.189$\pm$0.017 | 0.993$\pm$0.002
Table 11: Simulation results among simulations with correlation greater than
0.5 under different signal to noise ratios: “Cor” is the correlation between
the estimated linear hyperplane and the true hyperplane. “Mis” is the average
misclassification rate. “TP” is the true positive rate, i.e., the proportion
of non-zero coefficients that are correctly estimated to be non-zero.“TN” is
the true negative rate, i.e., the proportion of zero variables that are
correctly estimated to be zero. The margins of error (2* standard deviations
across 200 replicates) for each statistic are also listed following the $\pm$
symbol.

SNR | Methods | Cor | Mis | TP | TN | % Cor$>$0.5
---|---|---|---|---|---|---
0.1 | M-SDWD | 0.916$\pm$0.014 | 0.077$\pm$0.015 | 0.628$\pm$0.046 | 0.869$\pm$0.040 | 0.645
| M-SDWD ($\lambda_{1}=0$) | 0.858$\pm$0.017 | 0.088$\pm$0.016 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.66
| M-DWD | 0.878$\pm$0.018 | 0.080$\pm$0.016 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.62
| Full SDWD | 0.718$\pm$0.023 | 0.096$\pm$0.017 | 0.124$\pm$0.014 | 0.993$\pm$0.002 | 0.52
0.2 | M-SDWD | 0.935$\pm$0.010 | 0.048$\pm$0.010 | 0.697$\pm$0.038 | 0.831$\pm$0.042 | 0.905
| M-SDWD ($\lambda_{1}=0$) | 0.899$\pm$0.012 | 0.051$\pm$0.010 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.88
| M-DWD | 0.926$\pm$0.010 | 0.042$\pm$0.009 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.825
| Full SDWD | 0.760$\pm$0.017 | 0.060$\pm$0.011 | 0.156$\pm$0.016 | 0.994$\pm$0.001 | 0.765
0.3 | M-SDWD | 0.948$\pm$0.008 | 0.033$\pm$0.008 | 0.745$\pm$0.034 | 0.795$\pm$0.046 | 0.925
| M-SDWD ($\lambda_{1}=0$) | 0.921$\pm$0.010 | 0.035$\pm$0.009 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.905
| M-DWD | 0.941$\pm$0.008 | 0.030$\pm$0.008 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.88
| Full SDWD | 0.781$\pm$0.016 | 0.048$\pm$0.010 | 0.170$\pm$0.016 | 0.994$\pm$0.002 | 0.825
0.5 | M-SDWD | 0.963$\pm$0.006 | 0.021$\pm$0.007 | 0.792$\pm$0.030 | 0.767$\pm$0.049 | 0.96
| M-SDWD ($\lambda_{1}=0$) | 0.943$\pm$0.009 | 0.022$\pm$0.007 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.94
| M-DWD | 0.958$\pm$0.007 | 0.018$\pm$0.006 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.92
| Full SDWD | 0.816$\pm$0.015 | 0.030$\pm$0.009 | 0.202$\pm$0.018 | 0.994$\pm$0.001 | 0.89

##  Appendix D Assessment of convergence

We conducted more simulations to explore how the signal to noise ratio (SNR)
affects the convergence of the proposed methods and other alternatives. Figure
7 and 8 show the distributions of correlations between the true hyperplane and
the estimated hyperplane for SNR = 0.1 and SNR=0.3, under the more sparsity
case based on high dimensional data ($30\times 15\times 15$) with sample size
$N=100$. Combining with results for SNR=0.2 shown in the main article we can
see as the SNR increases, the number of correlations with small values
decreases, and the convergence issue becomes less severe.

Figure 7: Histogram of correlations between true hyperplane and estimates with
four classification methods under high-dimensional scenario with SNR=0.1 and
sample size N= 100. Figure 8: Histogram of correlations between true
hyperplane and estimates with four classification methods under high-
dimensional scenario with SNR=0.3 and sample size N= 100.

##  Appendix E Additional simulations with Rank $>$ 1

We conducted more simulations to evaluate the performances of the rank-$R$
multiway sparse model. These simulations were analogous to those in Section
5.3 of the main article, with different manipulated conditions. Table 12 gives
the results under a lower dimensional setting ($\mathbb{X}:15\times 4\times
5$). Table 13 gives results for both the lower-dimensional and higher-
dimensional settings for a model with rank $R=5$. From these results we can
conclude that the rank-$R$ model can generally perform well when the data were
generated under a true rank-$R$ model, although lower-rank approximations
perform comparably well when the dimension is small and the signal is sparse.

Table 14 shows the mean correlations between the estimated linear hyperplane
and the true hyperplane using the proposed method (M-SDWD) by assumed rank
($\hat{R}=1,2,3,4,5$) for different true ranks under the high dimensional
scenario with more sparsity when N = 40. In general, the correlation is
maximized when the assumed rank is equal to the true rank; however,
performance tends to be robust to deviations from the precise value of the
true rank, especially when the true rank is higher.

Table 12: Simulation results under the low dimensional scenario ($15\times
4\times 5$) when the true model is rank-2. In the Sparsity column, the numbers
in parentheses indicate the number of non-zero variables in each dimension.
“Cor” is the correlation between the estimated linear hyperplane and the true
hyperplane. “Mis” is the average misclassification rate. “TP” is the true
positive rate, i.e., the proportion of non-zero coefficients that are
correctly estimated to be non-zero.“TN” is the true negative rate, i.e., the
proportion of zero coefficients that are correctly estimated to be zero. The
margins of error (2* standard deviations across 200 replicates) for each
statistic are also listed following the $\pm$ symbol.

N | Sparsity | Methods | Cor | Mis | TP | TN
---|---|---|---|---|---|---
40 | More ($5\times 2\times 2$) | M-SDWD (R=2) | 0.716$\pm$0.034 | 0.191$\pm$0.027 | 0.659$\pm$0.052 | 0.678$\pm$0.048
|  | M-SDWD ($\lambda_{1}=0$, R=2 ) | 0.643$\pm$0.032 | 0.183$\pm$0.023 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.746$\pm$0.038 | 0.157$\pm$0.024 | 0.676$\pm$0.052 | 0.711$\pm$0.056
|  | M-SDWD ($\lambda_{1}=0$, R=1 ) | 0.713$\pm$0.036 | 0.167$\pm$0.024 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.724$\pm$0.041 | 0.168$\pm$0.025 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.618$\pm$0.033 | 0.188$\pm$0.023 | 0.410$\pm$0.034 | 0.861$\pm$0.033
| Less ($5\times 4\times 5$) | M-SDWD (R=2) | 0.875$\pm$0.014 | 0.038$\pm$0.010 | 0.786$\pm$0.036 | 0.440$\pm$0.052
|  | M-SDWD ($\lambda_{1}=0$, R=2 ) | 0.870$\pm$0.015 | 0.036$\pm$0.009 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.810$\pm$0.017 | 0.043$\pm$0.010 | 0.698$\pm$0.042 | 0.530$\pm$0.054
|  | M-SDWD ($\lambda_{1}=0$, R=1 ) | 0.818$\pm$0.015 | 0.040$\pm$0.009 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.824$\pm$0.017 | 0.039$\pm$0.010 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.703$\pm$0.019 | 0.064$\pm$0.014 | 0.334$\pm$0.032 | 0.847$\pm$0.033
| No ($15\times 4\times 5$) | M-SDWD (R=2) | 0.958$\pm$0.007 | 0.003$\pm$0.003 | 1.000$\pm$0.000 | -
|  | M-SDWD ($\lambda_{1}=0$, R=2 ) | 0.959$\pm$0.006 | 0.003$\pm$0.002 | 1.000$\pm$0.000 | -
|  | M-SDWD (R=1) | 0.840$\pm$0.014 | 0.006$\pm$0.003 | 0.839$\pm$0.034 | -
|  | M-SDWD ($\lambda_{1}=0$, R=1 ) | 0.849$\pm$0.012 | 0.005$\pm$0.003 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.843$\pm$0.014 | 0.005$\pm$0.003 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.773$\pm$0.012 | 0.011$\pm$0.005 | 0.469$\pm$0.038 | -
100 | More ($5\times 2\times 2$) | M-SDWD (R=2) | 0.853$\pm$0.021 | 0.144$\pm$0.023 | 0.764$\pm$0.046 | 0.702$\pm$0.051
|  | M-SDWD ($\lambda_{1}=0$, R=2 ) | 0.791$\pm$0.026 | 0.156$\pm$0.023 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.862$\pm$0.022 | 0.141$\pm$0.020 | 0.735$\pm$0.052 | 0.763$\pm$0.059
|  | M-SDWD ($\lambda_{1}=0$, R=1 ) | 0.840$\pm$0.022 | 0.148$\pm$0.021 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.834$\pm$0.031 | 0.148$\pm$0.022 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.798$\pm$0.027 | 0.160$\pm$0.025 | 0.395$\pm$0.037 | 0.934$\pm$0.025
| Less ($5\times 4\times 5$) | M-SDWD (R=2) | 0.941$\pm$0.009 | 0.027$\pm$0.007 | 0.858$\pm$0.029 | 0.410$\pm$0.050
|  | M-SDWD ($\lambda_{1}=0$, R=2 ) | 0.938$\pm$0.009 | 0.027$\pm$0.007 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.860$\pm$0.012 | 0.030$\pm$0.007 | 0.772$\pm$0.033 | 0.518$\pm$0.053
|  | M-SDWD ($\lambda_{1}=0$, R=1 ) | 0.862$\pm$0.012 | 0.031$\pm$0.007 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.868$\pm$0.012 | 0.031$\pm$0.007 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.838$\pm$0.011 | 0.037$\pm$0.009 | 0.371$\pm$0.025 | 0.893$\pm$0.019
| No ($15\times 4\times 5$) | M-SDWD (R=2) | 0.976$\pm$0.004 | 0.004$\pm$0.003 | 0.964$\pm$0.012 | -
|  | M-SDWD ($\lambda_{1}=0$, R=2 ) | 0.979$\pm$0.004 | 0.004$\pm$0.003 | 1.000$\pm$0.000 | -
|  | M-SDWD (R=1) | 0.849$\pm$0.013 | 0.007$\pm$0.003 | 0.880$\pm$0.026 | -
|  | M-SDWD ($\lambda_{1}=0$, R=1 ) | 0.854$\pm$0.012 | 0.007$\pm$0.003 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.846$\pm$0.015 | 0.007$\pm$0.003 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.860$\pm$0.009 | 0.009$\pm$0.005 | 0.543$\pm$0.037 | -

Table 13: Simulation results when the true model is rank-5 (R=5) under high
and low dimensional scenarios and sample size $N=40$. In the Sparsity column,
the numbers in parentheses indicate the number of non-zero variables in each
dimension. “Cor” is the correlation between the estimated linear hyperplane
and the true hyperplane. “Mis” is the average misclassification rate. “TP” is
the true positive rate, i.e., the proportion of non-zero coefficients that are
correctly estimated to be non-zero.“TN” is the true negative rate, i.e., the
proportion of zero coefficients that are correctly estimated to be zero. The
margins of error (2* standard deviations across 200 replicates) for each
statistic are also listed following the $\pm$ symbol.

Dimensions | Sparsity | Methods | Cor | Mis | TP | TN
---|---|---|---|---|---|---
High ($30\times 15\times 15$) | More ($5\times 5\times 5$) | M-SDWD (R=5) | 0.845$\pm$0.013 | 0.002$\pm$0.002 | 0.927$\pm$0.016 | 0.623$\pm$0.048
|  | M-SDWD ($\lambda_{1}=0$,R=5) | 0.774$\pm$0.017 | 0.007$\pm$0.004 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.719$\pm$0.015 | 0.005$\pm$0.003 | 0.718$\pm$0.037 | 0.736$\pm$0.051
|  | M-SDWD ($\lambda_{1}=0$,R=1) | 0.704$\pm$0.018 | 0.007$\pm$0.006 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.702$\pm$0.020 | 0.010$\pm$0.008 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.732$\pm$0.015 | 0.007$\pm$0.004 | 0.262$\pm$0.016 | 0.993$\pm$0.002
| Less ($10\times 15\times 15$) | M-SDWD (R=5) | 0.965$\pm$0.006 | 0.000$\pm$0.000 | 0.965$\pm$0.014 | 0.374$\pm$0.062
|  | M-SDWD ($\lambda_{1}=0$,R=5) | 0.981$\pm$0.002 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.598$\pm$0.011 | 0.000$\pm$0.000 | 0.854$\pm$0.031 | 0.450$\pm$0.063
|  | M-SDWD ($\lambda_{1}=0$,R=1) | 0.603$\pm$0.010 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.574$\pm$0.013 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.813$\pm$0.005 | 0.000$\pm$0.000 | 0.420$\pm$0.024 | 0.906$\pm$0.020
| No ($30\times 15\times 15$) | M-SDWD (R=5) | 0.976$\pm$0.003 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-SDWD ($\lambda_{1}=0$,R=5) | 0.981$\pm$0.003 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-SDWD (R=1) | 0.573$\pm$0.010 | 0.000$\pm$0.000 | 0.794$\pm$0.037 | -
|  | M-SDWD ($\lambda_{1}=0$,R=1) | 0.590$\pm$0.008 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.557$\pm$0.012 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.828$\pm$0.006 | 0.000$\pm$0.000 | 0.556$\pm$0.037 | -
Low ($15\times 4\times 5$) | More ($5\times 2\times 2$) | M-SDWD (R=5 | 0.733$\pm$0.021 | 0.128$\pm$0.022 | 0.816$\pm$0.042 | 0.472$\pm$0.045
|  | M-SDWD ($\lambda_{1}=0$,R=5) | 0.664$\pm$0.022 | 0.120$\pm$0.017 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.755$\pm$0.022 | 0.093$\pm$0.015 | 0.671$\pm$0.047 | 0.703$\pm$0.053
|  | M-SDWD ($\lambda_{1}=0$,R=1) | 0.731$\pm$0.024 | 0.103$\pm$0.016 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.739$\pm$0.026 | 0.102$\pm$0.017 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.729$\pm$0.023 | 0.096$\pm$0.016 | 0.442$\pm$0.025 | 0.909$\pm$0.023
| Less ($5\times 4\times 5$) | M-SDWD (R=5) | 0.892$\pm$0.011 | 0.010$\pm$0.007 | 0.922$\pm$0.021 | 0.269$\pm$0.045
|  | M-SDWD ($\lambda_{1}=0$, R=5) | 0.889$\pm$0.010 | 0.005$\pm$0.004 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-SDWD (R=1) | 0.717$\pm$0.015 | 0.010$\pm$0.004 | 0.739$\pm$0.040 | 0.515$\pm$0.057
|  | M-SDWD ($\lambda_{1}=0$,R=1) | 0.722$\pm$0.014 | 0.011$\pm$0.004 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | M-DWD | 0.717$\pm$0.015 | 0.010$\pm$0.004 | 1.000$\pm$0.000 | 0.000$\pm$0.000
|  | Full SDWD | 0.795$\pm$0.010 | 0.008$\pm$0.004 | 0.468$\pm$0.031 | 0.841$\pm$0.029
| No ($15\times 4\times 5$) | M-SDWD (R=5) | 0.857$\pm$0.009 | 0.000$\pm$0.000 | 0.950$\pm$0.018 | -
|  | M-SDWD ($\lambda_{1}=0$, R=5) | 0.864$\pm$0.007 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-SDWD (R=1) | 0.693$\pm$0.013 | 0.000$\pm$0.000 | 0.836$\pm$0.033 | -
|  | M-SDWD ($\lambda_{1}=0$,R=1) | 0.701$\pm$0.012 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | M-DWD | 0.687$\pm$0.014 | 0.000$\pm$0.000 | 1.000$\pm$0.000 | -
|  | Full SDWD | 0.841$\pm$0.005 | 0.000$\pm$0.000 | 0.573$\pm$0.035 | -

Table 14:  Mean correlations between the estimated linear hyperplane and the true hyperplane by assumed rank for different true ranks for N=40 under high dimensional scenarios $30\times 15\times 15$ with more sparsity ($5\times 5\times 5$ of variables has signals have signal discriminating the classes).The assumed rank that return the largest correlation for each true rank is in bold. | $\hat{R}=1$ | $\hat{R}=2$ | $\hat{R}=3$ | $\hat{R}=4$ | $\hat{R}=5$
---|---|---|---|---|---
$R=1$ | 0.503 | 0.490 | 0.432 | 0.486 | 0.417
$R=2$ | 0.737 | 0.803 | 0.766 | 0.779 | 0.726
$R=3$ | 0.731 | 0.804 | 0.811 | 0.806 | 0.791
$R=4$ | 0.712 | 0.794 | 0.828 | 0.835 | 0.834
$R=5$ | 0.703 | 0.808 | 0.828 | 0.819 | 0.827

##  Appendix F Simulations with correlated predictors

In this section, we considered a simulation study in which the predictors are
correlated. The covariance matrix for the predictor array is assumed to have a
separable structure $\Sigma=\Sigma_{1}\otimes\Sigma_{2}\otimes\Sigma_{3}$,
where $\Sigma_{k}:P_{k}\times P_{k}$ is the covariance along dimension $k$.
Each of $\Sigma_{1}$, $\Sigma_{2}$, and $\Sigma_{3}$ have an AR(1) structure,
with correlation determined by a shared parameter $\rho$:

$\Sigma_{k}=\begin{bmatrix}1&\rho&\rho^{2}&\dots\\\ &1&\rho&\dots&\\\
&&\ddots&&\\\ &&&1\end{bmatrix}\;\;\text{ for }k=1,2,3.$

Thus, $\rho$ controls the overall level of correlation in the predictors. The
samples for class -1 are generated via
$\mbox{vec}(\mathbb{X}_{i})=\text{Normal}(\mathbf{0},\Sigma)$ and the samples
for class $+1$ are generated via
$\mbox{vec}(\mathbb{X}_{i})=\text{Normal}(\mbox{vec}(\sqrt{\alpha}\mu_{1}),\Sigma)$,
with $\sqrt{\alpha}\mu_{1}$ generated as in Section 5.1 of the main article.
As a representative scenario we consider the high-dimensional ($30\times
15\times 15$), more sparsity, small sample size (N=40) case of Section 5.1 in
the main article, and generate $\sqrt{\alpha}\mu_{1}$ under those conditions.
Table 15 shows the results over $200$ replications across different methods
and different correlation levels $\rho$. This shows that the M-SDWD method
performs relatively well for scenarios with no correlation or mild
correlation, but less well for scenarios with high correlation. We find that
the multiway approaches are more prone to convergence to a local optima with
high levels of correlation. Table 16 shows the results for those replications
that converge to a reasonable solution (with correlation greater than $0.5$
with the true vector), and the multiway approaches perform better under all
levels of residual correlation $\rho$ for these replications.

Table 15: Simulation results under the high dimensional ($30\times 15\times 15$), more sparsity, small sample size (N=40) with different residual correlation levels $\rho$. “Cor” is the correlation between the estimated linear hyperplane and the true hyperplane. “Mis” is the average misclassification rate. “TP” is the true positive rate, i.e., the proportion of non-zero coefficients that are correctly estimated to be non-zero.“TN” is the true negative rate, i.e., the proportion of zero coefficients that are correctly estimated to be zero. The margins of error (2* standard deviations across 200 replicates) for each statistic are also listed following the $\pm$ symbol. $\rho$ | Methods | Cor | Mis | TP | TN
---|---|---|---|---|---
0 | M-SDWD | 0.562$\pm$0.061 | 0.204$\pm$0.033 | 0.523$\pm$0.049 | 0.724$\pm$0.051
| M-SDWD ($\lambda_{1}=0$) | 0.519$\pm$0.058 | 0.205$\pm$0.031 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| M-DWD | 0.516$\pm$0.060 | 0.226$\pm$0.034 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| Full SDWD | 0.347$\pm$0.037 | 0.263$\pm$0.027 | 0.197$\pm$0.033 | 0.884$\pm$0.034
0.3 | M-SDWD | 0.454$\pm$0.062 | 0.255$\pm$0.033 | 0.576$\pm$0.052 | 0.649$\pm$0.058
| M-SDWD ($\lambda_{1}=0$) | 0.436$\pm$0.059 | 0.260$\pm$0.032 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| M-DWD | 0.310$\pm$0.058 | 0.337$\pm$0.032 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| Full SDWD | 0.355$\pm$0.036 | 0.253$\pm$0.026 | 0.193$\pm$0.032 | 0.884$\pm$0.034
0.6 | M-SDWD | 0.323$\pm$0.062 | 0.325$\pm$0.035 | 0.456$\pm$0.057 | 0.680$\pm$0.058
| M-SDWD ($\lambda_{1}=0$) | 0.292$\pm$0.060 | 0.356$\pm$0.033 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| M-DWD | 0.067$\pm$0.032 | 0.474$\pm$0.019 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| Full SDWD | 0.375$\pm$0.039 | 0.267$\pm$0.029 | 0.177$\pm$0.034 | 0.900$\pm$0.035
0.9 | M-SDWD | 0.434$\pm$0.069 | 0.256$\pm$0.040 | 0.393$\pm$0.060 | 0.791$\pm$0.054
| M-SDWD ($\lambda_{1}=0$) | 0.427$\pm$0.072 | 0.288$\pm$0.040 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| M-DWD | 0.044$\pm$0.027 | 0.477$\pm$0.017 | 1.000$\pm$0.000 | 0.000$\pm$0.000
| Full SDWD | 0.561$\pm$0.043 | 0.177$\pm$0.032 | 0.193$\pm$0.040 | 0.897$\pm$0.043
Table 16: Simulation results among replications with “Cor”$>$0.5 under the
high dimensional ($30\times 15\times 15$), more sparsity, small sample size
(N=40) with different residual correlation levels $\rho$. “Cor” is the
correlation between the estimated linear hyperplane and the true hyperplane.
“Mis” is the average misclassification rate. “TP” is the true positive rate,
i.e., the proportion of non-zero coefficients that are correctly estimated to
be non-zero.“TN” is the true negative rate, i.e., the proportion of zero
coefficients that are correctly estimated to be zero. The margins of error (2*
standard deviations across 200 replicates) for each statistic are also listed
following the $\pm$ symbol.

$\rho$ | Methods | Cor | Mis | TP | TN |  $\%$Cor$>0.5$
---|---|---|---|---|---|---
0 | M-SDWD | 0.897$\pm$0.014 | 0.027$\pm$0.007 | 0.582$\pm$0.050 | 0.832$\pm$0.051 | 0.62
| M-SDWD ($\lambda_{1}=0$) | 0.841$\pm$0.018 | 0.036$\pm$0.010 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.60
| M-DWD | 0.864$\pm$0.018 | 0.036$\pm$0.010 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.59
| Full SDWD | 0.673$\pm$0.029 | 0.043$\pm$0.015 | 0.131$\pm$0.018 | 0.993$\pm$0.002 | 0.32
0.3 | M-SDWD | 0.901$\pm$0.016 | 0.022$\pm$0.008 | 0.601$\pm$0.058 | 0.840$\pm$0.056 | 0.49
| M-SDWD ($\lambda_{1}=0$) | 0.850$\pm$0.021 | 0.036$\pm$0.011 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.49
| M-DWD | 0.893$\pm$0.022 | 0.028$\pm$0.011 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.33
| Full SDWD | 0.674$\pm$0.028 | 0.042$\pm$0.016 | 0.121$\pm$0.015 | 0.994$\pm$0.001 | 0.32
0.6 | M-SDWD | 0.903$\pm$0.022 | 0.009$\pm$0.008 | 0.468$\pm$0.079 | 0.912$\pm$0.057 | 0.35
| M-SDWD ($\lambda_{1}=0$) | 0.906$\pm$0.023 | 0.026$\pm$0.015 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.30
| M-DWD | 0.956$\pm$0.032 | 0.012$\pm$0.024 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.05
| Full SDWD | 0.679$\pm$0.028 | 0.050$\pm$0.016 | 0.123$\pm$0.018 | 0.995$\pm$0.001 | 0.36
0.9 | M-SDWD | 0.877$\pm$0.026 | 0.002$\pm$0.002 | 0.475$\pm$0.075 | 0.934$\pm$0.042 | 0.49
| M-SDWD ($\lambda_{1}=0$) | 0.935$\pm$0.019 | 0.018$\pm$0.013 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.45
| M-DWD | 0.885$\pm$0.131 | 0.083$\pm$0.110 | 1.000$\pm$0.000 | 0.000$\pm$0.000 | 0.04
| Full SDWD | 0.730$\pm$0.018 | 0.056$\pm$0.016 | 0.117$\pm$0.011 | 0.996$\pm$0.001 | 0.69

##  Appendix G Tuning parameters selection for MRS data application

Here we expand on the selection of the tuning parameters for the application
of multi-way sparse DWD to classify dox treated mice (dox group) and controls
(no dox group) in Section 6.1 of the main article. We computed the predicted
DWD scores for all subjects by 10-fold cross-validation for a grid of
$\lambda_{1}$ ($0,0.0001,0.001,0.005,0.01,0.025,0.05)$ and
$\lambda_{2}$($0.25,0.5,0.75,1,3,5$), and select the best combinations of
parameters with the maximum t-test statistics that testing the differences
between the scores and the two classes. Table 17 shows the t-test statistics
between dox and no-dox groups for a grid of $\lambda_{1}$ and $\lambda_{2}$.
The best pair of tuning parameters are $\lambda_{1}=0.01$ and $\lambda_{2}=5$
with the maximum test statistic. Table 18 shows the misclassification rates
that classify the TG-dox group and other mice. Using t-test statistic as index
for selecting parameters gives a unique optimal pair of parameters, which also
located in the region of parameters with minimal misclassification rates.

Table 17: T-test statistics that test the differences of predicted DWD scores between dox and no-dox groups for a grid of $\lambda_{1}$ and $\lambda_{2}$ $\lambda_{2}$/$\lambda_{1}$ | 0 | 1e-04 | 0.001 | 0.005 | 0.01 | 0.025 | 0.05 | 0.1 | 0.25
---|---|---|---|---|---|---|---|---|---
0.25 | 2.098 | 2.105 | 2.194 | 2.537 | 2.789 | 2.759 | 2.823 | 3.189 | 3.549
0.50 | 1.904 | 1.906 | 1.932 | 2.060 | 2.303 | 2.712 | 2.778 | 2.948 | 3.090
0.75 | 1.556 | 1.564 | 1.644 | 1.950 | 2.350 | 2.678 | 2.663 | 2.658 | 3.020
1.00 | 2.676 | 2.677 | 2.684 | 2.694 | 2.762 | 3.032 | 3.024 | 3.069 | 3.455
3.00 | 3.251 | 3.253 | 3.266 | 3.308 | 3.377 | 3.483 | 3.525 | 3.573 | 3.660
5.00 | 3.172 | 3.172 | 3.180 | 3.216 | 3.264 | 3.327 | 3.384 | 3.511 | 3.427
Table 18: Misclassification rates that classify TG-dox and other groups for a grid of $\lambda_{1}$ and $\lambda_{2}$ $\lambda_{2}$/$\lambda_{1}$ | 0 | 1e-04 | 0.001 | 0.005 | 0.01 | 0.025 | 0.05 | 0.1 | 0.25
---|---|---|---|---|---|---|---|---|---
0.25 | 0.333 | 0.333 | 0.333 | 0.333 | 0.286 | 0.286 | 0.286 | 0.286 | 0.190
0.50 | 0.286 | 0.286 | 0.286 | 0.286 | 0.286 | 0.143 | 0.238 | 0.143 | 0.190
0.75 | 0.333 | 0.333 | 0.333 | 0.286 | 0.286 | 0.238 | 0.238 | 0.190 | 0.143
1.00 | 0.143 | 0.143 | 0.095 | 0.143 | 0.095 | 0.095 | 0.095 | 0.095 | 0.095
3.00 | 0.048 | 0.048 | 0.048 | 0.048 | 0.048 | 0.048 | 0.048 | 0.095 | 0.048
5.00 | 0.095 | 0.095 | 0.095 | 0.048 | 0.048 | 0.048 | 0.048 | 0.095 | 0.095

## References

  * (1)
  * Baranzini et al. (2004) Baranzini, S. E., Mousavi, P., Rio, J., Caillier, S. J., Stillman, A., Villoslada, P., Wyatt, M. M., Comabella, M., Greller, L. D., Somogyi, R. et al. (2004), ‘Transcription-based prediction of response to ifn$\beta$ using supervised computational methods’, PLoS Biol 3(1), e2.
  * Bauckhage (2007) Bauckhage, C. (2007), Robust tensor classifiers for color object recognition, in ‘International Conference Image Analysis and Recognition’, Springer, pp. 352–363.
  * Bradley & Mangasarian (1998) Bradley, P. S. & Mangasarian, O. L. (1998), Feature selection via concave minimization and support vector machines., in ‘ICML’, Vol. 98, pp. 82–90.
  * Breiman (2001) Breiman, L. (2001), ‘Random forests’, Machine learning 45(1), 5–32.
  * Cai & Liu (2011) Cai, T. & Liu, W. (2011), ‘A direct estimation approach to sparse linear discriminant analysis’, Journal of the American statistical association 106(496), 1566–1577.
  * Cortes & Vapnik (1995) Cortes, C. & Vapnik, V. (1995), ‘Support-vector networks’, Machine learning 20(3), 273–297.
  * Cover & Hart (1967) Cover, T. & Hart, P. (1967), ‘Nearest neighbor pattern classification’, IEEE transactions on information theory 13(1), 21–27.
  * Fan & Fan (2008) Fan, J. & Fan, Y. (2008), ‘High dimensional classification using features annealed independence rules’, Annals of statistics 36(6), 2605.
  * Fan et al. (2012) Fan, J., Feng, Y. & Tong, X. (2012), ‘A road to classification in high dimensional space: the regularized optimal affine discriminant’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 74(4), 745–771.
  * França et al. (2009) França, M. C., D’Abreu, A., Yasuda, C. L., Bonadia, L. C., Da Silva, M. S., Nucci, A., Lopes-Cendes, I. & Cendes, F. (2009), ‘A combined voxel-based morphometry and 1 h-mrs study in patients with friedreich’s ataxia’, Journal of neurology 256(7), 1114–1120.
  * Gramegna et al. (2017) Gramegna, L. L., Tonon, C., Manners, D. N., Pini, A., Rinaldi, R., Zanigni, S., Bianchini, C., Evangelisti, S., Fortuna, F., Carelli, V. et al. (2017), ‘Combined cerebellar proton mr spectroscopy and dwi study of patients with friedreich’s ataxia’, The Cerebellum 16(1), 82–88.
  * Harshman (1970) Harshman, R. A. (1970), ‘Foundations of the PARAFAC procedure: Models and conditions for an “explanatory” multi-modal factor analysis’, UCLA Working Papers in Phonetics 16, 1–84.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R. & Friedman, J. (2009), The elements of statistical learning: data mining, inference, and prediction, Springer Science & Business Media.
  * Huang et al. (2013) Huang, H., Liu, Y., Du, Y., Perou, C. M., Hayes, D. N., Todd, M. J. & Marron, J. S. (2013), ‘Multiclass distance-weighted discrimination’, Journal of Computational and Graphical Statistics 22(4), 953–969.
  * Hunter & Lange (2004) Hunter, D. R. & Lange, K. (2004), ‘A tutorial on mm algorithms’, The American Statistician 58(1), 30–37.
  * Iltis et al. (2010) Iltis, I., Hutter, D., Bushara, K. O., Clark, H. B., Gross, M., Eberly, L. E., Gomez, C. M. & Öz, G. (2010), ‘1h mr spectroscopy in friedreich’s ataxia and ataxia with oculomotor apraxia type 2’, Brain research 1358, 200–210.
  * Liu et al. (2011) Liu, Y., Zhang, H. H. & Wu, Y. (2011), ‘Hard or soft classification? large-margin unified machines’, Journal of the American Statistical Association 106(493), 166–177.
  * Lock (2018) Lock, E. F. (2018), ‘Tensor-on-tensor regression’, Journal of Computational and Graphical Statistics 27(3), 638–647.
  * Lyu et al. (2017) Lyu, T., Lock, E. F. & Eberly, L. E. (2017), ‘Discriminating sample groups with multi-way data’, Biostatistics 18(3), 434–450.
  * Mai et al. (2012) Mai, Q., Zou, H. & Yuan, M. (2012), ‘A direct approach to sparse discriminant analysis in ultra-high dimensions’, Biometrika 99(1), 29–42.
  * Marron et al. (2007) Marron, J. S., Todd, M. J. & Ahn, J. (2007), ‘Distance-weighted discrimination’, Journal of the American Statistical Association 102(480), 1267–1271.
  * Pan et al. (2018) Pan, Y., Mai, Q. & Zhang, X. (2018), ‘Covariate-adjusted tensor classification in high dimensions’, Journal of the American statistical association .
  * Pandolfo (2008) Pandolfo, M. (2008), ‘Friedreich Ataxia’, Archives of Neurology 65(10), 1296–1303.
  * Tao et al. (2005) Tao, D., Li, X., Hu, W., Maybank, S. & Wu, X. (2005), Supervised tensor learning, in ‘Fifth IEEE International Conference on Data Mining (ICDM’05)’, IEEE, pp. 8–pp.
  * Wang & Zou (2015) Wang, B. & Zou, H. (2015), sdwd: Sparse Distance Weighted Discrimination.  R package version 1.0.2.
  * Wang & Zou (2016) Wang, B. & Zou, H. (2016), ‘Sparse distance weighted discrimination’, Journal of Computational and Graphical Statistics 25(3), 826–838.
  * Wang et al. (2008) Wang, L., Zhu, J. & Zou, H. (2008), ‘Hybrid huberized support vector machines for microarray classification and gene selection’, Bioinformatics 24(3), 412–419.
  * Wimalawarne et al. (2016) Wimalawarne, K., Tomioka, R. & Sugiyama, M. (2016), ‘Theoretical and experimental analyses of tensor-based regression and classification’, Neural computation 28(4), 686–715.
  * Witten & Tibshirani (2011) Witten, D. M. & Tibshirani, R. (2011), ‘Penalized classification using Fisher’s linear discriminant’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73(5), 753–772.
  * Ye et al. (2004) Ye, J., Janardan, R. & Li, Q. (2004), ‘Two-dimensional linear discriminant analysis’, Advances in neural information processing systems 17, 1569–1576.
  * Zacharoff et al. (2012) Zacharoff, L., Tkac, I., Song, Q., Tang, C., Bolan, P. J., Mangia, S., Henry, P.-G., Li, T. & Dubinsky, J. M. (2012), ‘Cortical metabolites as biomarkers in the r6/2 model of huntington’s disease’, Journal of Cerebral Blood Flow & Metabolism 32(3), 502–514.
  * Zhang et al. (2006) Zhang, H. H., Ahn, J., Lin, X. & Park, C. (2006), ‘Gene selection using support vector machines with non-convex penalty’, Bioinformatics 22(1), 88–95.
  * Zou (2019) Zou, H. (2019), ‘Classification with high dimensional features’, Wiley Interdisciplinary Reviews: Computational Statistics 11(1), e1453.