Sturm: Sparse Tubal-Regularized Multilinear Regression for fMRI

TL;DR

This paper introduces Sturm, a novel sparse multilinear regression method utilizing t-SVD and tubal tensor nuclear norm for improved classification of high-dimensional fMRI data with limited samples.

Contribution

The paper proposes a new tensor regression approach based on t-SVD that enhances sparsity and classification accuracy in fMRI analysis, overcoming limitations of previous tensor methods.

Findings

Sturm outperforms existing methods in fMRI classification tasks.

Effective with a small number of voxels.

Demonstrates superior accuracy on multiple datasets.

Abstract

While functional magnetic resonance imaging (fMRI) is important for healthcare/neuroscience applications, it is challenging to classify or interpret due to its multi-dimensional structure, high dimensionality, and small number of samples available. Recent sparse multilinear regression methods based on tensor are emerging as promising solutions for fMRI, yet existing works rely on unfolding/folding operations and a tensor rank relaxation with limited tightness. The newly proposed tensor singular value decomposition (t-SVD) sheds light on new directions. In this work, we study t-SVD for sparse multilinear regression and propose a Sparse tubal-regularized multilinear regression (Sturm) method for fMRI. Specifically, the Sturm model performs multilinear regression with two regularization terms: a tubal tensor nuclear norm based on t-SVD and a standard L1 norm. We further derive the algorithm under the alternating direction method of multipliers framework. We perform experiments on four classification problems, including both resting-state fMRI for disease diagnosis and task-based fMRI for neural decoding. The results show the superior performance of Sturm in classifying fMRI using just a small number of voxels.

Figures

Tables

Table 1. Table 1 : Important notations.

Notation	Description
$a$	Lowercase letter denotes scalar
$𝐚$	Bold lowercase letter denotes vector
$𝐀$	Bold uppercase letter denotes matrix
$𝒜$	Calligraphic uppercase letter denotes tensor
$𝒜\ast ℬ$	t-product between tensors $𝒜$ and $ℬ$
${𝒜}^{\top }$	Tensor conjugate transpose of $𝒜$
${𝐀}^{\left(i_3\right)}$	The $i_3$th mode-3 (frontal) slice of $𝒜$
${𝒜}_{ℱ}$	The discrete Fourier transform of $𝒜$

Table 2. Table 2 : Summary of the four classification problems and respective datasets. # denotes the number of volumes (samples). Pos. and Neg. are short for the positive and negative classes, respectively. For diagnosis problems on ABIDE/ADHD-200, patients and health subjects are considered as positive and negative classes, respectively. The two neural decoding problems are formed by using six OpenfMRI datasets listed as Pos. vs Neg.

Classification Problem	# Pos.	# Neg.	Input data size $I_1\times I_2\times I_3$
Resting 1 – ABIDENYU&UM	101	131	$61\times 73\times 61$
Resting 2 – ADHD-200NYU	118	98	$53\times 64\times 46$
Task 1 – Balloon vs Mixed	32	32	$91\times 109\times 91$
Task 2 – Simon vs Flanker	42	52	$91\times 109\times 91$

Table 3. Table 3 : Classification accuracy (mean ± plus-or-minus \pm standard deviation in %). Resting 1 and Resting 2 denote two disease diagnosis problems on ABIDE NYU&UM and ADHD-200, respectively. Task 1 and Task 2 denote two neural decoding problems on OpenfMRI datasets for Balloon vs Mixed gamble and Simon vs Flanker, respectively. The best accuracy among all of the compared algorithms for each column is highlighted in bold and the second best is underlined .

Method	Resting 1	Resting 2	Task 1	Task 2	Average
					Resting	Task	All
SVM	60.78 $\pm$ 0.09	63.97 $\pm$ 0.09	87.38 $\pm$ 0.12	82.56 $\pm$ 0.17	62.38	84.97	73.67
Lasso	61.16 $\pm$ 0.08	64.84 $\pm$ 0.11	87.38 $\pm$ 0.12	85.22 $\pm$ 0.07	63.00	86.30	74.65
ENet	61.21 $\pm$ 0.10	64.38 $\pm$ 0.10	81.19 $\pm$ 0.15	82.56 $\pm$ 0.17	62.80	81.87	72.34
Remurs	60.72 $\pm$ 0.08	62.13 $\pm$ 0.09	87.14 $\pm$ 0.13	84.67 $\pm$ 0.15	61.43	85.90	73.67
Sturm	62.05 $\pm$ 0.11	63.47 $\pm$ 0.07	89.10 $\pm$ 0.09	86.89 $\pm$ 0.16	62.76	88.00	75.38
Lasso + SVM	63.37 $\pm$ 0.08	62.56 $\pm$ 0.09	74.05 $\pm$ 0.20	72.11 $\pm$ 0.16	62.97	73.08	68.02
ENet + SVM	64.20 $\pm$ 0.07	61.61 $\pm$ 0.08	76.43 $\pm$ 0.14	72.00 $\pm$ 0.14	62.91	74.21	68.56
Remurs + SVM	64.67 $\pm$ 0.10	60.23 $\pm$ 0.10	81.19 $\pm$ 0.12	83.56 $\pm$ 0.19	62.45	82.37	72.41
Sturm+ SVM	64.66 $\pm$ 0.12	66.24 $\pm$ 0.06	78.10 $\pm$ 0.22	82.44 $\pm$ 0.16	65.45	80.27	72.86

Table 4. Table 4 : Sparsity (mean ± plus-or-minus \pm standard deviation) for respective results in Table 3 with the best and second best highlighted.

Method	Resting 1	Resting 2	Task 1	Task 2	Average
					Resting	Task	All
Lasso	0.52 $\pm$ 0.09	0.23 $\pm$ 0.32	0.74 $\pm$ 0.12	0.73 $\pm$ 0.01	0.38	0.73	0.55
ENet	0.60 $\pm$ 0.01	0.01 $\pm$ 0.01	0.96 $\pm$ 0.05	0.95 $\pm$ 0.03	0.31	0.96	0.63
Remurs	0.69 $\pm$ 0.03	0.73 $\pm$ 0.17	0.81 $\pm$ 0.08	0.81 $\pm$ 0.07	0.71	0.81	0.76
Sturm	0.86 $\pm$ 0.18	0.86 $\pm$ 0.24	0.72 $\pm$ 0.24	0.60 $\pm$ 0.15	0.86	0.66	0.76
Lasso + SVM	0.57 $\pm$ 0.05	0.19 $\pm$ 0.40	0.77 $\pm$ 0.10	0.75 $\pm$ 0.06	0.38	0.76	0.57
ENet + SVM	0.58 $\pm$ 0.09	0.02 $\pm$ 0.01	0.96 $\pm$ 0.04	0.95 $\pm$ 0.04	0.30	0.96	0.63
Remurs + SVM	0.70 $\pm$ 0.13	0.74 $\pm$ 0.17	0.80 $\pm$ 0.04	0.79 $\pm$ 0.13	0.72	0.79	0.76
Sturm + SVM	0.87 $\pm$ 0.07	0.99 $\pm$ 0.01	0.85 $\pm$ 0.14	0.56 $\pm$ 0.11	0.93	0.71	0.82

Equations

$${\mathring{\mathbf{c}}}_{i_1{j}_4}=\mathcal{C}(i_1,j_4,:)=\sum _{i_2=1}^{I_2}\mathcal{A}(i_1,i_2,:)\ast \mathcal{B}(i_2,j_4,:),$$

$${\mathcal{Q}}^{\top }\ast \mathcal{Q}=\mathcal{Q}\ast {\mathcal{Q}}^{\top }=\mathcal{I},$$

$$\mathcal{A}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^{\top },$$

$$\Vert \mathcal{A}{\Vert }_{TNN}=\frac{1}{I_3}\sum _{i_3=1}^{I_3}\Vert {\mathbf{A}}_{\mathcal{F}}^{(i_3)}{\Vert }_{\ast },$$

$$\underset{\mathcal{W}}{\min }\frac{1}{M}\sum _{m=1}^{M}L(⟨{\mathcal{X}}_m,\mathcal{W}⟩,y_m)+\lambda \Omega (\mathcal{W}),$$

$$⟨\mathcal{X},\mathcal{W}⟩:=\sum _{i_1}\sum _{i_2}\sum _{i_3}\mathcal{X}(i_1,i_2,i_3)\cdot \mathcal{W}(i_1,i_2,i_3),$$

$$\underset{\mathcal{W}}{\min }\frac{1}{2}\sum _{m=1}^M(y_m-⟨{\mathcal{X}}_m,\mathcal{W}⟩{)}^2+\tau \Vert \mathcal{W}{\Vert }_{TNN}+\gamma \Vert \mathcal{W}{\Vert }_1,$$

$$\Vert \mathcal{W}{\Vert }_1=\sum _{i_1}\sum _{i_2}\sum _{i_3}\left|\mathcal{W}(i_1,i_2,i_3)\right|,$$

$$\underset{\mathcal{W}}{\min }\frac{1}{2}\sum _{m=1}^M(y_m-⟨{\mathcal{X}}_m,\mathcal{A}⟩{)}^2+\tau \Vert \mathcal{B}{\Vert }_{TNN}+\gamma \Vert \mathcal{W}{\Vert }_1$$

$$s.t.\mathcal{A}=\mathcal{W}\text{and }\mathcal{B}=\mathcal{W}.$$

$$\begin{split}&L_{\rho}(\mathcal{A},\mathcal{B},\mathcal{W},\mathcal{P},\mathcal{Q})=\frac{1}{2}\sum_{m=1}^{M}(y_{m}-\langle\mathcal{X}_{m},\mathcal{A}\rangle)^{2}\\ &+\tau\|\mathcal{B}\|_{TNN}+\gamma\|\mathcal{W}\|_{1}\\ &+\Big{\langle}\mathcal{P},\mathcal{A}-\mathcal{W}\Big{\rangle}+\frac{\rho}{2}\|\mathcal{A}-\mathcal{W}\|_{F}^{2}\\ &+\Big{\langle}\mathcal{Q},\mathcal{B}-\mathcal{W}\Big{\rangle}+\frac{\rho}{2}\|\mathcal{B}-\mathcal{W}\|_{F}^{2}.\end{split}$$

$$\begin{array}{rl}& {\mathcal{A}}^{k+1}=\arg \underset{\mathcal{A}}{\min }{L}_{\rho }(\mathcal{A},{\mathcal{B}}^k,{\mathcal{W}}^k,{{\mathcal{P}}^{\prime }}^k,{{\mathcal{Q}}^{\prime }}^k)\\ & =\arg \underset{\mathcal{A}}{\min }\frac{1}{2}\sum _{m=1}^M(y_m-⟨{\mathcal{X}}_m,\mathcal{A}⟩{)}^2+\frac{\rho }{2}\Vert \mathcal{A}-{\mathcal{W}}^k+{{\mathcal{P}}^{\prime }}^k{\Vert }_F.\end{array}$$

$${\mathbf{a}}^{k+1}=({\mathbf{X}}^{\top }\mathbf{X}+\rho \mathbf{I}{)}^{-1}({\mathbf{X}}^{\top }\mathbf{y}+\rho ({\mathbf{w}}^k-{{\mathbf{p}}^{\prime }}^k)),$$

$$\begin{array}{rl}& {\mathcal{B}}^{k+1}=\arg \underset{\mathcal{B}}{\min }{L}_{\rho }({\mathcal{A}}^{k+1},\mathcal{B},{\mathcal{W}}^k,{{\mathcal{P}}^{\prime }}^k,{{\mathcal{Q}}^{\prime }}^k)\\ & =\arg \underset{\mathcal{B}}{\min }\tau \Vert \mathcal{B}{\Vert }_{TNN}+\frac{\rho }{2}\Vert \mathcal{B}-{\mathcal{W}}^k+{{\mathcal{Q}}^{\prime }}^k{\Vert }_F^2\\ & ={\mathtt{prox}}_{\frac{\tau }{\rho }\Vert \cdot {\Vert }_{TNN}}({\mathcal{W}}^k-{{\mathcal{Q}}^{\prime }}^k).\end{array}$$

$${\mathtt{prox}}_{\mu \Vert \cdot {\Vert }_{TNN}}(\mathcal{T}):=\arg \underset{\mathcal{W}}{\min }\mu \Vert \mathcal{W}{\Vert }_{TNN}+\frac{1}{2}\Vert \mathcal{W}-\mathcal{T}{\Vert }_F^2,$$

$$\begin{split}&\mathcal{W}^{k+1}=\arg\min_{\mathcal{W}}L_{\rho}(\mathcal{A}^{k+1},\mathcal{B}^{k+1},\mathcal{W},\mathcal{P^{\prime}}^{k},\mathcal{Q^{\prime}}^{k})\\ &={\tt prox}_{\frac{\gamma}{2\rho}\|\cdot\|_{1}}\Big{(}\frac{\mathcal{A}^{k+1}+\mathcal{P^{\prime}}^{k}+\mathcal{B}^{k+1}+\mathcal{Q^{\prime}}^{k}}{2}\Big{)}.\end{split}$$

$${\mathtt{prox}}_{\mu \Vert \cdot {\Vert }_1}(\mathcal{T})=(\mathcal{T}-\mu {)}_{+}.$$

$${{\mathcal{P}}^{\prime }}^{k+1}={{\mathcal{P}}^{\prime }}^k+{\mathcal{A}}^{k+1}-{\mathcal{W}}^{k+1},$$

$${{\mathcal{Q}}^{\prime }}^{k+1}={{\mathcal{Q}}^{\prime }}^k+{\mathcal{B}}^{k+1}-{\mathcal{W}}^{k+1}.$$

$${\tt prox}_{\frac{\gamma}{2\rho}\|\cdot\|_{1}}\Big{(}\frac{\mathcal{A}^{k+1}+\mathcal{P^{\prime}}^{k}+\mathcal{B}^{k+1}+\mathcal{Q^{\prime}}^{k}}{2}\Big{)};$$

Full text

1 Introduction

Brain diseases affect millions of people worldwide and impose significant challenges to healthcare systems. Functional magnetic resonance imaging (fMRI) is a key medical imaging technique for diagnosis, monitoring and treatment of brain diseases. Beyond healthcare, fMRI is also an indispensable tool in neuroscience studies (Faro & Mohamed, 2010).

fMRI records the Blood Oxygenation Level Dependent (BOLD) signal caused by changes in blood flow (Ogawa et al., 1990), as depicted in Fig. 1. fMRI scan is a four-dimensional (4-D) sequence composed of magnetic resonance imaging (MRI) volumes sampled every few seconds, with over 100,000 voxels each MRI volume. It is often represented as a three-dimensonal (3-D) volume, with the statistics of each voxel summarized along the time axis. Figure 1 visualizes example values of such statistics, showing very different characteristics from natural images.

Unlike natural images, fMRI data is expensive to obtain. Thus, the number of fMRI samples in a study is typically limited to dozens. This makes fMRI challenging to analyze, particularly in its full (i.e., whole brain). Moreover, in healthcare and neuroscience, prediction accuracy is not the only concern. It is also important to interpret the learned features to domain experts such as clinicians or neuroscientists. This makes sparse learning models (Ryali et al., 2010; Simon et al., 2013; Rao et al., 2013; Hastie et al., 2015) attractive because they can reveal direct dependency of a response with a small portion of input features. Therefore, tensor-based, sparse multilinear regression methods are emerging recently, where tensor refers to multidimensional array. For simplicity, we consider only third-order (i.e., 3-D) tensor in this paper.

Sparse multilinear regression models relate a predictor/feature tensor with a univariate response via a coefficient tensor, generalizing Lasso-based models (Tibshirani, 1996; Hastie et al., 2015) to tensor data. Regularization that promotes sparsity and low rankness is also generalized to the coefficient tensor. For example, the regularized multilinear regression and selection (Remurs) model (Song & Lu, 2017) incorporates a sparse regularization term, via an $\ell_{1}$ norm, and a Tucker rank-minimization term, via a summation of the nuclear norms (SNN) of unfolded matrices. There are also CANDECOMP/PARAFAC (CP) rank-based methods (Tan et al., 2012; He et al., 2018). For example, the fast stagewise unit-rank tensor factorization (SURF) (He et al., 2018) enforces the low rankness via the CP decomposition model (Hitchcock, 1927; Harshman, 1970; Carroll & Chang, 1970), with the rank to be incrementally estimated via residual computation rather than direct minimization. And it proposes a divide-and-conquer strategy to improve efficiency and scalability with a greedy algorithm (Cormen et al., 2009). Although Remurs minimizes the Tucker rank directly, it involves unfolding/folding operations that transform tensor to/from matrix, which break some tensor structure, and SNN is not a tight convex relaxation of its target rank (Romera-Paredes & Pontil, 2013).

A new Tubal Tensor Nuclear Norm (TNN) (Zhang et al., 2014; Zhang & Aeron, 2017; Lu et al., 2016, 2018a) is recently proposed based on the tubal rank, which originates from the tensor singular value decomposition (t-SVD) (Kilmer et al., 2013). In tensor recovery problems, TNN guarantees the exact recovery under certain conditions. Improved performance has been observed on image/video data on tensor completion/recovery and robust tensor principal component analysis (PCA) (Zhang & Aeron, 2017; Lu et al., 2016, 2018a). In addition, TNN does not require unfolding/folding of tensor in its optimization.

In this work, we study sparse multilinear regression under the t-SVD framework for fMRI classification. The success of TNN is limited to unsupervised learning settings such as completion/recovery and robust PCA. To our knowledge, TNN has not been studied in a supervised setting yet, such as multilinear regression where the problem is not recovery given samples of a tensor but prediction of a response with a set of training tensor samples. Moreover, the targeted fMRI classification tasks have additional challenge of small sample size (relative to the feature dimension).

Our contributions are twofold: Firstly, we propose a Sparse tubal-regularized multilinear regression (Sturm) method that incorporates TNN regularization. Specifically, we formulate the Sturm model by incorporating both a TNN regularization and a sparsity regularization on the coefficient tensor. We solve the resulted Sturm problem using the alternating direction method of multipliers (ADMM) framework (Boyd et al., 2011). TNN-based formulation allows efficient parameter update in the Fourier domain, which is highly parallelizable.

Our second contribution is that we evaluate Sturm and related methods on both resting-state and task-based fMRI classification problems, instead of only one of them as in previous works (Zhou et al., 2013b, 2014; Shi et al., 2014; Zhou et al., 2017; He et al., 2017). We use public datasets with identifiable subsets for repeatibility, and examine both the classification accuracy and sparsity. The results show Sturm outperforming other state-of-the-art methods on the whole.

2 Related Work

Tensor decomposition and rank. There are three popular tensor decomposition methods and associated tensor rank definitions: 1) The CP decomposition of a tensor is written as the summation of $R$ rank-one tensors (Hitchcock, 1927; Harshman, 1970; Carroll & Chang, 1970), where $R$ is the CP rank. 2) For a 3-D tensor, the Tucker decomposition decomposes it into a core tensor and three factor matrices (De Lathauwer et al., 2000), and the Tucker rank is 3-tuple consisting of the rank for each mode-$n$ unfolded matrix. 3) The t-SVD views a 3-D tensor as a matrix of tubes (mode-$3$ vectors) oriented along the third dimension and decomposes it as circular convolutions on three tensors (Braman, 2010; Kilmer et al., 2013; Kilmer & Martin, 2011; Gleich et al., 2013). This leads to a new tensor tubal rank defined as the number of non-zero singular tubes. Please refer to (Zhang & Aeron, 2017; Zhang et al., 2014; Yuan & Zhang, 2016; Lu et al., 2016, 2018a; Mu et al., 2014) for more detailed discussion, e.g., on their pros and cons.

Low-rank tensor completion/recovery. Tensor completion/recovery takes a tensor with missing/noisy entries as the input and aims to recover/complete those entries. Low-rank assumption makes solving such problems feasible, and often with provable theoretical guarantees under mild conditions. The three types of tensor decomposition have their corresponding tensor completion/recovery approaches that minimize respective tensor ranks. Direct minimization of tensor rank is NP-hard (Hillar & Lim, 2013). Therefore, CP rank, Tucker rank, and tubal rank are relaxed to CP-based nuclear norm (Shi et al., 2017), the sum of matrix nuclear norm (Liu et al., 2013), and the tubal tensor nuclear norm.

Sparse and low-rank multilinear regression. Sparse and low-rank constraints have also been applied to multilinear regression problems. Multilinear regression models (Signoretto et al., 2010; Su et al., 2012; Guo et al., 2012; Zhou et al., 2013a) relate a predictor/feature tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ with a univariate response $y$, via a coefficient tensor $\mathcal{W}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$. The Remurs model (Song & Lu, 2017) introduces regularization with a sparsity constraint, via $\ell_{1}$ norm, and a low-rank constraint via Tucker rank-based SNN. Instead of minimizing the rank, the fast stagewise unit-rank tensor factorization (SURF) (He et al., 2018) is an efficient and scalable method that imposes a low CP rank constraint by setting a maximum rank and increasing the rank estimate from 1 to the maximum, stopping upon zero residual.

Tensor methods on fMRI. Due to the high dimensionality of fMRI data, regions of interest (ROIs) are typically used rather than all voxels in the original 3-D spatial domain (Chen et al., 2015; He et al., 2018). ROI analysis requires strong prior knowledge to determine the regions. In contrast, whole-brain fMRI analysis (Poldrack et al., 2013) is more data-driven. Thus, tensor-based machine learning methods (Cichocki, 2013; Cichocki et al., 2009) have been developed for fMRI, with promising results reported (Acar et al., 2017; Zhou et al., 2013a; Song et al., 2015; He et al., 2017; Song & Lu, 2017; Ozdemir et al., 2017; Barnathan et al., 2011). However, in these works, the learning methods are only evaluated on either resting-state fMRI for disease diagnosis (Zhou et al., 2017) or task-based fMRI for neural decoding (He et al., 2017; Chen et al., 2015; Song & Lu, 2017), but not both.

3 Tubal Tensor Nuclear Norm

Notations. Table 1 summarizes important notations used in this paper. We use lowercase, bold lowercase, bold uppercase, calligraphic uppercase letters to denote scalar, vector, matrix, and tensor, respectively. We denote indices by lowercase letters spanning the range from 1 to the uppercase letter of the index, e.g., $m=1,\ldots,M$. A third-order tensor $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is addressed by three indices $\{i_{n}\}$, $n=1,2,3$. Each $i_{n}$ usually addresses the $n$th mode of $\mathcal{A}$, while such convention may not be strictly followed when the context is clear. The $i_{3}$th mode-3 slice, a.k.a. the frontal slice, of $\mathcal{A}$ is denoted as $\mathbf{A}^{(i_{3})}$, a matrix obtained by fixing the mode-3 index $i_{3}$, i.e., $\mathbf{A}^{(i_{3})}=\mathcal{A}(:,:,i_{3})$. The $(i_{1},i_{2})$th tube of $\mathcal{A}$, denoted as $\mathring{\mathbf{a}}_{i_{1}i_{2}}$, is a mode-$3$ vector obtained by fixing the first two mode indices, i.e., $\mathcal{A}(i_{1},i_{2},:)$.

t-SVD and tubal rank. We first review the t-SVD framework following the definitions in (Kilmer et al., 2013). The t-product (tensor-tensor product) between tensor $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ and $\mathcal{B}\in\mathbb{R}^{I_{2}\times J_{4}\times I_{3}}$ is defined as $\mathcal{A}*\mathcal{B}=\mathcal{C}\in\mathbb{R}^{I_{1}\times J_{4}\times I_{3}}$. The $(i_{1},j_{4})$th tube $\mathring{\mathbf{c}}_{i_{1}j_{4}}$ of $\mathcal{C}$ is computed as

[TABLE]

where $*$ denotes the circular convolution (Rabiner & Gold, 1975) between two tubes (vectors) of the same size and the respective t-product between tensors. The tensor conjugate transpose of $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is denoted as $\mathcal{A}^{\top}\in\mathbb{R}^{I_{2}\times I_{1}\times I_{3}}$, obtained by conjugate transposing each of the frontal slice and then reversing the order of transposed frontal slices $\mathcal{A}(:,:,i_{3})$. A tensor $\mathcal{I}\in\mathbb{R}^{I\times I\times I_{3}}$ is an identity tensor if its first mode-3 slice $\mathcal{I}^{(1)}$ is an $I\times I$ identity matrix and all the rest mode-3 slices, i.e. $\mathcal{I}^{(i_{3})}$ for $i_{3}=2,...,I_{3}$, are zero matrices. An orthogonal tensor is a tensor $\mathcal{Q}\in\mathbb{R}^{I\times I\times I_{3}}$ that satisfies the following condition,

[TABLE]

where $\mathcal{I}\in\mathbb{R}^{I\times I\times I_{3}}$ is an identity tensor and $*$ is the t-product defined above. If $\mathcal{A}$’s all mode-3 slices $\mathbf{A}^{(i_{3})}$, $i_{3}=1,...,I_{3}$ are diagonal matrices, it is called an f-diagonal tensor. Based on these definitions, the t-SVD of $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is defined as

[TABLE]

where $\mathcal{U}\in\mathbb{R}^{I_{1}\times I_{1}\times I_{3}}$, $\mathcal{V}\in\mathbb{R}^{I_{2}\times I_{2}\times I_{3}}$ are orthogonal tensors, and $\mathcal{S}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is an f-diagonal tensor. This t-SVD definition leads to a new tensor rank, the tubal rank, which is defined as the number of nonzero singular tubes of $\mathcal{S}$, i.e. $\#\{i:\mathcal{S}(i_{2},i_{2},:)\neq 0\}$, assuming $I_{1}\geq I_{2}$. Figure 2 is an illustration of t-SVD.

t-SVD via Fourier transform. t-SVD can be computed via the discrete Fourier transform (DFT) for better efficiency. We denote the Fourier transformed tensor $\mathcal{A}$ as $\mathcal{A}_{\mathscr{F}}$, obtained via fast Fourier transform (FFT) along mode-3, i.e., $\mathcal{A}_{\mathscr{F}}={\tt fft}(\mathcal{A},[\hskip 2.84526pt],3)$. The connection between t-SVD and DFT is detailed in (Kilmer et al., 2013).

Tubal Tensor nuclear norm. TNN is a convex relaxation for the tubal rank, as an average tubal multi-rank within the unit tensor spectral norm ball (Lu et al., 2018a). It can be defined via DFT, similar to t-SVD computation above. The TNN of $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is defined as

[TABLE]

where $\|\cdot\|_{*}$ denotes the matrix nuclear norm. Please refer to (Lu et al., 2018a) for a detailed derivation and complete theoretical analysis, e.g., on the tightness of the relaxation.

Note that when $I_{3}=1$, the above definitions for tensors will be equivalent to the counterparts for matrices.

4 Sparse Multilinear Regression with Tubal Rank Regularization

Tubal rank-based TNN has shown to be superior to Tucker rank-based SNN in tensor completion/recovery (Zhang & Aeron, 2017), and tensor robust PCA (Lu et al., 2016), which are all unsupervised learning settings. To our knowledge, there is no study on TNN in a supervised learning setting yet. In this work, we explore the supervised learning with TNN and study whether TNN can improve supervised learning, e.g., multilinear regression. In the following, we propose the Sturm model and derive the Sturm algorithm under the ADMM framework.

4.1 The Sturm model

We incorporate TNN in the multilinear regression problem, which trains a model from $M$ pairs of feature tensors and their response labels $(\mathcal{X}_{m}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}},y_{m})$ with $m=1,...,M$ to relate them via a coefficient tensor $\mathcal{W}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$. This can be achieved by minimizing a loss function, typically with certain regularization:

[TABLE]

where $L(\cdot)$ is a loss function, $\Omega(\cdot)$ is a regularization function, $\lambda$ is a balancing hyperparameter, and $\langle\mathcal{X},\mathcal{W}\rangle$ denotes the inner product (a.k.a. the scalar product) of two tensors of the same size defined as

[TABLE]

The Remurs (Song & Lu, 2017) model uses a conventional least square loss function and assumes $\mathcal{W}$ to be both sparse and low rank. The sparsity of $\mathcal{W}$ is regularized by an $\ell_{1}$ norm and the rank by an SNN norm. However, the SNN requires unfolding $\mathcal{W}$ into matrices, susceptible to losing some higher-order structural information. Moreover, it has been pointed out in (Romera-Paredes & Pontil, 2013) that SNN is not a tight convex relaxation of its target rank.

This motivates us to propose a Sparse tubal-regularized multilinear regression (Sturm) model which replaces SNN in Remurs with TNN. This leads to the following objective function

[TABLE]

where $\tau$ and $\gamma$ are hyperparameters, and $\|\mathcal{W}\|_{1}$ is the $\ell_{1}$ norm of tensor $\mathcal{W}$, defined as

[TABLE]

which is equivalent to the $\ell_{1}$ norm of its vectorized representation $\mathbf{w}$. Here, the TNN regularization term $\|\mathcal{W}\|_{TNN}$ enforces low tubal rank in $\mathcal{W}$. The trade-off between $\tau$ and $\gamma$ as well as the degenerated versions follow the analysis for the Remurs (Song & Lu, 2017).

4.2 The Sturm algorithm via ADMM

ADMM (Boyd et al., 2011) is a standard solver for Problem (7). Thus, we derive an ADMM algorithm to optimize the Sturm objective function. We begin with introducing two auxiliary variables, $\mathcal{A}$ and $\mathcal{B}$ to disentangle the TNN and the $\ell_{1}$-norm regularization:

[TABLE]

Then, we introduce two Lagrangian dual variables $\mathcal{P}$ (for $\mathcal{A}$) and $\mathcal{Q}$ (for $\mathcal{B}$). With a Lagrangian constant $\rho$, the augmented Lagrangian becomes,

[TABLE]

We further introduce two scaled dual variables $\mathcal{P^{\prime}}=\frac{1}{\rho}\mathcal{P}$ and $\mathcal{Q^{\prime}}=\frac{1}{\rho}\mathcal{Q}$ only for notation convenience. Next, we derive the update from iteration $k$ to $k+1$ by taking an alternating strategy, i.e., minimizing one variable with all other variables fixed.

Updating $\mathcal{A}^{k+1}$:

[TABLE]

This can be rewritten as a linear-quadratic objective function by vectorizing all the tensors. Specifically, let $\mathbf{a}={\tt vec}(\mathcal{A})$, $\mathbf{w}^{k}={\tt vec}(\mathcal{W}^{k})$, $\mathbf{p^{\prime}}^{k}={\tt vec}(\mathcal{P^{\prime}}^{k})$, $\mathbf{y}=[y_{1}\cdots y_{M}]^{\top}$, $\mathbf{x}_{m}={\tt vec}(\mathcal{X}_{m})$, and $\mathbf{X}=[\mathbf{x}_{1}\cdots\mathbf{x}_{M}]^{\top}$. Then we get an equivalent objective function with the following solution:

[TABLE]

where $\mathbf{I}$ is an identity matrix. Note that this does not break/lose any structure because Eq. (11) and Eq. (12) are equivalent. $\mathcal{A}^{k+1}$ is obtained by folding (reshaping) $\mathbf{a}^{k+1}$ into a third-order tensor, denoted as $\mathcal{A}^{k+1}={\tt tensor}_{3}(\mathbf{a}^{k+1})$. Here, for a fixed $\rho$, we can avoid high per-iteration complexity of updating $\mathbf{a}^{k+1}$ by pre-computing a Cholesky decomposition of $(\mathbf{X}^{\top}\mathbf{X}+\rho\mathbf{I})$, which does not change over iterations.

Updating $\mathcal{B}^{k+1}$:

[TABLE]

This means that $\mathcal{B}^{k+1}$ can be solved by passing parameter $\frac{\tau}{\rho}$ to the proximal operator of the TNN (Zhang et al., 2014; Zhang & Aeron, 2017). The proximal operator for the TNN at tensor $\mathcal{T}$ with parameter $\mu$ is denoted by ${\tt prox}_{\mu\|\cdot\|_{TNN}}(\mathcal{T})$ and defined as

[TABLE]

where $\|\cdot\|_{F}$ is the Frobenius norm defined as $\|\mathcal{T}\|_{F}=\sqrt{\langle\mathcal{T},\mathcal{T}\rangle}$ using Eq. (6). The proximal operator for TNN can be more efficiently computed in the Fourier domain, as in Algorithm 1, where in Step 4, $(\mathbf{s}-\mu)_{+}=\max\{\mathbf{s}-\mu,0\}$ and ${\tt diag}(\mathbf{a})$ denotes a diagonal matrix whose diagonal elements are from $\mathbf{a}$.

Updating $\mathcal{W}^{k+1}$:

[TABLE]

It can be solved by calling the proximal operator of the $\ell_{1}$ norm with parameter $\frac{\gamma}{2\rho}$, which is simply the element-wise soft-thresholding, i.e.

[TABLE]

Updating $\mathcal{P}^{k+1}$ and $\mathcal{Q}^{k+1}$: The updates of $\mathcal{P}$ and $\mathcal{Q}$ are simply dual ascent steps:

[TABLE]

The complete procedure is summarized in Algorithm 2. The code will be made publicly available via GitHub.

4.3 Computational complexity

Finally, we analyze the per-iteration computational complexity of Algorithm 2. Let $I=I_{1}I_{2}I_{3}$. Step 3 takes $O(IM+\min\{M^{2},I^{2}\})$. Step  4 takes $O(\min\{I_{1},I_{2}\}I)$ for the singular value thresholding in the Fourier domain, plus $O(I\log(I_{3}))$ for ${\tt fft}$ and ${\tt ifft}$. Step 5 takes $O(I)$ because the proximal operation for $\ell_{1}$ norm is element-wise. Step 6 and Step 7 take $O(I)$. As a result, in a high dimensional (or small sample) setting where $I\gg M$, the per-iteration complexity is $O(I(\log(I_{3})+M))$.

5 Experiments

We evaluate our Sturm algorithm on four binary classification problems with six datasets from three public fMRI repositories. While known works typically focus on either resting-state or task-based fMRI only and rarely both, here we study both types. We test the efficacy of the Sturm approach against five state-of-the-art algorithms and three additional variations, in terms of both classification accuracy and sparsity.

5.1 Classification problems and datasets

Resting-state fMRI for disease diagnosis. Resting-state fMRI is scanned when the subject is not doing anything, i.e., at rest. It is commonly used for brain disease diagnosis, i.e., classifying clinical population. In this paper, we consider only the binary classification of patients (positive) and health control subjects (negative).

Task-based fMRI for neural decoding. Task-based fMRI is scanned when the subject is performing certain tasks, such as viewing pictures or reading sentences (Wang et al., 2003). It is commonly used for studies decoding brain cognitive states, or neural decoding. The objective is to classify (decode) the tasks performed by subjects using the fMRI information. In this paper, we consider only binary classification of two different tasks.

Chosen datasets. We study four fMRI classification problems on six datasets from three public repositories, with the key information summarized in Table 2. Two are disease diagnosis problems on resting-state fMRI, and the other two are neural decoding problems on task-based fMRI, as described below.

• •

Resting 1 – ABIDENYU&UM: The Autism Brain Imaging Data Exchange (ABIDE)111http://fcon_1000.projects.nitrc.org/indi/abide (Craddock et al., 2013) consists of patients with autism spectrum disorder (ASD) and healthy control subjects. We chose the largest two subsets contributed by New York University (NYU) and University of Michigan (UM). The fMRI data has been preprocessed by the pipeline of Configurable Pipeline for the Analysis of Connectomes (CPAC). Quality control was performed by selecting the functional images with quality ‘OK’ reported in the phenotype data.

• •

Resting 2 – ADHD-200NYU: We chose the NYU subset from the Attention Deficit Hyperactivity Disorder (ADHD) 200 (ADHD-200) dataset222http://neurobureau.projects.nitrc.org/ADHD200/Data.html (Bellec et al., 2017), with ADHD patents and healthy controls. The raw data is preprocessed by the pipeline of Neuroimaging Analysis Kit (NIAK).

• •

**Task 1 – Balloon vs Mixed gamble **: We chose two gamble-related datasets from the OpenfMRI repository333Data used in this paper are available at OpenfMRI: https://legacy.openfmri.org, now known as OpenNeuro: https://openneuro.org. (Poldrack et al., 2013) project to form a classification problem. They are 1) Balloon analog risk-taking task (BART) and 2) Mixed gambles task.

• •

Task 2 – Simon vs Flanker: We chose another two recognition and response related tasks from OpenfMRI for binary classification. They are 1) Simon task and 2) Flanker task.

Resting-state fMRI preprocessing. The raw resting-state brain fMRI data is 4-D. We follow typical approaches to reduce the 4-D data to 3-D by either taking the average (He et al., 2017) or the amplitude (Yu-Feng et al., 2007) of low frequency fluctuation of voxel values along the time dimension. We perform experiments on both and report the best results.

Task-based fMRI preprocessing. Following (Poldrack et al., 2013), we re-implemented a preprocessing pipeline to process the OpenfMRI data to obtain the 3D statistical parametric maps (SPMs) for each brain condition with a standard template. We used the same criteria as in (Poldrack et al., 2013) to selected one contrast (one specific brain condition over experimental conditions) per task for classification.

The tubal mode of fMRI. Figure 1 illustrates how fMRI scan is obtained along the axial direction, which is mode 3 in tensor representation. Each image along the diagonal is a mode-3 (frontal) slice. Therefore, it is a natural choice to consider mode 3 as the tubal mode to apply Sturm.

5.2 Algorithms

We evaluate Sturm and Sturm + SVM (support vector machine) against the following five algorithms and three additional algorithms via combination with SVM.

• •

SVM: We chose linear SVM for both speed and prediction accuracy. (We studied both the linear and Gaussain RBF kernel SVM and found the linear one performs better on the whole.)

• •

Lasso: It is a linear regression method with the $\ell_{1}$ norm regularization.

• •

Elastic Net (ENet): It is a linear regression method with $\ell_{1}$ and $\ell_{2}$ regularization.

• •

Remurs (Song & Lu, 2017): It is a multilinear regression model with $\ell_{1}$ norm and Tucker rank-based SNN regularization.

• •

Multi-way Multi-level Kernel Modeling (MMK) (He et al., 2017): It is a kernelized CP tensor factorization method to learn nonlinear features from tensors. Gaussain RBF kernel MMK is used with pre-computed kernel SVM.

SVM, Lasso, and ENet take vectorized fMRI data as input while Remurs and MMK directly take 3-D fMRI tensors as input. Lasso, ENet, Remurs, and Sturm can also be used for (embedded) feature selection. Therefore, we can add an SVM after each of them to obtain Lasso + SVM, ENet + SVM, Remurs + SVM and Sturm + SVM. The code for Sturm is built on the software library from (Lu, 2016; Lu et al., 2018b). Remurs, Lasso, and ENet are implemented with the SLEP package (Liu et al., 2009). MMK code is kindly provided by the first author of (He et al., 2017).

5.3 Algorithm and evaluation settings

Model hyperparameter tuning. Default settings are used for all existing algorithms. For Sturm, we follow the Remurs default setting (Song & Lu, 2017) to set $\rho$ to 1 and use the same set $\{10^{-3},5\times 10^{-3},10^{-2},\dots,5\times 10^{2},10^{3}\}$ for $\tau$ and $\gamma$, while scaling the first term in Eq. (7) by a factor $\alpha=\sqrt{(\max(I_{1},I_{2})\times I_{3})}$ to better balance the scales of the loss function and regularization terms (Lu et al., 2016, 2018a).

Image resizing. To improve computational efficiency and reduce the small sample size problem (and overfitting), the input 3-D tensors are further re-sized into three different sizes with a factor $\beta$, choosing from $\{0.3,0.5,0.7\}$.

Feature selection. In Lasso + SVM, ENet + SVM, Remurs + SVM, and Sturm + SVM, we rank the selected features by their associated absolute values of $\mathcal{W}$ in the descending order and feed the top $\eta\%$ of the features to SVM. We study five values of $\eta$: $\{1,5,10,50,100\}$.

Evaluation metric and method. The classification accuracy is our primary evaluation metric, and we also examine the sparsity of the obtained solution for all algorithms except SVM and MMK. For a particular binary classification problem, we perform ten-fold cross validation and report the mean and standard deviation of the classification accuracy and sparsity over ten runs. For each of the ten (test) folds, we perform an inner nine-fold cross validation using the remaining nine folds to determine $\tau$ and $\gamma$ (jointly for ENet, Remurs, and Sturm on a $13\times 13$ grid), $\beta$, and $\eta$ above that give the highest classification accuracy, with the corresponding sparsity recorded. The sparsity is calculated as the ratio of the number of zeros in the output coefficient tensor $\mathcal{W}$ to its size $I_{1}\times I_{2}\times I_{3}$. In general, higher sparsity implies better interpretability (Hastie et al., 2015).

5.4 Results and discussion

Table 3 reports the classification accuracy for all algorithms except MMK. This is because all MMK results are below 60% on all four problems, possibly due to the default settings of the CP rank and SVM kernel. (We have tried a few alternative settings without improvement, though further tuning may still lead to better results.) Table 4 presents the respective sparsity values except SVM, which uses all features so the sparsity is zero. In both tables, the best results are highlighted in bold, with the second best ones underlined.

Performance on Resting 1 & 2. For these two resting-state problems, Sturm + SVM has the highest accuracy of 65.45%, and Lasso is the second best with 2.45% lower accuracy. In terms of sparsity, Sturm + SVM and Sturm are the top two. Specifically, on Resting 1, Remurs + SVM and Sturm + SVM are the top two algorithms with almost identical accuracy of 64.67% and 64.66%, respectively. Moreover, Sturm + SVM also has the highest sparsity of 0.87 and Sturm has the second-highest sparsity of 0.86. For Resting 2, Sturm + SVM has outperformed all other algorithms on both accuracy (66.24%) and sparsity (0.99).

Performance on Task 1 & 2. For these two task-based problems, Sturm has outperformed all other algorithms in accuracy, with 89.10% on Task 1 and 86.89% on Task 2. Lasso is again the second best in accuracy. ENet and ENet + SVM has the best sparsity of 0.96 while their accuracy values are only 81.87% and 74.21%, respectively. Sturm + SVM has significant drop in accuracy compared with Sturm alone, and Lasso + SVM, ENet + SVM and Remurs + SVM all have lower accuracy compared to without SVM.

Summary. There are four key observations on the whole:

• •

Sturm has the best overall accuracy of 75.38%. Sturm has outperformed Remurs in accuracy for all four classification problems. The only difference between Sturm and Remurs is replacing SNN with TNN. Therefore, this superiority indicates that tubal rank-based TNN is superior to Tucker rank-based SNN in the supervised, regression setting.

• •

The reported sparsity corresponds to the best solution determined via nine-fold cross validation. Lasso and Lasso + SVM have the lowest, i.e., poorest, sparsity. ENet and ENet + SVM also have much lower sparsity than Remurs/Sturm and their + SVM versions, more than 0.10 (10%) lower. On the other hand, ENet and ENet + SVM have the highest sparsity on task-based fMRI while getting a solution closely to zero sparsity on Resting 2, showing high variation.

• •

Performing SVM after the four regression methods can improve the classification accuracy (though not always) on resting-state fMRI, while it has degraded their classification performance on task-based fMRI in all cases. In particularly, Lasso + SVM and Sturm + SVM on Task 1, and Lasso + SVM and ENet + SVM on Task 2 have dropped more than 10% in accuracy.

• •

Disease diagnosis on resting-state fMRI is significantly more challenging than neural decoding on task-based fMRI. Although it should be aware from Table 2 that the number of samples is different for the resting-state and task-based fMRI, engaged brain activities are generally easier to classify than those not engaged and the difference is consistent with results reported in the literature.

5.5 Analysis

Hyperparameter sensitivity. Figure 3 illustrates the classification performance sensitivity of Sturm on its two hyperparameters $\tau$ and $\gamma$ for two problems: Resting 2 and Task 2. In general, the lower right, i.e., large $\tau$ and small $\gamma$ values, shows higher accuracy. This implies that the tubal rank regularization helps more in improving the accuracy than sparsity regularization. However, Fig. 3(a) shows poorer smoothness than Fig. 3(b), another indication of resting-state fMRI being more challenging than task-based fMRI. This makes hyperparameter tuning more difficult for resting-state fMRI, (partly) causing the poorer classification performance than task-based fMRI.

Convergence analysis. Figure 4 shows the convergence of $\mathcal{W}$ and the Sturm objective function value in (7) on Resting 2 and Task 2. It can be seen that $\mathcal{W}$ has a fast convergence speed in both cases, though the objective function converges at a slower rate.

6 Conclusion

In this paper, we proposed a sparse tubal rank regularized multilinear regression model (Sturm). It performs regression with penalties on the tubal tensor nuclear norm, a convex relaxation of the tubal rank, and the standard $\ell_{1}$ norm. To the best of our knowledge, this is the first supervised learning method using TNN. The Sturm algorithm is derived via a standard ADMM framework.

We evaluated Sturm (and Sturm + SVM) in terms of classification accuracy and sparsity against eight other methods on four binary classification problems from three public fMRI repositories. The datasets include both resting-state and task-based fMRI, unlike most existing works focusing on only one of them. The results showed the superior overall performance of Sturm (and Sturm + SVM in some cases) over other methods and confirmed the benefits of TNN and tubal rank regularization in a supervised setting.