Title: SUBSPACE NETWORK: DEEP MULTI-TASK CENSORED REGRESSION FOR MODELING NEURODEGENERATIVE DISEASES

Abstract: Over the past decade a wide spectrum of machine learning models have been developed to model the neurodegenerative diseases, associating biomarkers, especially non-intrusive neuroimaging markers, with key clinical scores measuring the cognitive status of patients. Multi-task learning (MTL) has been extensively explored in these studies to address challenges associated to high dimensionality and small cohort size. However, most existing MTL approaches are based on linear models and suffer from two major limitations: 1) they cannot explicitly consider upper/lower bounds in these clinical scores; 2) they lack the capability to capture complicated non-linear effects among the variables. In this paper, we propose the Subspace Network, an efficient deep modeling approach for non-linear multi-task censored regression. Each layer of the subspace network performs a multi-task censored regression to improve upon the predictions from the last layer via sketching a low-dimensional subspace to perform knowledge transfer among learning tasks. Under mild assumptions, for each layer the parametric subspace can be recovered using only one pass of training data. Empirical results demonstrate that the proposed subspace network quickly picks up correct parameter subspaces, and outperforms state-of-the-arts in predicting neurodegenerative clinical scores using information in brain imaging.

1. INTRODUCTION:

Recent years have witnessed increasing interests on applying machine learning (ML) techniques to analyze biomedical data. Such data-driven approaches deliver promising performance improvements in many challenging predictive problems. For example, in the field of neurodegenerative diseases such as Alzheimer's disease and Parkinson's disease, researchers have exploited ML algorithms to predict the cognitive functionality of the patients from the brain imaging scans, e.g., using the magnetic resonance imaging (MRI) as in Adeli-Mosabbeb et al. (2015) ; Zhang et al. (2012) ; Zhou et al. (2011b) . As a key finding, there are typically various types of prediction targets (e.g., cognitive scores), and they can be jointly learned using multi-task learning (MTL), e.g., Caruana (1998) ; Evgeniou & Pontil (2004) ; Zhang et al. (2012) , where the predictive information is shared and transferred among related models to reinforce their generalization performance.

Two challenges persist despite the progress of applying MTL in disease modeling problems. First, it is important to notice that clinical targets, different from typical regression targets, are often naturally bounded. For example, the result from Mini-Mental State Examination (MMSE) test, a key reference for deciding cognitive impairments, ranges from 0 to 30 (a healthy subject): a smaller score indicates a higher level of cognitive dysfunction (please refer to Tombaugh & McIntyre (1992) ). Other cognitive scores, such as Clinical Dementia Rating Scale (CDR) Hughes et al. (1982) and Alzheimer's Disease Assessment Scale-Cog (ADAS-Cog) Rosen et al. (1984) , also have specific upper and lower bounds. Most existing approaches, e.g., Zhang et al. (2012) ; Zhou et al. (2011b) ; Poulin et al. (2011) , relied on linear regression without considering the range constraint, partially due to the fact that mainstream MTL models for regression, e.g., Jalali et al. (2010) ; Argyriou et al. (2007) ; Zhang et al. (2012) ; Zhou et al. (2011a) , are developed using the least squares loss and cannot be directly extended to censored regressions. As the second challenge, a majority of MTL research focused on linear models because of computational efficiency and theoretical guarantees. However, linear models cannot capture the complicated non-linear relationship between features and clinical A critical difference on subspace decoupling between previous studies such as Mardani et al. (2015 ) Shen et al. (2016) and our method lies in our assumption of a low-rank structure on the parameter space among tasks rather than the original feature space. • By leveraging the recent advances in online subspace sensing Mardani et al. (2015) ; Shen et al. (2016) , we show that the parametric subspace can be recovered for each layer with feeding only one pass of the training data, allowing for more efficient layer-wise training. Synthetic experiments verify the technical claims of the proposed SN, and it outperforms various state-of-the-arts methods in modeling neurodegenerative diseases on real datasets.

2. MULTI-TASK CENSORED REGRESSION VIA PARAMETER SUBSPACE:



SKETCHING AND REFINEMENT:

In censored regression, we are given a set of

N observations D = {(x i , y i )} N i=1 of D dimensional fea- ture vectors {x i ∈ R D } and T corresponding outcomes {y i ∈ R T + }, where each outcome y i,t ∈ R + , t ∈ {1, • • • , T },

can be cognitive scores (e.g., MMSE and ADAS-Cog) or other biomarkers of interest such as proteomics1 . For each outcome, the censored regression assumes a nonlinear relationship between the features and the outcome through a rectified linear unit (ReLU) transformation, i.e., y i,t = ReLU W t x i + where W t ∈ R D is the coefficient for input features, is i.i.d. noise, and ReLU is defined by ReLU(z) = max(z, 0). We can thus collectively represent the censored regression for multiple tasks by:

EQUATION

where W = [W 1 , . . . , W T ] ∈ R T ×D is the coefficient matrix. We consider the regression problem for each outcome as a learning task, and one commonly used task relationship assumption is that the transformation matrix W ∈ R T ×D belongs to a linear low-rank subspace U. The subspace allows us to represent W as product of two matrices, W = U V , where columns of U ∈ R T ×R = [U 1 , . . . , U T ] span the linear subspace U, and V ∈ R R×D is the embedding coefficient. We note that the output y can be entry-wise decoupled, such that for each component y i,t = ReLU(U t V x i + ).

By assuming Gaussian noise ∼ N (0, σ 2 ), we derive the following likelihood function:

Pr(y i,t , x i |U t , V ) = φ y i,t -U t V x i σ I(y i,t ∈ (0, ∞)) + 1 -Q 0 -U t V x i σ I(y i,t = 0),

Algorithm 1 Single-layer parameter subspace sketching and refinement.

Require: Training data D = {(x i , y i )} N i=1 , rank parameters λ and R, Ensure: parameter subspace U , parameter sketch V Initialize U -at random for i = 1, . . . , N do // 1. Sketching parameters in the current subspace

V + = arg min V -log Pr(y i , x i |U -, V ) + λ 2 V 2 F // 2. Parallel subspace refinement {U t } T t=1 for t = 1, . . . , T do U + t = arg min Ut -log Pr(y i,t , x i |U t , V + ) + λ 2 U t 2 2 end for Set U -= U + , V -= V + end for

where φ is the probabilistic density function of the standardized Gaussian N (0, 1) and Q is the standard Gaussian tail. σ controls how accurately the low-rank subspace assumption can fit the data. Note that other noise models can be assumed here too. The likelihood of (x i , y i ) pair is thus given by:

Pr(y i , x i |U, V ) = T t=1 φ y i,t -U t V x i σ I(y i,t ∈ (0, ∞)) + 1 -Q - U t V x i σ I(y i,t = 0) .

The likelihood function allows us to estimate subspace U and coefficient V from data D. To enforce a low-rank subspace, one common approach is to impose a trace norm on U V , where trace norm of a matrix A is defined by A * = j s j (A) and s j (A) is the jth singular value of A.

Since U V * = min U,V 1 2 ( U 2 F + V 2 F

), e.g., see Srebro et al. (2005) ; Mardani et al. (2015) , the objective function of multi-task censored regression problem is given by:

min U,V - N i=1 log Pr(y i , x i |U, V ) + λ 2 ( U 2 F + V 2 F ).

(2)

2.1. AN ONLINE ALGORITHM:

We propose to solve the objective in (2) via the block coordinate descent approach which is reduced to iteratively updating the following two subproblems:

V + = arg min V - N i=1 log Pr(y i , x i |U -, V ) + λ 2 V 2 F , (P:V) U + = arg min U - N i=1 log Pr(y i , x i |U, V + ) + λ 2 U 2 F . (P:U)

Define the instantaneous cost of the i-th datum:

g(x i , y i , U, V ) = -log Pr(x i , y i |U, V ) + λ 2 U 2 F + λ 2 V 2

F , and the online optimization form of (2) can be recast as an empirical cost minimization given below:

min U,V 1 N N i=1 g(x i , y i , U, V ).

According to the analysis in Section 2.2, one pass of the training data can warrant the subspace learning problem. We outline the solver for each subproblem as follows:

Problem (P:V) sketches parameters in the current space. We solve (P:V) using gradient descent. The parameter sketching couples all the subspace dimensions in V (not decoupled as in Shen et al. (2016) ), and thus we need to solve this collectively. The update of V (V + ) can be obtained by solving the online problem given below:

min V g(x i , y i ; U -, V ) ≡ - T t=1 log Pr(y i,t , x|U - t , V ) + λ 2 V 2 F = - T t=1 log φ y i,t -U - t V x σ I(y i,t ∈ (0, ∞)) + 1 -Q -U - t V x σ I(y i,t = 0) + λ 2 V 2 F .

V + can be computed by the following gradient update:

V + = V --η∇ V g(x i , y i ; U -, V + )

, where the gradient is given by:

∇ V g(x i , y i ; U -, V + ) = λV + T t=1    - yi,t-(U - t ) V xi σ 2 U - t x i y i,t ∈ (0, ∞) φ(zt) σ[1-Q(zi,t)] U - t x T i y i,t = 0 where z i,t = σ -1 (-U - t V x

). The Alg. 3 for solving (P:V) is summarized in the Appendix.

Problem (P:U) refines the subspace U + based on sketching. We solve (P:U) using stochastic gradient descent (SGD). We note that the problem is decoupled for different subspace dimensions t = 1, . . . , T (i.e., rows of U ). With careful parallel design, this procedure can be done very efficiently. Given a training data point (x i , y i ), the problem related to the t-th subspace basis is:

min Ut g t (x i , y i,t ; U t , V + ) ≡ -log Pr(y i,t , x i |U t , V + ) + λ 2 U t 2 2 = -log φ y i,t -U t V + x i σ I(y i,t ∈ (0, ∞)) + 1 -Q -U t V + x i σ I(y i,t = 0) + λ 2 U t 2 2 .

We can revise subspace by the following gradient update:

U + t = U - t -µ t ∇ Ut g t (x i , y i,t ; U t , V + )

, where the gradient is given by:

∇ Ut g t (x i , y i,t ; U i,t , V + ) = λU t + - yi,t-U t V + x σ 2 V + x i y i,t ∈ (0, ∞) φ(zi,t) σ[1-Q(zi,t)] V + x i y i,t = 0 where z i,t = σ -1 (-U t V + x i ).

We summarize the procedure in Algorithm 1 and show in Section 2.2 that under mild assumptions this procedure will be able to capture the underlying subspace structure in the parameter space with just one pass of the data.

2.2. THEORETICAL RESULTS:

We establish both asymptotic and non-asymptotic convergence properties for Algorithm 1. The proof scheme is inspired by a series of previous works: Mardani et al. (2013; 2015) ; Shen et al. (2016) . We briefly present the proof sketch, and more proof details can be found in Appendix. At each iteration i = 1, 2, ..., N , we sample (x i , y i ), and let U i , V i denote the intermediate U and V , to be differentiated from U t , V t which are the t-th columns of U, V . For the proof feasibility, we assume that {(x i , y i )} N i=1 are sampled i.i.d., and the subspace sequence {U i } N i=1 lies in a compact set. Asymptotic Case: To estimate U , the Stochastic Gradient Descent (SGD) iterations can be seen as minimizing the approximate cost 1 N N i=1 g (x i , y i , U, V ), where g is a tight quadratic surrogate for g based on the second-order Taylor approximation around U N -1 . Furthermore, g can be shown to be smooth, by bounding its first-order and second-order gradients w.r.t. each U t (similar to Appendix 1 of Shen et al. (2016) ).

Following Mairal et al. (2010) ; Mardani et al. (2015) , it can then be established that, as N → ∞, the subspace sequence {U i } N i=1 asymptotically converges to a stationary-point of the batch estimator, under a few mild conditions. We can sequentially show: 1)

N i=1 g (x i , y i , U i , V i ) asymptotically converges to N i=1 g(x i , y i , U i , V i ),

according to the quasi-martingale property in the almost sure sense, owing to the tightness of g ; 2) the first point implies convergence of the associated gradient sequence, due to the regularity of g; 3) g t (x i , y i , U, V ) is bi-convex for block variables U t and V .

Non-Asymptotic Case: When N is finite, Mardani et al. (2013) asserts that the distance between successive subspace estimates will vanish as fast as o(1/i): U i -U i-1 F ≤ B i , for some constant B that is independent of i and N . Following Shen et al. (2016) to exploit the unsupervised formulation of regret analysis as in Kasiviswanathan et al. (2012) ; Shalev-Shwartz et al. (2012) , we can similarly obtain a tight regret bound that will again vanish if N → ∞.

3. SUBSPACE NETWORK VIA HIERARCHICAL SKETCHING AND REFINEMENT:

The single layer model in (1) has very limited power to capture the highly nonlinear regression relationships, as the parameters are linearly linked to the subspace except for a ReLU operation.

Algorithm 2 Network expansion via hierarchical parameter subspace sketching and refinement.

Require::

Training data D = {(x i , y i )}, target network depth K. Ensure: The deep subspace network f

Set f [0] (x) = y and solve f [0] using Algorithm 1. for k = 1, . . . , K -1 do // 1. Subspace sketching based on the current subspace using Algorithm 1:

U * [k] , V * [k] = arg min U [k] ,V [k] E (x,y)∼D (y, ReLU U [k] V [k] f [k-1] (x) ) ,

// 2. Expand the layer using the refined subspace as our new network:

f [k] (x) = ReLU U * [k] V * [k] f [k-1] (x) end for return f = f [K]

However, the single-layer procedure in Algorithm 1 has provided a building block, based on which we can develop an efficient algorithm to train a deep subspace network (SN) in a greedy fashion. We thus propose a network expansion procedure to overcome such limitation.

After we obtain the parameter subspace U and sketch V for the single-layer case (1), we project the data points by x = ReLU(U V x). A straightforward idea of the expansion is to use (x, y) as the new samples to train another layer. Let f [k-1] denote the network structure we obtained before the k-th expansion starts, k = 1, 2, ..., K -1, the expansion can recursively stack more ReLU layers:

EQUATION

However, we observe that simply stacking layers by repeating (3) many times can cause substantial information loss and degrade the generalization performance, especially since our training is layerby-layer without "looking back" (i.e., top-down joint tuning). Inspired by deep residual networks by He et al. (2016) that exploit "skip connections" to pass lower-level data and features to higher levels, we concatenate the original samples with the newly transformed, censored outputs after each time of expansion, i.e., reformulating x = [ReLU(U V x); x] (similar manners could be found in Zhou & Feng (2017) ).The new form after the expansion is given below:

f [k] (x) = ReLU U [k] V [k] [f [k-1] (x); x] + .

We summarize the network expansion process in Algorithm 2. The architecture of the resulting SN is illustrated in Figure 1 . Compared to the single layer model (1), SN gradually refines the parameter subspaces by multiple stacked nonlinear projections. It is expected to achieve superior performance due to the higher learning capacity, and can also be viewed as a gradient boosting method. Meanwhile, the layer-wise low-rank subspace structural prior would improve generalization further compared to the naive multi-layer networks.

4. EXPERIMENT:



4.1. SIMULATIONS ON SYNTHETIC DATA:

Subspace recovery in a single layer model. We first evaluate the subspace recovered by the proposed Algorithm 1 using synthetic data. We generated X ∈ R N ×D , U ∈ R T ×R and V ∈ R R×D , all as i.i.d. random Gaussian matrices. The target matrix Y ∈ R N ×T was then synthesized using (1). We set N = 5, 000, D = 200, T = 100 R = 10, and use random noise ∼ N (0, 3 2 ).

Figure 2a shows the plot of subspace difference between the ground-truth U and the learned subspace U i throughout the iterations, i.e., U -U i F / U F w.r.t. i. This result verifies that Algorithm 1 is able to correctly find and smoothly converge to the underlying low-rank subspace of the synthetic data.

The objective values throughout the online training process of Algorithm 1 are plotted in Figure 2b . We further show the plot of iteration-wise subspace differences, defined as U i -U i-1 F / U F , in Figure 2c , which complies with the o(1/t) result in our non-asymptotic analysis. Moreover, the distribution of correlation between recovered weights and true weights for all tasks is given in Figure 3 in Appendix, with most predicted weights having correlations with ground truth of above 0.9.

Subspace recovery in a multi-layer subspace network. We re-generated synthetic data by repeatedly applying (1) for three times, each time following the same setting as the single-layer model. A three-layer SN was then learned using Algorithm 2. As one simple baseline, a multi-layer perceptron (MLP) is trained, whose three hidden layers have the same dimensions as the three ReLU layers of the SN. Inspired by Xue et al. (2013) ; Sainath et al. (2013) ; Wang et al. (2015) , we then applied low-rank matrix factorization to each layer of MLP, with the same desired rank R, creating the factorized MLP (f-MLP) baseline that has the identical architecture (including both ReLU hidden layers and linear bottleneck layers) to SN. We further re-trained the f-MLP on the same data from end to end, leading to another baseline, the retrained factorized MLP (rf-MLP).

Table 1 evaluates the subspace recovery fidelity in three layers, using three different metrics: (1) the maximum mutual coherence of all column pairs from two matrices, defined in (Candes & Romberg (2007) ) as a classical measurement on how correlated the two matrices' column subspaces are;

(2) the mean mutual coherence of all column pairs from two matrices; (3) the subspace difference defined the same as in the single-layer case2 . Note that the two mutual coherence-based metrics are immune to linear transformations of subspace coordinates, to which the 2 -based subspace difference might become fragile. SN achieves clear overall advantages under all three measurements, over f-MLP and rf-MLP. More notably, while the performance margin of SN in subspace difference seems to be small, the much sharper margins, in two (more robust) mutual coherence-based measurements, suggest that the recovered subspaces by SN are significantly better aligned with the groundtruth.

Benefits of Going Deep. We re-generate synthetic data again in the same way as the first single-layer experiment; yet differently, we now aim to show that a deep SN will boost performance over single-layer subspace recovery, even the data generation does not follow a known multi-layer model. We compare SN (both 1-layer and 3-layer) with two carefully chosen sets of state-of-art approaches: (1) single and multi-task "shallow" models; (2) deep models. For the first set, the least squares (LS) is treated as a naive baseline, while ridge (LS + 2 ) and lasso (LS + 1 ) regressions are considered for shrinkage or variables selection purpose; Censor regression, also known as the Tobit model, is a non-linear method to predict bounded targets , e.g., Berberidis et al. (2016) . Multi-task models with regularizations on trace norm (Multi Trace) and 2,1 norm (Multi 2,1 ) have been demonstrated to be successful on simultaneous structured/sparse learning, e.g., Yang et al. (2010) ; Zhang et al. (2013) . 3 We also verify the benefits of accounting for boundedness of targets (Uncensored vs. Censored) in both single-task and multi-task settings, with best performance reported for each scenario (LS + 1 for single-task and Multi Trace for multi-task). For the set of deep model baselines, we construct three DNNs for fair comparison: i) A 3-layer fully connected DNN with the same architecture as SN, with a plain MSE loss; ii) A 3-layer fully connected DNN as i) with ReLU added for output layer before feeding into the MSE loss, which naturally implements non-negativity censored training and evaluation; iii) A factorized and re-trained DNN from ii), following the same procedure of rf-MLP in the multi-layer synthetic experiment. Apparently, ii) and iii) are constructed to verify if DNN also benefits from the censored target and the low-rank assumption, respectively.

We performed 10-fold random-sampling validation on the same dataset, i.e., randomly splitting into training and validation data 10 times. For each split, we fit model on training data and evaluated performance on validation data. Average normalized mean square error (ANMSE) across all tasks is obtained as the overall performance for each split. For methods without hyper parameters (least square and censor regression), an average of ANMSE for 10 splits is regarded as the final performance; for methods with tunable parameters, e.g., λ in lasso, we perform a grid search on λ values and choose the optimal ANMSE result. We consider different splitting size with training samples containing [40%, 50%, 60%, 70%, 80%] of all the samples.

Table 2 further compares the performance of all approaches. We see that: (1) all censored models significantly outperform their uncensored counterparts, verifying the necessity of adding censoring targets for regression. Therefore, we will use censored baselines hereinafter, unless otherwise specified;

(2) the more structured MTL models tend to outperform single task models by capturing task relatedness. That is also evidenced by the performance margin of DNN iii over DNN i; (3) the nonlinear models are undoubtedly more favorable: we even see the single-task Tobit model to outperform MTL models; (4) As a nonlinear, censored MTL model, SN combines the best of them all, accounting for its superior performance over all competitors. In particular, even a 1-layer SN already produces comparable performance to the 3-layer DNN iii (which also a nonlinear, censored MTL model trained with back-propagation, with three times the parameter amount of SN), thanks to SN's theoretically solid online algorithm in characterizing subspaces.

Furthermore, we grow the number of layers in SN from 2 to 20, finding that as a feed-forward model without end-to-end optimization, SN can also benefit from its increasing depth. As Table 3 reveals, SN steadily improves with more layers, until reaching a plateau at around 5 layers (as the underlying data distribution is relatively simple here). The observation is consistent among different splits.

4.2. EXPERIMENTS ON REAL DATA:

We evaluated SN in a real clinical setting with the goal to build models to predict important clinical scores representing a subject's cognitive status and signaling the progression of Alzheimer's disease (AD), from structural Magnetic Resonance Imaging (sMRI) data. AD is one major neurodegenerative disease that accounts for 60 to 80 percent of dementia. The National Institutes of Health has thus Calibration. In MTL formulations we typically assume that noise variance σ 2 is the same across all tasks, which may not be true in many cases. To deal with heterogeneous σ 2 among tasks, we design a calibration step in our optimization process, where we estimate task-specific σ2 t using y -ŷ 2 2 /N before ReLU, as the input for next layer and repeat on layer-wise. We compare performance of both non-calibrated and calibrated methods.

Performance. We adopt the same two sets of comparison baselines for the real data experiments, as the last synthetic experiment. Different from synthetic data where the low-rank structure is predefined, for real data, there is no groundtruth rank available and we have to try different rank assumptions. Table 4 compares the performances between σ 2 non-calibrated versus calibrated models. We observe a clear improvement by assuming different σ 2 across tasks. Table 5 shows the results for all comparison methods, with SN outperforming all else. Table 7 in Appendix shows the SN performance growth with increasing the number of layers. Table 8 in Appendix further reveals the performance of DNNs and SN using varying rank estimations in real data. As expected, the U-shape curve suggests that an overly low rank may not be informative enough to recover the original weight space, while a high-rank structure cannot enforce as strong a structural prior. However, the overall robustness of SN to rank assumptions is fairly remarkable: its performance under all ranks is competitive, consistently outperforming DNNs under the same rank assumptions and other baselines.

5. CONCLUSIONS AND FUTURE WORK:

In this paper we proposed a Subspace Network (SN), an efficient deep modeling approach for non-linear multi-task censored regression, where each layer of the subspace network performs a multi-task censored regression to improve upon the predictions from the last layer via sketching a low-dimensional subspace to perform knowledge transfer among learning tasks. We show that under mild assumptions, for each layer we can recover the parametric subspace using only one pass of training data. We demonstrate empirically that the subspace network can quickly capture correct parameter subspaces, and outperforms state-of-the-arts in predicting neurodegenerative clinical scores from brain imaging. Based on similar formulations, the proposed method can be easily extended to cases where the targets have nonzero bounds, or both lower and upper bounds. iteration t (t ≤ N ), we use the previous U t-1 to span the partial data at i = 1, 2, ..., t. Prompted by the alternating nature of iterations, we adopt a variant of the unsupervised regret to assess the goodness of online subspace estimates in representing the partially available data. With g t (x t , y t , U t-1 , V ) being the loss incurred by the estimate U t-1 for predicting the t-th datum, the cumulative online loss for a stream of size T is given by:

EQUATION

Further, we will assess the cost of the last estimate U T using:

EQUATION

We define C T as the batch estimator cost. For the sequence {U t } T t=1 , we define the online regret:

EQUATION

We investigate the convergence rate of the sequence {R T } to zero as T grows. 

i,t = σ -1 (-U - t V x i ). Compute the gradient for y t : ∇g t (x i , y i,t ; U -, V + ) =    - yi,t-(U - t ) V xi σ 2 U - t x i y i,t ∈ (0, ∞) φ(zi,t) σ[1-Q(zi,t)] U - t x i y i,t = 0 end for // 2. Update the current sketch V - V + = V --η T t=1 ∇g t (x, y t ; U -, V + ) + λV Set V -= V +

:

Without loss of generality, in this paper we assume that outcomes are lower censored at 0. By using variants of Tobit models, e.g., as inShen et al. (2016), the proposed algorithms and analysis can be extended to other censored models with minor changes in the loss function.

The higher in terms of the two mutual coherence-based metrics, the better subspace recovery is achieved.That is different from the subspace difference case where the smaller the better,

Least squares, ridge, lasso, and censor regression are implemented by MATLAB optimization toolbox. MTLs are implemented throughMALSAR Zhou et al. (2011a) with parameters carefully tuned.

Figures and Tables:

Figure 1: The proposed subspace network via hierarchical subspace sketching and refinement. targets. For example, Association et al. (2013) showed the early onset of Alzheimer's disease to be related to single-gene mutations on chromosomes 21, 14, and 1, and the effects of such mutations on the cognitive impairment are hardly linear (please refer to Martins et al. (2005); Sweet et al. (2012)). Recent advances in multi-task deep neural networks Seltzer & Droppo (2013); Zhang et al. (2014); Wu et al. (2015) provide a promising direction, but their model complexity and high demands of training data prohibit their broader usages in clinical cohort studies.

Figure 2: (a) Subspace differences, w.r.t. the index i; (b) Convergence of Algorithm 1, w.r.t. the index i; (c) Iteration-wise subspace differences, w.r.t. the index i.

Figure 3: (a) Predicted weight vs true weight for task 1; (b) Predicted weight vs true weight for task 2; (c) Distribution of correlation between predicted weight and true weight for all tasks

Comparison of subspace differences for each layer of SN, f-MLP, and rf-MLP.

Table content: <table><tr><td>Metric</td><td colspan="3">Subspace Difference</td><td colspan="3">Maximum Mutual Coherence</td><td colspan="3">Mean Mutual Coherence</td></tr><tr><td>Method</td><td>SN</td><td>f-MLP</td><td>rf-MLP</td><td>SN</td><td>f-MLP</td><td>rf-MLP</td><td>SN</td><td>f-MLP</td><td>rf-MLP</td></tr><tr><td>Layer 1</td><td>0.0313</td><td>0.0315</td><td>0.0317</td><td>0.7608</td><td>0.7727</td><td>0.7895</td><td>0.2900</td><td>0.2725</td><td>0.2735</td></tr><tr><td>Layer 2</td><td>0.0321</td><td>0.0321</td><td>0.0321</td><td>0.8283</td><td>0.7603</td><td>0.7654</td><td>0.2882</td><td>0.2820</td><td>0.2829</td></tr><tr><td>Layer 3</td><td>0.0312</td><td>0.0315</td><td>0.0313</td><td>0.8493</td><td>0.7233</td><td>0.7890</td><td>0.2586</td><td>0.2506</td><td>0.2485</td></tr></table>

Average normalized mean square error under different approaches for synthetic data. Standard deviation of 10 trials is given in parenthesis.

Table content: <table><tr><td>Percent</td><td>Uncensored (LS + 1 )</td><td>Single Task (Shallow) Censored (LS + 1)</td><td>Nonlinear Censored (Tobit)</td><td colspan="3">Multi Task (Shallow) Uncensored (Multi Trace) Censored (Multi Trace)</td></tr><tr><td>40%</td><td>0.1412 (0.0007)</td><td>0.1127 (0.0010)</td><td>0.0428 (0.0003)</td><td colspan="2">0.1333 (0.0009)</td><td>0.1053 (0.0027)</td></tr><tr><td>50%</td><td>0.1384 (0.0005)</td><td>0.1102 (0.0010)</td><td>0.0408 (0.0004)</td><td colspan="2">0.1323 (0.0010)</td><td>0.1054 (0.0042)</td></tr><tr><td>60%</td><td>0.1365 (0.0005)</td><td>0.1088 (0.0009)</td><td>0.0395 (0.0003)</td><td colspan="2">0.1325 (0.0012)</td><td>0.1031 (0.0046)</td></tr><tr><td>70%</td><td>0.1349 (0.0005)</td><td>0.1078 (0.0010)</td><td>0.0388 (0.0004)</td><td colspan="2">0.1315 (0.0013)</td><td>0.1024 (0.0042)</td></tr><tr><td>80%</td><td>0.1343 (0.0011)</td><td>0.1070 (0.0012)</td><td>0.0383 (0.0006)</td><td colspan="2">0.1308 (0.0008)</td><td>0.1040 (0.0011)</td></tr><tr><td>Percent</td><td>DNN i (naive)</td><td colspan="2">Deep Neural Network DNN ii (censored) DNN iii (censored + low-rank)</td><td>Layer 1</td><td colspan="2">Subspace Net (SN)</td><td>Layer 3</td></tr><tr><td>40%</td><td>0.0623 (0.0041)</td><td>0.0489 (0.0035)</td><td>0.0431 (0.0041)</td><td colspan="2">0.0390 (0.0004)</td><td>0.0369 (0.0002)</td></tr><tr><td>50%</td><td>0.0593 (0.0048)</td><td>0.0462 (0.0042)</td><td>0.0400 (0.0039)</td><td colspan="2">0.0389 (0.0007)</td><td>0.0366 (0.0003)</td></tr><tr><td>60%</td><td>0.0587 (0.0053)</td><td>0.0455 (0.0054)</td><td>0.0395 (0.0050)</td><td colspan="2">0.0388 (0.0006)</td><td>0.0364 (0.0003)</td></tr><tr><td>70%</td><td>0.0590 (0.0071)</td><td>0.0447 (0.0043)</td><td>0.0386 (0.0058)</td><td colspan="2">0.0388 (0.0006)</td><td>0.0363 (0.0003)</td></tr><tr><td>80%</td><td>0.0555 (0.0057)</td><td>0.0431 (0.0053)</td><td>0.0380 (0.0057)</td><td colspan="2">0.0390 (0.0008)</td><td>0.0364 (0.0005)</td></tr></table>

Average normalized mean square error at each layer for subspace network (R = 10) for synthetic data. Standard deviation of 10 trials is given in parenthesis.

Table content: <table><tr><td>Percent</td><td>Layer 1</td><td>Layer 2</td><td>Layer 3</td><td>Layer 10</td><td>Layer 20</td></tr><tr><td>40%</td><td>0.0390 (0.0004)</td><td>0.0381 (0.0005)</td><td>0.0369 (0.0002)</td><td>0.0368 (0.0002)</td><td>0.0368 (0.0002)</td></tr><tr><td>50%</td><td>0.0389 (0.0007)</td><td>0.0379 (0.0005)</td><td>0.0366 (0.0003)</td><td>0.0366 (0.0003)</td><td>0.0365 (0.0003)</td></tr><tr><td>60%</td><td>0.0388 (0.0006)</td><td>0.0378 (0.0004)</td><td>0.0364 (0.0003)</td><td>0.0364 (0.0003)</td><td>0.0363 (0.0003)</td></tr><tr><td>70%</td><td>0.0388 (0.0006)</td><td>0.0378 (0.0005)</td><td>0.0363 (0.0003)</td><td>0.0363 (0.0003)</td><td>0.0362 (0.0003)</td></tr><tr><td>80%</td><td>0.0390 (0.0008)</td><td>0.0378 (0.0006)</td><td>0.0364 (0.0005)</td><td>0.0363 (0.0005)</td><td>0.0363 (0.0005)</td></tr></table>

Average normalized mean square error under different approaches for real data. Standard deviation of 10 trials is given in parenthesis. on studies investigating brain and fluid biomarkes of the disease, and supported the long running project called Alzheimer's Disease Neuroimaging Initiative (ADNI) from 2003. We used the ADNI-1 cohort (http://adni.loni.usc.edu/). In the experiments, we used the 1.5 Tesla structural MRI collected at the baseline, and performed cortical reconstruction and volumetric segmentations with the FreeSurfer following the procotol in Jack et al.(2008). For each MRI image, we extracted 138 features representing the cortical thickness and surface areas of region-of-interests (ROIs) using the Desikan-Killiany cortical atlasDesikan et al. (2006). After preprocessing, we obtained a dataset containing 670 samples and 138 features. These imaging features were used to predict a set of 30 clinical scores including ADAS scoresRosen et al. (1984) at baseline and future (6 months from baseline), baseline Logical Memory from Wechsler Memory Scale IV Scale-Fourth (2009), Neurobattery scores (i.e. immediate recall total score and Rey Auditory Verbal Learning Test scores), and the Neuropsychiatric Inventory Cummings (1997) at baseline and future.

Table content: <table><tr><td>Percent</td><td>Least Square</td><td>Single Task (Censored) LS + 1</td><td>Tobit</td><td colspan="2">Multi Task (Censored) Multi Trace Multi 2,1</td></tr><tr><td>40%</td><td>0.3874 (0.0203)</td><td>0.2393 (0.0056)</td><td>0.3870 (0.0306)</td><td>0.2572 (0.0156)</td><td>0.2006 (0.0099)</td></tr><tr><td>50%</td><td>0.3119 (0.0124)</td><td>0.2202 (0.0049)</td><td>0.3072 (0.0144)</td><td>0.2406 (0.0175)</td><td>0.2002 (0.0132)</td></tr><tr><td>60%</td><td>0.2779 (0.0123)</td><td>0.2112 (0.0055)</td><td>0.2719 (0.0114)</td><td>0.2596 (0.0233)</td><td>0.2072 (0.0204)</td></tr><tr><td>70%</td><td>0.2563 (0.0108)</td><td>0.2037 (0.0042)</td><td>0.2516 (0.0108)</td><td>0.2368 (0.0362)</td><td>0.2017 (0.0116)</td></tr><tr><td>80%</td><td>0.2422 (0.0112)</td><td>0.2005 (0.0054)</td><td>0.2384 (0.0099)</td><td>0.2176 (0.0171)</td><td>0.2009 (0.0050)</td></tr><tr><td>Percent</td><td>DNN i (naive)</td><td colspan="2">Deep Neural Network DNN ii (censored) DNN iii (censored + low-rank)</td><td colspan="2">Subspace Net (SN) Layer 1 Layer 3</td></tr><tr><td>40%</td><td>0.2549 (0.0442)</td><td>0.2388 (0.0121)</td><td>0.2113 (0.0063)</td><td>0.2016 (0.0057)</td><td>0.1981 (0.0031)</td></tr><tr><td>50%</td><td>0.2236 (0.0066)</td><td>0.2208 (0.0062)</td><td>0.2127 (0.0118)</td><td>0.1992 (0.0040)</td><td>0.1971 (0.0038)</td></tr><tr><td>60%</td><td>0.2215 (0.0076)</td><td>0.2200 (0.0076)</td><td>0.2087 (0.0102)</td><td>0.1990 (0.0061)</td><td>0.1967 (0.0038)</td></tr><tr><td>70%</td><td>0.2149 (0.0077)</td><td>0.2141 (0.0079)</td><td>0.2093 (0.0137)</td><td>0.1981 (0.0046)</td><td>0.1953 (0.0039)</td></tr><tr><td>80%</td><td>0.2132 (0.0138)</td><td>0.2090 (0.0079)</td><td>0.2069 (0.0135)</td><td>0.1970 (0.0034)</td><td>0.1956 (0.0040)</td></tr><tr><td/><td/><td/><td colspan="3">Table 4: ANMSE for non-calibrated vs. cal-</td></tr><tr><td/><td/><td/><td colspan="3">ibrated SN for real data (6 layers). Standard</td></tr><tr><td/><td/><td/><td colspan="3">deviation of 10 trials is given in parenthesis.</td></tr><tr><td/><td/><td/><td>Percent</td><td>Non-calibrate</td><td>Calibrate</td></tr><tr><td/><td/><td/><td>40%</td><td>0.1993 (0.0034)</td><td>0.1977 (0.0031)</td></tr><tr><td/><td/><td/><td>50%</td><td>0.1987 (0.0043)</td><td>0.1967 (0.0036)</td></tr><tr><td/><td/><td/><td>60%</td><td>0.1991 (0.0044)</td><td>0.1964 (0.0039)</td></tr><tr><td/><td/><td/><td>70%</td><td>0.1982 (0.0042)</td><td>0.1951 (0.0038)</td></tr><tr><td/><td/><td/><td>80%</td><td>0.1984 (0.0041)</td><td>0.1954 (0.0039)</td></tr></table>

Due to the nonconvexity of the online subspace iterates, it is challenging to directly analyze how fast the online cumulative loss Ct approaches the optimal batch cost C t . AsShen et al. (2016) advocates, we instead investigate whether Ĉt converges to Ct . That is established by first referring to the Lemma 2 ofMardani et al. (2013): the distance between successive subspace estimates will vanish as fast as o(1/t): U t -U t-1 F ≤ B t , for some constant B that is independent of t and N . Next, following the proof of Proposition 2 inShen et al. (2016), we can similarly show that: if {U t } T t=1 and {V t x t } T t=1 are uniformly bounded, i.e., U t F ≤ B 1 , and V t x t 2 ≤ B 2 , ∀t ≤ T , then with constants B 1 , B 2 > 0 and by choosing a constant step size µ t = µ, we have a bounded regret as:

Table content: <table><tr><td>R T ≤</td><td>B 2 (ln(T ) + 1) 2 2µT</td><td>+</td><td>5B 2 6µT</td><td>.</td></tr><tr><td colspan="3">Algorithm 3 Gradient descent algorithm for problem P:V.</td><td/><td/></tr></table>

Running time and platforms for different methods for synthetic data.Computation speed. All experiments run on the same machine (1 x Six-core Intel Xeon E5-1650 v3 [3.50GHz], 12 logic cores, 128 GB RAM). GPU accelerations are enabled for DNN baselines, while SN has not exploited the same accelerations yet. The running time for a single round training on synthetic data (N=5000, D=200, T=100) is given in Table6. Training each layer of SN will cost 109 seconds on average. As we can see, SN improves generalization performance without significant loss of speed. Furthermore, we will explore to accelerate SN further, by reading data in batch mode and performing parallel updates.

Table content: <table><tr><td>Method</td><td>Time (s)</td><td>Platform</td></tr><tr><td>Least Square</td><td>0.02</td><td>Matlab</td></tr><tr><td>LS+ 2</td><td>0.02</td><td>Matlab</td></tr><tr><td>LS+ 1</td><td>18.4</td><td>Matlab</td></tr><tr><td>Multi-trace</td><td>32.3</td><td>Matlab</td></tr><tr><td>Multi-21</td><td>27.0</td><td>Matlab</td></tr><tr><td>Censor</td><td>1680</td><td>Matlab</td></tr><tr><td>SN (per layer)</td><td>109</td><td>Python</td></tr><tr><td>DNN</td><td>659</td><td>Tensorflow</td></tr></table>

Average normalized mean square error at each layer for subspace network (R = 5) for real data. Standard deviation of 10 trials is given in parenthesis.

Table content: <table><tr><td>Percent</td><td>Layer 1</td><td>Layer 2</td><td>Layer 3</td><td>Layer 5</td><td>Layer 10</td></tr><tr><td>40%</td><td colspan="5">0.2016 (0.0057) 0.2000 (0.0039) 0.1981 (0.0031) 0.1977 (0.0031) 0.1977 (0.0031)</td></tr><tr><td>50%</td><td colspan="5">0.1992 (0.0040) 0.1992 (0.0053) 0.1971 (0.0038) 0.1968 (0.0036) 0.1967 (0.0035)</td></tr><tr><td>60%</td><td colspan="5">0.1990 (0.0061) 0.1990 (0.0047) 0.1967 (0.0038) 0.1964 (0.0039) 0.1964 (0.0038)</td></tr><tr><td>70%</td><td colspan="5">0.1981 (0.0046) 0.1966 (0.0052) 0.1953 (0.0039) 0.1952 (0.0039) 0.1951 (0.0038)</td></tr><tr><td>80%</td><td colspan="5">0.1970 (0.0034) 0.1967 (0.0044) 0.1956 (0.0040) 0.1955 (0.0039) 0.1953 (0.0039)</td></tr></table>

Average normalized mean square error under different rank assumptions for real data. Standard deviation of 10 trials is given in parenthesis.

Table content: <table><tr><td>Method</td><td>Percent -Rank</td><td>R = 1</td><td>R = 3</td><td>R = 5</td><td>R = 10</td></tr><tr><td/><td>40%</td><td>0.2052 (0.0030)</td><td>0.1993 (0.0036)</td><td>0.1983 (0.0033)</td><td>0.2010 (0.0044)</td></tr><tr><td/><td>50%</td><td>0.2047 (0.0029)</td><td>0.1983 (0.0034)</td><td>0.1981 (0.0043)</td><td>0.2001 (0.0046)</td></tr><tr><td>SN</td><td>60%</td><td>0.2052 (0.0033)</td><td>0.1988 (0.0047)</td><td>0.1978 (0.0044)</td><td>0.1996 (0.0052)</td></tr><tr><td/><td>70%</td><td>0.2043 (0.0044)</td><td>0.1975 (0.0042)</td><td>0.1966 (0.0038)</td><td>0.1990 (0.0057)</td></tr><tr><td/><td>80%</td><td>0.2058 (0.0051)</td><td>0.1977 (0.0042)</td><td>0.1965 (0.0046)</td><td>0.1990 (0.0058)</td></tr><tr><td/><td>40%</td><td>0.2322 (0.0146)</td><td>0.2360 (0.0060)</td><td>0.2113 (0.0063)</td><td>0.2196 (0.0124)</td></tr><tr><td>DNN iii (censored + low-rank)</td><td>50% 60% 70%</td><td>0.2298 (0.0093) 0.2244 (0.0132) 0.2178 (0.0129)</td><td>0.2256 (0.0127) 0.2277 (0.0099) 0.2177 (0.0115)</td><td>0.2127 (0.0118) 0.2087 (0.0102) 0.2093 (0.0137)</td><td>0.2235 (0.0142) 0.2145 (0.0208) 0.2083 (0.0127)</td></tr><tr><td/><td>80%</td><td>0.2256 (0.0117)</td><td>0.2250 (0.0079)</td><td>0.2069 (0.0135)</td><td>0.2158 (0.0183)</td></tr></table>

Peer review (iclr2018):
 This work proposes a multi task learning framework for the modeling of clinical data in neurodegenerative diseases. Differently from previous applications of machine learning in neurodegeneration modeling, the proposed approach models the clinical data accounting for the bounded nature of cognitive tests scores. The framework is represented by a feed-forward deep architecture analogous to a residual network. At each layer a low-rank constraint is enforced on the linear transformation, while the cost function is specified in order to differentially account for the bounds of the predicted variables. The idea of explicitly accounting for the boundedness of clinical scores is interesting, although the assumption of the proposed model is still incorrect: clinical scores are defined on discrete scales. For this reason the Gaussian assumption for the cost function used in the method is still not appropriate for the proposed application. Furthermore, while being the main methodological drive of this work, the paper does not show evidence about improved predictive performance and generalisation when accounting for the boundedness of the regression targets. The proposed algorithm is also generally compared with respect to linear methods, and the authors could have provided a more rigorous benchmark including standard non-linear prediction approaches (e.g. random forests, NN, GP, …). Overall, the proposed methods seems to provide little added value to the large amount of predictive methods proposed so far for prediction in neurodegenerative disorders. Moreover, the proposed experimental paradigm appears flawed. What is the interest of predicting baseline (or 6 months at best) cognitive scores (relatively low-cost and part of any routine clinical assessment) from brain imaging data (high-cost and not routine)? Other remarks. - In section 2.2 and 4 there is some confusion between iteration indices and samples indices “i”. - Contrarily to what is stated in the introduction, the loss functions proposed in page 3 (first two formulas) only accounts for the lower bound of the predicted variables. - Figure 2, synthetic data. The scale of the improvement of the subspace difference is quite tiny, in the order of 1e-2 when compared to U, and of 1e-5 across iterations. The loss function of Figure 2.b also does not show a strong improvement across iterations, while indicating a rather large instability of the optimisation procedure. These aspects may be a sign of convergence issues. - The dimensionality of the subspace representation importantly depends on the choice of the rank R of U and V. This is a crucial parameters that is however not discussed nor analysed in the paper. - The synthetic example of page 7 is quite misleading and potentially biased towards the proposed model. The authors are generating the synthetic data according to the model, and it is thus not surprising that they managed to obtain the best performance. In particular, due to the nonlinear nature of (1), all the competing linear models are expected to perform poorly in this kind of setting. - The computation time for the linear model shown in Table 3 is quite surprising (~20 minutes for linear regression of 5k observations). Is there anything that I am missing?