Individualized Multilayer Tensor Learning with An Application in Imaging Analysis

TL;DR

This paper introduces a novel multilayer tensor learning approach tailored for heterogeneous multimodal imaging data, effectively capturing spatial features and improving disease prediction accuracy in breast cancer imaging.

Contribution

It develops an individualized multilayer tensor decomposition method that incorporates heterogeneity and multimodality, with scalable algorithms and theoretical guarantees.

Findings

Outperforms existing methods in simulated data

Achieves higher prediction accuracy in breast cancer imaging

Provides theoretical analysis of algorithm convergence and consistency

Abstract

This work is motivated by multimodality breast cancer imaging data, which is quite challenging in that the signals of discrete tumor-associated microvesicles (TMVs) are randomly distributed with heterogeneous patterns. This imposes a significant challenge for conventional imaging regression and dimension reduction models assuming a homogeneous feature structure. We develop an innovative multilayer tensor learning method to incorporate heterogeneity to a higher-order tensor decomposition and predict disease status effectively through utilizing subject-wise imaging features and multimodality information. Specifically, we construct a multilayer decomposition which leverages an individualized imaging layer in addition to a modality-specific tensor structure. One major advantage of our approach is that we are able to efficiently capture the heterogeneous spatial features of signals that are not characterized by a population structure as well as integrating multimodality information simultaneously. To achieve scalable computing, we develop a new bi-level block improvement algorithm. In theory, we investigate both the algorithm convergence property, tensor signal recovery error bound and asymptotic consistency for prediction model estimation. We also apply the proposed method for simulated and human breast cancer imaging data. Numerical results demonstrate that the proposed method outperforms other existing competing methods.

Tables

Table 1. Table 1: Prediction results (standard errors in subscripts) of the proposed method (IMTL) compared with five other competing methods for simulation study, based on 100 replications. Equal sample size for training set ( N t ​ r subscript 𝑁 𝑡 𝑟 N_{tr} ), validation set ( N v ​ l subscript 𝑁 𝑣 𝑙 N_{vl} ) and testing set ( N t ​ s subscript 𝑁 𝑡 𝑠 N_{ts} ).

${𝑵}_{𝒕𝒓}$		VPL	TR	MPCA	HOCPD	CNN	IMTL
$60$	Pred Accuracy	${0.593}_{0.051}$	${0.573}_{0.054}$	${0.598}_{0.078}$	${0.731}_{0.175}$	${0.714}_{0.083}$	${0.858}_{0.069}$
	Sensitivity	${0.000}_{0.000}$	${0.415}_{0.103}$	${0.443}_{0.143}$	${0.526}_{0.229}$	${0.379}_{0.209}$	${0.687}_{0.119}$
	Specificity	${1.000}_{0.000}$	${0.674}_{0.097}$	${0.691}_{0.080}$	${0.915}_{0.094}$	${0.955}_{0.023}$	${0.985}_{0.064}$
$100$	Pred Accuracy	${0.608}_{0.043}$	${0.603}_{0.045}$	${0.634}_{0.057}$	${0.808}_{0.168}$	${0.884}_{0.069}$	${0.941}_{0.041}$
	Sensitivity	${0.001}_{0.005}$	${0.453}_{0.090}$	${0.490}_{0.113}$	${0.682}_{0.226}$	${0.774}_{0.160}$	${0.878}_{0.064}$
	Specificity	${0.998}_{0.008}$	${0.697}_{0.091}$	${0.728}_{0.070}$	${0.931}_{0.064}$	${0.985}_{0.010}$	${0.993}_{0.030}$

Table 2. Table 2: Prediction results of the proposed method (IMTL) compared with five other competing methods for human breast cancer imaging data, based on 30 random replications. Sample size for the training set, the validation set and the testing set are 60, 40 and 60, respectively.

Model	VPL	TR	MPCA	HOCPD	CNN	IMTL
Pred Accuracy	0.681	0.656	0.766	0.803	0.871	0.921
Sensitivity	0.086	0.628	0.535	0.656	0.729	0.903
Specificity	0.962	0.671	0.875	0.878	0.928	0.925

Equations

$${\mathcal{X}}_i={\mathbf{\Theta }}_i+{\mathcal{N}}_i,i=1,\dots ,N,$$

$$\underset{\mathbf{\beta }}{\min }\frac{1}{N}\sum _{i=1}^N\mathcal{L}\left(y_i,f({\mathcal{X}}_i),\mathbf{\beta }\right),$$

$${\mathcal{X}}_i\approx \sum _{r=1}^R{w}_{ir}{\mathcal{B}}_r,i=1,\dots ,N,$$

$$\mathbf{\Theta }=\mathbf{W}\circ {\mathbf{B}}^1\circ \cdots \circ {\mathbf{B}}^D=\sum _{r=1}^R{\mathbf{W}}_{\cdot r}\circ {\mathbf{B}}_{\cdot r}^1\circ \cdots \circ {\mathbf{B}}_{\cdot r}^D,$$

$${\mathcal{X}}_i={\mathbf{\Theta }}_i+{\mathcal{N}}_i=\sum _{r=1}^R{\bar{w}}_{ir}{\mathcal{B}}_r+{\mathcal{N}}_i,$$

$$\bm{\Theta}_{i}=\sum\limits_{r=1}^{R}\bar{w}_{ir}\bm{\mathcal{B}}_{r}+\bm{\mathcal{S}}_{i},\quad\mbox{ s.t.}\quad\big{\langle}\bm{\mathcal{B}}_{r},\bm{\mathcal{S}}_{i}\big{\rangle}=0,\;1\leq r\leq R,$$

$$\min_{\bm{\bm{w}}_{i},\bm{\mathcal{B}}_{r},\mu_{i},\bar{\bm{\mathcal{S}}}_{i}}\sum_{i=1}^{N}\|\bm{\mathcal{X}}_{i}-\sum\limits_{r=1}^{R}\bar{w}_{ir}\bm{\mathcal{B}}_{r}-\mu_{i}\bar{\bm{\mathcal{S}}}_{i}\|^{2}_{F},\quad\mbox{ s.t.}\quad\big{\langle}\bm{\mathcal{B}}_{r},\bar{\bm{\mathcal{S}}}_{i}\big{\rangle}=0,\;1\leq r\leq R.$$

$${\mathbf{\Theta }}^{(m)}=\mathbf{W}\circ {\mathbf{B}}^{(m),1}\circ \cdots \circ {\mathbf{B}}^{(m),D}+\mathcal{S},$$

$$\bm{\Theta}^{(m)}_{i}=\sum\limits_{r=1}^{R}\bar{w}^{(m)}_{ir}\bm{\mathcal{B}}^{(m)}_{r}+\mu_{i}\bar{\bm{\mathcal{S}}}_{i},\quad\mbox{s.t.}\quad\big{\langle}\bm{\mathcal{B}}^{(m)}_{r},\bar{\bm{\mathcal{S}}}_{i}\big{\rangle}=0,\;1\leq r\leq R,$$

$${\mathbf{\Theta }}_i=\sum _{r=1}^R{w}_{ir}{\mathbf{b}}_r^1\circ \cdots \circ {\mathbf{b}}_r^D\circ {\mathbf{b}}_r^{modal}+{\mathbf{s}}_i^1\circ \cdots \circ {\mathbf{s}}_i^D\circ {\mathbf{s}}_i^{modal},$$

$${\mathbf{\Theta }}_i^{(m)}=\sum _{r=1}^R{w}_{ir}{b}_{mr}^{modal}{\mathbf{b}}_r^1\circ \cdots \circ {\mathbf{b}}_r^D+s_{im}^{modal}{\mathbf{s}}_i^1\circ \cdots \circ {\mathbf{s}}_i^D.$$

$$L\left(\mathbf{W},\{{\mathbf{B}}^{(m),d}{\}}_{m^{\prime }s,d^{\prime }s},\{{\mathbf{s}}_i^d{\}}_{i^{\prime }s,d^{\prime }s}|\mathcal{X}\right)=\sum _{m=1}^M\Vert {\mathcal{X}}^{(m)}-{\mathbf{\Theta }}^{(m)}{\Vert }_F^2+{\lambda }_s\sum _{m,i,r}⟨{\mathcal{B}}_r^{(m)},{\mathcal{S}}_i{⟩}^2.$$

$$\hat{\bm{W}}^{[new]}=\arg\min_{\bm{W}}\sum_{m=1}^{M}\left\|\tilde{\bm{\mathcal{X}}}^{(m)}_{(1)}-\bm{W}\big{(}\hat{\bm{B}}^{(m),D}\odot\cdots\odot\hat{\bm{B}}^{(m),1}\big{)}^{T}\right\|^{2}_{F},$$

$$L^{(m)}\big{(}\bm{B}^{(m),d}|\{\bm{B}^{(m),d^{\prime}}\}_{d^{\prime}\neq d}\big{)}=\left\|\tilde{\bm{\mathcal{X}}}^{(m)}_{(d+1)}-\bm{B}^{(m),d}\big{(}\bm{B}_{-d}^{(m)}\odot\hat{\bm{W}}\big{)}^{T}\right\|^{2}_{F}+\lambda_{s}\sum_{i}\Big{\langle}\bm{B}^{(m),d},\;\hat{\bm{\mathcal{S}}}_{i,(d)}\bm{B}_{-d}^{(m)}\Big{\rangle}^{2}$$

$$\begin{array}{rl}& {\mathbf{B}}^{\ast (m),d}=\arg \underset{{\mathbf{B}}^{(m),d}}{\min }{L}^{(m)}({\mathbf{B}}^{(m),d}|\{{\mathbf{B}}^{(m),d^{\prime },[t]}{\}}_{d^{\prime }\ne d}),\\ & d_m^{\ast }=\arg \underset{1\le d\le D}{\min }{L}^{(m)}({\mathbf{B}}^{\ast (m),d}|\{{\mathbf{B}}^{(m),d^{\prime },[t]}{\}}_{d^{\prime }\ne d}),\end{array}$$

$$L_{i}\big{(}\bm{s}_{i}^{d}\;|\;\{\bm{s}_{i}^{d^{\prime}}\}_{d^{\prime}\neq d}\big{)}=\sum_{m=1}^{M}\|\bar{\bm{\mathcal{X}}}^{(m)}_{i,(d)}-\bm{s}_{i}^{d}(\bm{S}_{i,-d})^{T}\|^{2}_{F}+\lambda_{s}\sum_{m,r}\langle\bm{s}_{i}^{d},\;\hat{\bm{\mathcal{B}}}^{(m)}_{r,(d)}\bm{S}_{i,-d}\rangle^{2}$$

$${\mathbf{s}}_i^{\ast d}=\arg \underset{{\mathbf{s}}_i^d}{\min }{L}_i({\mathbf{s}}_i^d|\{{\mathbf{s}}_i^{d^{\prime },[t]}{\}}_{d^{\prime }\ne d}),\text{and}{d}_i^{\ast }=\arg \underset{1\le d\le D}{\min }{L}_i({\mathbf{s}}_i^{\ast d}|\{{\mathbf{s}}_i^{d^{\prime },[t]}{\}}_{d^{\prime }\ne d}),$$

$${\mathbf{\Theta }}_{[1:n]}=\sum _{r=1}^R{\mathbf{w}}_r^{[1:n]}\circ {\mathbf{B}}_{\cdot r}^1\circ \cdots \circ {\mathbf{B}}_{\cdot r}^D+\sum _{i=1}^n{\tilde{\mathbf{u}}}_i^{[1:n]}\circ {\mathbf{s}}_i^1\circ \cdots \circ {\mathbf{s}}_i^D,$$

$$\sum _{d=1}^D{\mathcal{K}}_{{\tilde{\mathbf{B}}}_{[1:n]}^d}\ge 2R+n+D$$

$$l(\mathbf{\Theta }|X_{i,j_1\cdots j_D}^{(m)}=x_{i,j_1\cdots j_D}^{(m)})=(x_{i,j_1\cdots j_D}^{(m)}-{\theta }_{i,j_1\cdots j_D}^{(m)}{)}^2=(x_{i,j_1\cdots j_D}^{(m)}-\sum _{r=1}^R{w}_{ir}{b}_{j_{1}r}^{(m),1}\cdots b_{j_{D}r}^{(m),D}-s_{i{j}_1}^1\cdots s_{i{j}_D}^D{)}^2.$$

$$L(\mathbf{\Theta }|\mathcal{X})=\sum _{i=1}^N\sum _{m=1}^M\sum _{j_1=1}^{p_1}\cdots \sum _{j_D=1}^{p_D}l(\mathbf{\Theta }|x_{i,j_1\cdots j_D}^{(m)})+{\lambda }_{|\mathcal{I}|}p(\mathbf{\Theta }),$$

$$\hat{\mathbf{\Theta }}=\arg \underset{\mathbf{\Theta }\in \mathcal{S}}{\min }L(\mathbf{\Theta }).$$

$$K(\mathbf{\Theta },{\mathbf{\Theta }}_0)=\frac{1}{NM{p}_1\cdots p_D}\sum _{i=1}^N\sum _{m=1}^M\sum _{j_1=1}^{p_1}\cdots \sum _{j_D=1}^{p_D}\mathbf{E}\{l_{\Delta }(\mathbf{\Theta }|X_{i,j_1\cdots j_D}^{(m)})\},$$

$$V(\mathbf{\Theta },{\mathbf{\Theta }}_0)=\frac{1}{NM{p}_1\cdots p_D}\sum _{i=1}^N\sum _{m=1}^M\sum _{j_1=1}^{p_1}\cdots \sum _{j_D=1}^{p_D}\text{mbox}Var\{l_{\Delta }(\mathbf{\Theta }|X_{i,j_1\cdots j_D}^m)\}.$$

$$\mathbf{P}(\frac{1}{|\mathcal{I}{|}^{1/2}}\Vert \hat{\mathbf{\Theta }}-{\mathbf{\Theta }}_0{\Vert }_F\ge {\tau }_{|\mathcal{I}|})\le 7\exp (-c_1|\mathcal{I}|{\tau }_{|\mathcal{I}|}^2)$$

$$G_N(\mathbf{y},F({\mathbf{\Theta }}_{|\mathcal{I}|}),\mathbf{\beta })=\frac{1}{N}\sum _{i=1}^N\mathcal{L}(y_i,f({\mathbf{\Theta }}_i),\mathbf{\beta })+{\lambda }_{\mathbf{\beta }}\cdot p^{\ast }(\mathbf{\beta }),$$

$$|\mathcal{L}(y_i,{\eta }_i)-\mathcal{L}(y_i,{\tilde{\eta }}_i)|\le K_i(y_i)|{\eta }_i-{\tilde{\eta }}_i|,1\le i\le N,$$

$$\bigg{|}G_{N}(\bm{y},F(\bm{\mathcal{X}}_{|\mathcal{I}|}),\bm{\beta})-G_{N}(\bm{y},F(\bm{\Theta}_{|\mathcal{I}|}),\bm{\beta})\bigg{|}\rightarrow_{p}0.$$

$$|\hat{\mathbf{\beta }}({\mathcal{X}}_{|\mathcal{I}|})-\tilde{\mathbf{\beta }}({\mathbf{\Theta }}_{|\mathcal{I}|})|{\to }_{p}0\text{and}\hat{\mathbf{\beta }}({\mathcal{X}}_{|\mathcal{I}|}){\to }_p{\mathbf{\beta }}_0,$$

$$\big{|}\mathcal{R}(\hat{\mathcal{L}}_{N})-\mathcal{R}(\tilde{\mathcal{L}}_{N})\big{|}\rightarrow_{p}0.$$

Full text

Key words: Cancer imaging, dimension reduction, heterogeneous modeling, high-order tensor decomposition, multimodality integration, spatial information.

1 Introduction

In recent years imaging analysis has encountered explosive growth due to high demand and applications in biomedical imaging for diagnosing disease status and assessing treatment outcomes. The biomedical applications of imaging analyses are especially powerful in cancer radiotherapy and neuroimaging (Fass, 2008; Martino et al., 2008; Caffo et al., 2010; Zhou et al., 2013; Li et al., 2015). In addition, new advanced technologies in imaging analysis can be utilized to better understand the structures of the human body associated with biological, psychological and clinical traits, (Lindquist, 2008; Lazar, 2008; Friston, 2009; Hinrichs et al., 2009; Kang et al., 2012), so cancer and chronic diseases can be diagnosed earlier and intervention treatments can be implemented.

This paper is motivated by utilizing multimodality optical imaging data for diagnosing early-stage breast cancer prior to tumor formation. Figure 1 illustrates four-modality optical imaging data for cancerous and normal rat breast tissues, showing that a large number of tumor-associated microvesicles (TMVs; circled in red) appear spatially aligned in the cancerous tissue. One unique aspect of cancer imaging is that individuals might have imaging at various locations, and TMVs as an important biomarker are also randomly located leading to heterogeneous signal patterns. This is quite different from the brain imaging where the same regions of the brain usually share similar characteristics and functions over different subjects. Moreover, TMVs have relatively weak signals in terms of smaller size and pixel contrast values compared to modality-specific background and noise. Both great heterogeneity and weak strength make the signals difficult to be captured by a conventional population model which mainly characterizes the population-wise variation. In addition, it is also crucial to integrate information from multiple imaging modalities effectively.

In contrast to standard vector data, imaging data are usually in the form of a multi-dimensional array (also known as a multi-order tensor), which preserves higher-order structures containing rich spatial information. Traditional statistical models mostly treat covariates as vectors. However, fitting a model with the vectorized imaging array becomes infeasible when the array size is large. For instance, a four-modality breast cancer imaging size of $3,600\times 3,600$ for each modality implicitly requires $4\times 3,600^{2}=51,840,000$ regression parameters. This leads to an ultra-high dimensional problem which could be unestimable even with additional regularization techniques. Most importantly, vectorization of an array is not capable of preserving the important spatial structure of imaging. Under the vectorization framework, Bayesian variable selection approaches have been developed for high-dimensional imaging regression models to identify important regions, by applying Markov random field priors to account for the spatial correlation between voxels (Penny et al., 2005; Bowman, 2007; Bowman et al., 2008; Derado et al., 2010; Li et al., 2015). However, this would rapidly increase the complexity of the prior when the tensor order increases, and could be computationally infeasible if the tensor order is high.

Alternatively, functional data analysis can be adopted to construct a two-dimensional image predictor (Reiss and Ogden, 2010) in a functional regression model. Nevertheless, the extension to three dimensions or beyond could be impractical due to high-dimensional parameters arising from higher-order imaging data (Zhou et al., 2013). Recent developments in imaging analyses employ a two-stage strategy, that is, perform a dimension reduction such as principle component analysis (PCA) first, and then fit a regression model based on the extracted principal components (Caffo et al., 2010), however, which still requires conversion of the image into a vector first.

To utilize the unique higher-order structure of the image covariates, Zhou et al. (2013) and Li et al. (2013) propose a tensor regression model, where the coefficients associated with imaging voxels are formulated as a tensor and assumed to be low-rank. But this potentially requires voxel-wise registration over different individuals, which is not applicable for cancer imaging. In addition, the tensor decomposition technique, such as the CANDECOMP/PARAFAC (CP) method (Beckmann and Smith, 2005), is applied for processing signals from the image covariates directly (Cong et al., 2012; Mørup et al., 2006; Karahan et al., 2015; Lee et al., 2007; Miranda et al., 2018). Due to its population-wise low-rank constraint, it is not capable of capturing heterogeneous imaging information.

Another prevalent tool in image recognition and classification is the deep learning method, in particular the convolutional neural network (CNN, Ciresan et al., 2011; Krizhevsky et al., 2012). In general, the CNN is powerful for image classification, which utilizes multiple hidden layers with a variety of convolution, pooling and activation methods. However, due to its complex architecture and the large number of parameters involved, the CNN usually requires massive training data to guarantee good performance (Keshari et al., 2018; Wagner et al., 2013; Abbasi-Asl and Yu, 2017a, b), especially when the imaging signals are heterogeneous and relatively weak. In our motivated cancer imaging study, such a large sample size for homogeneous samples is infeasible.

In this paper, we propose an individualized multilayer tensor learning (IMTL) utilizing additional individualized layers to incorporate heterogeneity to a higher-order tensor decomposition, which captures extra variation of the imaging covariates that cannot be characterized by a population structure. Furthermore, we extend it to multimodality data to integrate individual-specific information over different modalities and incorporate modality-specific features through different layers. We develop a scalable and efficient algorithm which can be implemented with parallel computing. In addition, we establish theoretical results regarding signal tensor recovery, identifiability and predictive model consistency as well as the algorithm convergence property.

From a feature extraction perspective, in a sense, the decomposed basis tensor layers in the proposed approach can be viewed as some full-size feature extraction filters analogues to those in the CNN. In contrast to the CNN which applies homogeneous filters over all subjects, the individualized layers in the proposed method serve as subject-specific filters to capture additional individual-wise characteristics and heterogeneity information, which enables us to effectively capture heterogeneous outcome-associated signals such as the tumor-associated microvesicles (TMVs) in breast cancer imaging data. Wang et al. (2017a) propose a similar idea for a multiple graph factorization model via decomposing the edge probability matrix to a population-shared baseline and an individual deviation. However, Wang et al. (2017a)’s model is not applicable to the problem considered in this paper, as their approach focuses on binary data and assumes a symmetric matrix structure, while we have a multimodal three-way tensor covariate with continuous measurements of pixels.

Another advantage of the proposed method is that the individualized layer efficiently integrates multimodality images from the same subject, as it extracts a common structure over different modalities while accounting for modality-specific background noise. Multimodal imaging analysis has drawn great attention in recent years, e.g., in cross-modality registration (Maes et al., 1997; Luan et al., 2008) and multiple image fusion (Wang et al., 2006). To the best of our knowledge, most existing methods model the inter-modality information either after selecting features from each imaging modality separately (Zhang and Shen, 2012; Liu et al., 2014), or by combining all modalities into a joint model but assuming independence between modalities (Hinrichs et al., 2011; Yuan et al., 2012). In contrast, through the individualized layers in tensor decomposition, the proposed model integrates multimodality information in both feature extraction and prediction modeling more effectively.

The paper is organized as follows. Section 2 introduces the main methodology and the background. Section 3 presents the estimation and computation framework. Section 4 establishes theoretical results. Section 5 provides simulation studies. Section 6 illustrates an application for human breast cancer optical imaging data. The last section provides concluding remarks and discussion.

2 Methodology

2.1 Notation and Operations

Tensor, which also refers to a multi-dimensional array, plays a central role in the proposed approach. We start with providing a brief summary of tensor notation. Extensive references can be found in the review reports (Kolda, 2006; Kolda and Bader, 2009). A $D$th-order (or $D$-way) tensor is a $D$-dimensional array $\bm{\mathcal{X}}\in\mathbb{R}^{p_{1}\times p_{2}\times\ldots\times p_{D}}$, where the dimensions of a tensor are also known as modes. We denote $p_{d}$ ($1\leq d\leq D$) as the marginal dimension of the $d$th mode. In this paper, we use the terms way, mode and order interchangeably. The mode-d matricization of $\bm{\mathcal{X}}$, denoted as $\bm{\mathcal{X}}_{(d)}$, is defined as a $(p_{d}\times\prod_{d^{\prime}\neq d}p_{d^{\prime}})$-dimensional matrix, where the $(i_{1},\ldots,i_{D})$th element of $\bm{\mathcal{X}}$ maps to the $(i_{d},j)$th element of $\bm{\mathcal{X}}_{(d)}$ and $j=1+\sum_{d^{\prime}\neq d}(i_{d^{\prime}}-1)\prod_{d^{\prime\prime}<d^{\prime},d^{\prime\prime}\neq d}p_{d^{\prime\prime}}$ (Kolda and Bader, 2009). Moreover, we denote $\mbox{vec}(\cdot)$ as a vectorizing operation which converts a tensor to a vector, where the element $x_{i_{1},\ldots,i_{D}}$ in a $D$-way tensor $\bm{\mathcal{X}}$ is turned to be the $\left(i_{1}+\sum_{d=2}^{D}[(i_{d}-1)\prod_{j=1}^{d-1}p_{j}]\right)$th element in the long vector $\mbox{vec}(\bm{\mathcal{X}})$.

Next we introduce some useful matrix and tensor operations. Given two matrices $\bm{A}\in\mathbb{R}^{m\times n}$ and $\bm{B}\in\mathbb{R}^{p\times q}$, the Kronecker product is an $mp\times nq$ matrix denoted as $\bm{A}\otimes\bm{B}$. If $\bm{A}$ and $\bm{B}$ have the same number of columns $n=q$, then the Khatri-Rao product (Rao and Mitra, 1971) is defined as an $mp\times n$ matrix by $\bm{A}\odot\bm{B}=[\bm{A}_{\cdot 1}\otimes\bm{B}_{\cdot 1}\;\cdots\;\bm{A}_{\cdot n}\otimes\bm{B}_{\cdot n}]$, where $\bm{A}_{\cdot k}$ is the $k$th column of $\bm{A}$. The inner product $\langle\cdot,\cdot\rangle$ of two tensors with the same dimension is defined as $\langle\bm{\mathcal{A}},\bm{\mathcal{B}}\rangle=\langle\mbox{vec}(\bm{\mathcal{A}}),\mbox{vec}(\bm{\mathcal{B}})\rangle$, which follows that $\langle\bm{\mathcal{A}},\bm{\mathcal{A}}\rangle=\parallel\bm{\mathcal{A}}\parallel^{2}_{F}$, where $\parallel\cdot\parallel_{F}$ is the Frobenius norm. Moreover, an outer product “$\circ$” operating on multiple vectors $\bm{b}^{1}\in\mathbb{R}^{p_{1}},\ldots,\bm{b}^{D}\in\mathbb{R}^{p_{D}}$ creates a rank-1 $D$-way tensor $\bm{\mathcal{B}}=\bm{b}^{1}\circ\bm{b}^{2}\circ\ldots\circ\bm{b}^{D}$, where the $(i_{1},i_{2},\ldots,i_{D})$th element of $\bm{\mathcal{B}}$ is defined as $b_{i_{1},\ldots,i_{D}}=b^{1}_{i_{1}}b^{2}_{i_{2}}\ldots b^{D}_{i_{D}}$, and $\bm{b}^{d}=(b^{d}_{1},\ldots,b^{d}_{p_{d}})^{\prime}$ is the factor vector at mode $d$.

Consequently, a $D$-way tensor $\bm{\mathcal{B}}\in\mathbb{R}^{p_{1}\times p_{2}\times\ldots\times p_{D}}$ is of rank $R$ if it can be represented as $\bm{\mathcal{B}}=\sum_{r=1}^{R}\bm{b}_{r}^{1}\circ\bm{b}_{r}^{2}\circ\ldots\circ\bm{b}_{r}^{D},$ where $\bm{b}^{d}_{r}$’s ($r=1,\ldots,R$) are $p_{d}$-dimensional vectors ($d=1,\ldots,D$). For ease of notation, we introduce the Kruskal operator (Kolda and Bader, 2009) for a rank-$R$ $D$-way tensor as $\bm{\mathcal{B}}=\bm{B}^{1}\circ\bm{B}^{2}\circ\cdots\circ\bm{B}^{D}=\llbracket\bm{B}^{1},\bm{B}^{2},\ldots,\bm{B}^{D}\rrbracket$, where $\bm{B}^{d}=[\bm{b}_{1}^{d},\bm{b}_{2}^{d},\ldots,\bm{b}_{R}^{d}]\in\mathbb{R}^{p_{d}\times R}$ denotes the factor matrix at mode $d$. This tensor rank decomposition is also known as CANDECOMP/PARAFAC (CP) decomposition, or Kruskal decomposition (Coppi and Bolasco, 1988). An alternative tensor decomposition is high-order singular value decomposition, also called Tucker decomposition (Tucker, 1966), which decomposes a tensor into a $D$-way core tensor associated with $D$ orthonormal bases-matrices. However, the core tensor in Tucker decomposition is not guaranteed to be diagonal, and thus the tensor rank is not very clear. Therefore we adopt the CP decomposition here.

2.2 Basic Framework and Background

The sample imaging data is usually a mixture of true signal structure and random noises. In general, we assume that

[TABLE]

where $\bm{\mathcal{X}}_{i}\in\mathbb{R}^{p_{1}\times\cdots\times p_{D}}$ is the observed $D$-way sample tensor covariate for the $i$th subject, $\bm{\Theta}_{i}$ is the true signal tensor associated with the outcome response, and $\bm{\mathcal{N}}_{i}$ is the corresponding dimensional noise background. The elements of $\bm{\mathcal{N}}_{i}$ are independent and identically distributed with mean of zero and variance of $\sigma^{2}$.

Let $f(\cdot)$ denote an appropriate feature extraction or transformation mapping on the tensor predictor $\bm{\Theta}_{i}\longmapsto f(\bm{\Theta}_{i}):\mathbb{R}^{p_{1}\times\cdots\times p_{D}}\longrightarrow\mathbb{R}^{q}$. In a sample version with recovered signal $\hat{\bm{\Theta}}_{i}$, let $f(\bm{\mathcal{X}}_{i})=f(\hat{\bm{\Theta}}_{i})$ denote the extracted sample feature. A general supervised learning stage to link the image tensor covariate $\bm{\mathcal{X}}_{i}$ and the observed outcome response $y_{i}$ can be modeled as

[TABLE]

where $\mathcal{L}(\cdot)$ is a selected loss function related to a predictive model based on the response types. For example, an $L_{2}$ loss, $\mathcal{L}(y_{i},f(\bm{\mathcal{X}}_{i}),\bm{\beta})=\|y_{i}-f(\bm{\mathcal{X}}_{i})^{T}\bm{\beta}\|^{2}$, is usually applied for a linear regression model with a continuous response variable, while for the binary response in real data, we can employ a logistic regression where the corresponding $\mathcal{L}(\cdot)$ refers to the negative loglikelihood function. Alternatively, it is also flexible to utilize a machine learning model such as the support vector machine (SVM, Vapnik et al., 1995) with a hinge loss for a classification problem. Furthermore, regularization and penalty can be also incorporated to the model in (2).

Apparently, that feature extraction and data transformation is the key part to utilize the imaging information. One naive transformation method is to directly unfold the imaging tensor $\bm{\mathcal{X}}_{i}$ to a long vector via $f(\bm{\mathcal{X}}_{i})=\mbox{vec}(\bm{\mathcal{X}}_{i})$. However, the number of unknown parameters in (2) using the vectorizing covariates is $\prod_{d=1}^{D}p_{d}$, which is ultra-high-dimensional and would lead to an unestimable model even applying common regularization methods. In addition, some measurement errors on the imaging data are also incorporated to the model in (2). A natural solution to solve these problems is to employ a dimension reduction technique to recover true signal structure and thus extract important features from the tensor covariates, and then fit the model in (2) based on the extracted information. Intuitively, we consider a low-rank approximation for tensor $\bm{\mathcal{X}}_{i}$ as

[TABLE]

where $\bm{\mathcal{B}}_{r}$’s ($r=1,\ldots,R$) are normalized $D$-way tensor bases shared by populations. This rank-$R$ structure can be estimated by minimizing the difference between the observed imaging and approximated values under the Frobenius norm: $\sum_{i=1}^{N}\parallel\bm{\mathcal{X}}_{i}-\sum_{r=1}^{R}w_{ir}\bm{\mathcal{B}}_{r}\parallel^{2}_{F}$, and thus let $f(\bm{\mathcal{X}}_{i})=(\hat{w}_{i1},\ldots,\hat{w}_{iR})^{\prime}$ be the extracted features, where the $\hat{w}_{ir}$’s are the projection loadings of the $\bm{\mathcal{X}}_{i}$ onto the estimated bases. However, without specifying any structures on tensor bases, the general approximation above is still likely to be unestimable and inefficient, as it is essentially a singular value decomposition problem corresponding to an $N\times\prod_{d=1}^{D}p_{d}$-dimensional matrix for the vectorized data.

2.3 Individualized Multilayer Tensor Decomposition

In this subsection, we propose a novel individualized multilayer tensor decomposition framework for true signal recovery and feature extraction. Following the model in (1), we have a $D$-way tensor covariate $\bm{\mathcal{X}}_{i}\in\mathbb{R}^{p_{1}\times\cdots\times p_{D}}$ for each individual ($i=1,\ldots,N$), and let $\bm{\mathcal{X}}=[\bm{\mathcal{X}}_{1}\;\cdots\;\bm{\mathcal{X}}_{N}]\in\mathbb{R}^{N\times p_{1}\times\cdots\times p_{D}}$ denote the aligned higher-order tensor covariate, $\bm{\Theta}=[\bm{\Theta}_{1}\;\cdots\;\bm{\Theta}_{N}]$ denote the aligned true signal tensor, and $\bm{\mathcal{N}}=[\bm{\mathcal{N}}_{1}\;\cdots\;\bm{\mathcal{N}}_{N}]$ denote the corresponding dimensional noise tensor. In general, the observed tensor covariate is mixed with a true signal structure and some random noises: $\bm{\mathcal{X}}=\bm{\Theta}+\bm{\mathcal{N}}.$

To motivate our model, we first start with considering a low-rank structure modeling of the grand signal tensor $\bm{\Theta}$ following a CP decomposition:

[TABLE]

where $\bm{B}^{d}_{\cdot r}$ ($r=1,\ldots,R$) denotes the $r$th column vector of mode-$d$ factor matrix $\bm{B}^{d}\in\mathbb{R}^{p_{d}\times R}$ ($d=1,\ldots,D$), and $\bm{W}=[\bm{w}_{1},\ldots,\bm{w}_{N}]^{\prime}\in\mathbb{R}^{N\times R}$ is the factor matrix corresponding to the dimension of individuals, and $\bm{w}_{i}=(w_{i1},\ldots,w_{iR})^{\prime}$ is the loading vector associated to the $i$th individual. To avoid the indeterminacy from scaling, we let $\bar{\bm{B}}_{\cdot r}^{d}=\frac{1}{\|\bm{B}_{\cdot r}^{d}\|}\bm{B}_{\cdot r}^{d}$ be the normalized factor vector at mode $d$. We further Let $\bm{\mathcal{B}}_{r}=\bar{\bm{B}}_{\cdot r}^{1}\circ\bar{\bm{B}}_{\cdot r}^{2}\circ\ldots\circ\bar{\bm{B}}_{\cdot r}^{D}$ and $\bar{w}_{ir}=w_{ir}\prod_{d=1}^{D}\|\bm{B}_{\cdot r}^{d}\|$, and thus $\|\bm{\mathcal{B}}_{r}\|_{F}=1$ (see supplementary material A.1 for details). Then each individual tensor can be represented as

[TABLE]

where $\bm{\mathcal{B}}_{r}$’s serve as population-wise rank-$1$ basis tensors, and thus the individual-specific feature based on the signal tensor $\bm{\Theta}_{i}$ is pooled as $f(\bm{\Theta}_{i})=\bar{\bm{w}}_{i}=(\bar{w}_{i1},\ldots,\bar{w}_{iR})^{\prime}$. In the sample version, the extracted feature based on the observed sample covariates is $f(\bm{\mathcal{X}}_{i})=f(\hat{\bm{\Theta}}_{i})=\hat{\bar{\bm{w}}}_{i}$, where the estimated signals $\hat{\bm{\Theta}}_{i}$’s as well as the $\hat{\bm{w}}_{i}$’s are obtained by minimizing an $L_{2}$-loss: $\|\bm{\mathcal{X}}-\bm{\Theta}\|^{2}_{F}$ under the model in (3). We refer this model to the higher-order CP decomposition (HOCPD) model. The HOCPD method is powerful for reducing the tensor covariates’ dimensionality, as the feature dimension is effectively reduced from $\prod_{d=1}^{D}p_{d}$ to $R$ for the supervised learning in (2); however, it highly depends on the low-rank assumption on the grand signal tensor, which could fail to capture complex tensor data information if there is significant heterogeneous variation arising from unique individuals. Moreover, it is infeasible to capture heterogeneity simply by increasing the rank in an HOCPD model, as it might lead to non-identifiability and unstable model estimation.

In the following, we propose an individualized multilayer tensor learning (IMTL) method to incorporate heterogeneous structures to tensor decomposition. For the $i$th individual, we assume

[TABLE]

where $\bm{\mathcal{B}}_{r}$’s are the population-shared rank-1 bases analogous to those in (4), and $\bm{\mathcal{S}}_{i}=\bm{s}_{i}^{1}\circ\bm{s}_{i}^{2}\circ\ldots\circ\bm{s}_{i}^{D}$ is an individualized rank-1 structure. With normalized vectors $\bar{\bm{s}}_{i}^{d}$’s and $u_{i}=\prod_{d=1}^{D}\|\bm{s}^{d}_{i}\|$, we can rewrite $\bm{\mathcal{S}}_{i}=u_{i}\bar{\bm{\mathcal{S}}}_{i}$ where $\bar{\bm{\mathcal{S}}}_{i}=\bar{\bm{s}}_{i}^{1}\circ\bar{\bm{s}}_{i}^{2}\circ\ldots\circ\bar{\bm{s}}_{i}^{D}$. Given observed sample tensor covariates, the latent parameters in model (5) are estimated through

[TABLE]

The orthogonal constraint is imposed to ensure the identifiability between the homogeneous layers and the heterogeneous layers, and also to improve the computational stability. Based on the recovered signal $\hat{\bm{\Theta}}_{i}$, the extracted individual feature is $f(\hat{\bm{\Theta}}_{i})=(\hat{\bar{\bm{w}}}^{\prime}_{i},\hat{u}_{i},\{\hat{\bar{\bm{s}}}^{{}^{\prime}d}_{i}\}^{D}_{d=1})^{\prime}$, containing the weights of the population-shared layers $\hat{\bar{\bm{w}}}_{i}=(\bar{w}_{i1},\ldots,\bar{w}_{iR})^{\prime}$, the weight of the individualized layer $\hat{\mu}_{i}$, and the decomposed factors of the individualized layers $\hat{\bm{\bar{s}}}^{d}_{i}$’s, respectively.

In contrast to the HOCPD in (4), each signal tensor $\bm{\Theta}_{i}$ in (5) is characterized by two different layers, one consisting of a linear combination of homogeneous structures ($\bm{\mathcal{B}}_{r}$’s), and another containing an individualized structure ($\bm{\mathcal{S}}_{i}$) capturing the heterogeneity of individual features. Intuitively, the individualized layer $\bm{\mathcal{S}}_{i}$ is the “best” rank-1 structure extracted from the “residual” tensor of a higher-order decomposition. This additional layer enables us to capture the individual-specific information of the tensor covariate which is not characterized by a population structure $\sum\limits_{r=1}^{R}\bar{w}_{ir}\bm{\mathcal{B}}_{r}$, while the rank-1 structure in (5) avoids overfitting on the individualized layers.

Additionally, In practice, the informative heterogeneous signals could be relatively weak compared to the population-wise background and noise. For example, in real data, the important biomarker TMVs are discrete and random spatially distributed, while the signal strength is also quite weak. Therefore, they are very likely to be washed out if only low-rank homogeneous bases are used. In some sense, the homogeneous bases in the population layer serve as some full-size filters to remove those homogeneous variation, so the subject-specific signals which deviate from the low-rank homogeneous background can be captured. Consequently, the population layer and the individual layer effectively coordinate to extract informative features from sample tensor covariates and thus improve model prediction power.

2.4 Multilayer Tensor Reconstruction for Multimodality Data

Multimodality imaging data is widely adopted and has drawn great attentions in recent years. In our motivating application, optical imaging produces multiple images using different wavelengths of light in one examination. Although a single modality model is applicable, there is a critical need to integrate all information collected from multiple modalities, so important features associated with clinical outcomes can be extracted more effectively and the prediction power can be greatly enhanced. In the following, we develop a multimodality imaging tensor model which extends the individualized multilayer approach in (5) to incorporate different sources of modality information.

We consider an M-modality tensor covariate $\bm{\Theta}=[\bm{\Theta}^{(1)}\;\cdots\;\bm{\Theta}^{(M)}]$, where $\bm{\Theta}^{(m)}\in\mathbb{R}^{N\times p_{1}\times\ldots\times p_{D}}$ ($1\leq m\leq M$). For a single modality $m$, we propose a multilayer tensor decomposition

[TABLE]

where $\bm{B}^{(m),d}$’s are modality-specific factors, and $\bm{\mathcal{S}}=[\bm{\mathcal{S}}_{1}\;\bm{\mathcal{S}}_{2}\;\ldots\;\bm{\mathcal{S}}_{N}]\in\mathbb{R}^{N\times p_{1}\times\ldots\times p_{D}}$ contains heterogenous layers. With normalized factors, each individual tensor can be rewritten as

[TABLE]

where $\{\bm{\mathcal{B}}_{r}^{(m)}=\bar{\bm{B}}^{(m),1}_{\cdot r}\circ\cdots\circ\bar{\bm{B}}^{(m),D}_{\cdot r}\}_{r=1}^{R}$ is the set of bases for the $m$th modality, and $\bm{\mathcal{S}}_{i}=\mu_{i}\bar{\bm{\mathcal{S}}}_{i}$ is an individualized layer shared by different modalities $\bm{\Theta}^{(m)}_{i}$ ($m=1,\ldots,M$) of the $i$th subject. Similarly, the proposed model in (6) can be estimated by minimizing a sum of squared loss: $\sum_{m=1}^{M}\parallel\bm{\mathcal{X}}^{(m)}-\bm{\Theta}^{(m)}\parallel^{2}_{F}$ with the constraints. Let $\bm{\Theta}_{i}=[\bm{\Theta}^{(1)}_{i}\ldots\bm{\Theta}^{(M)}_{i}]$ be the multimodality tensor signal for the $i$th subject, and denote $\hat{\bm{\Theta}}_{i}$ as the sample estimator based on (6). The extracted subject-wise feature is $f(\hat{\bm{\Theta}}_{i})=(\hat{\bar{\bm{w}}}_{i},\hat{\mu}_{i},\{\hat{\bar{\bm{s}}}^{d}_{i}\}_{d=1}^{D})^{\prime}$, containing the weights of the modality specific layers $\hat{\bar{\bm{w}}}_{i}=\{\hat{\bar{w}}^{(m)}_{ir}\}^{1\leq m\leq M}_{1\leq r\leq R}$, the weight of the individualized layer $\hat{\mu}_{i}$, and the decomposed factors of the individualized layer $\hat{\bm{\bar{s}}}^{d}_{i}$’s, respectively.

For the proposed multimodality model in (6), the individualized layer $\bm{\mathcal{S}}_{i}$ is a common low-rank structure capturing the individual-specific spatial pattern shared by different modalities. Similarly, for a given subject, if we consider removing the modality-specific layers $\sum\limits_{r=1}^{R}\bar{w}^{(m)}_{ir}\bm{\mathcal{B}}^{(m)}_{r}$ from each modality $\bm{\Theta}_{i}^{(m)}$ first, and then aligning the “residuals” from multiple modalities to a higher-order tensor ($\mathbf{R}^{p_{1}\times\cdots\times p_{D}\times M}$), then the individualized layer can be interpreted as a tensor slice along the modality mode of the “best” rank-1 structure extracted from the higher-order residual tensor.

In addition to characterizing the heterogeneous structures analogues to the singlemodality case, the individualized layer also makes it feasible to integrate the information of multiple modalities effectively. Although the individual-specific signals are randomly distributed across different individuals, their spatial pattern is shared by different modalities from the same individual. Therefore, the individualized layer extracted over multiple modalities aggregates the information regarding their common structure and thus captures this individual-specific feature more accurately. Figure 2 provides an illustration of the individualized layers and the modality-specific layers on the four-modality breast cancer images.

Alternatively, we can simply align the four modalities together to a higher-order tensor predictor with a size of $p_{1}\times...\times p_{D}\times M$, and then apply the proposed singlemodality model in (5). However, this is inadequate to capture the variation among different modalities. In the following, we compare the singlemodality model in (5) with the integrated tensor and multimodality model in (6). Let $\bm{\Theta}_{i}=[\bm{\Theta}^{(1)}_{i},\ldots,\bm{\Theta}^{(M)}_{i}]\in\mathbf{R}^{p_{1}\times\cdots\times p_{D}\times M}$ denote the multimodality integrated tensor. With a singlemodality model, we have

[TABLE]

where $\bm{b}_{r}^{modal}=(b^{modal}_{1r},\ldots,b^{modal}_{Mr})^{\prime}$ and $\bm{s}_{i}^{modal}=(s^{modal}_{i1},\ldots,b^{modal}_{iM})^{\prime}$ are both basis vectors on the mode of modality. Hence, for the $m$th modality ($m=1,2,3,4$), the corresponding tensor slice ($p_{1}\times\cdots\times p_{D}$-dimensional tensor) can be written as

[TABLE]

This implies that the population-shared basis layers $\bm{\mathcal{B}}_{r}=\bm{b}^{1}_{r}\circ\cdots\circ\bm{b}^{D}_{r}$ ($1\leq r\leq R$) are the same across all modalities up to scaling and permutation. In contrast, the multimodality model in (6) allows different modality-specific layers $\bm{\mathcal{B}}^{(m)}_{r}$ across different modalities to capture the modality-wise heterogeneity. In fact, the singlemodality model on the integrated data is a special case of the multimodality model (6), with additional constraints on all modalities sharing the same set of population bases.

3 Computation

In this section, we address the estimation problem of multilayer decomposition in (6) for the proposed IMTL model. In contrast to traditional CP decomposition, incorporating modality layers and individual layers of the proposed method significantly increases the computation cost, and conventional algorithms for tensor decomposition are not necessarily scalable in our situation. Therefore, we propose a bi-level block improvement algorithm which alternately updates different layers and apply a maximum block improvement (MBI) strategy for the estimation of each layer.

The proposed model in (6) yields a constrained optimization requiring orthogonality between the modality layers and the individual layers. Hence, we employ a penalization method (Nocedal and Wright, 2006) to achieve orthogonality constrained optimization, that is, to minimize the objective function

[TABLE]

Note that the orthogonal property of the decomposed layers is not critical here as it only serves to ensure identifiability. Although the penalization method cannot guarantee an exact orthogonality of the estimated layers, it is sufficient for pursuing the layers’ identifiability in practice.

The objective function in (7) is convex and differentiable with respect to each separate block of parameters, which makes it applicable to apply a block-wise updating procedure. However, the traditional tensor decomposition algorithms, e.g., alternating least squares (ALS, Carroll and Chang, 1970) and block relaxation algorithm (De Leeuw, 1994; Lange, 2010), are nearly infeasible and unscalable due to the large number ($(M+N)D+1$) of blocks, and updating one block at each iteration could lead to poor convergence performance. In addition, the estimation of different blocks (layers) utilizes different parts of the data, for example, $\bm{B}^{(m),d}$’s are estimated within each modality and $\bm{s}_{i}^{d}$’s are estimated within each subject. This allows for alternately estimating different layers, and makes parallel computing feasible, which significantly improves the computational scalability. Finally, in general, for $D>2$, the ALS-type cyclical algorithm is not guaranteed to converge to a stationary point (Chen et al., 2012). Therefore, we propose a bi-level algorithm, that is, within the estimation of each layer, we apply the MBI strategy (Chen et al., 2012) in estimating specific blocks, and then we update each layer alternately.

Specifically, given the individual-layer factors, we let $\bm{\mathcal{X}}^{(m)}=\bm{\mathcal{X}}^{(m)}-\hat{\bm{\mathcal{S}}}$ and alternately update blocks of $\bm{W}$ and $\{\bm{B}^{(m),d}\}$’s; that is,

[TABLE]

where $\tilde{\bm{\mathcal{X}}}^{(m)}_{(1)}$ is the mode-1 matricization. The sub-optimization turns to be a ridge-type problem and thus has a unique explicit solution. Next, given $\hat{\bm{W}}$, we update $\{\bm{B}^{(m),d}\}_{d^{\prime}s}$ through parallel computing across different modalities ($m=1,\ldots,M$), and apply the MBI strategy. Let

[TABLE]

be the reduced objective function with respect to $\bm{B}^{(m),d}$ while fixing all other blocks, where $\bm{B}_{-d}^{(m)}=\bm{B}^{(m),D}\odot\cdots\bm{B}^{(m),d+1}\odot\bm{B}^{(m),d-1}\cdots\odot\bm{B}^{(m),1}$ and $\hat{\bm{\mathcal{S}}}_{i,(d)}$ is the mode-d matricization of $\hat{\bm{\mathcal{S}}}_{i}$. Given the latest iterate $\{\bm{B}^{(m),1,[t]},\ldots,\bm{B}^{(m),D,[t]}\}$, we calculate

[TABLE]

and update only the block $\bm{B}^{(m),d,[t]}$ at the mode $d=d_{m}^{*}$ for maximum improvement.

Similarly, given the modality-layer factors, the individualized layers can be estimated within each subject using parallel computing. For each subject $i$, let $\bar{\bm{\mathcal{X}}}^{(m)}_{i}=\bm{\mathcal{X}}^{(m)}_{i}-\sum\limits_{r=1}^{R}\hat{w}_{ir}\hat{\bm{\mathcal{B}}}^{(m)}_{r}$, and let

[TABLE]

be the reduced objective function with respect to $\bm{s}_{i}^{d}$, where $\bm{S}_{i,-d}=\bm{s}_{i}^{1}\odot\cdots\bm{s}_{i}^{d-1}\odot\bm{s}_{i}^{d+1}\cdots\odot\bm{s}_{i}^{D}$ and $\hat{\bm{\mathcal{B}}}^{(m)}_{r,(d)}$ is the mode-d matricization. Given the latest iterate $\{\bm{s}_{i}^{1,[t]},\ldots,\bm{s}_{i}^{D,[t]}\}$, calculate

[TABLE]

and update only mode-$d_{i}^{*}$ block. The detailed algorithm is summarized in following Algorithm 1.

For any objective function $f(\bm{A}_{1},\ldots,\bm{A}_{n})$ with blocks of parameters $\{\bm{A}_{i}\}_{i=1}^{n}$, a point of $(\bm{A}^{*}_{1},\ldots,\bm{A}^{*}_{n})$ in the parameter space is defined as a block-wise stationary point of $f(\cdot)$ if for any block, there is $\bm{A}^{*}_{i}=\arg\min_{\bm{A}_{i}}f(\bm{A}^{*}_{1},\ldots,\bm{A}^{*}_{i-1},\bm{A}_{i},\bm{A}^{*}_{i+1},\ldots,\bm{A}^{*}_{n}).$

Let $\{\bm{W},\bm{B}^{(m),d},\bm{s}^{d}_{i}\}_{m,i,d}$ be an element-wise collection of all parameters and let $\Omega$ be the corresponding parameter space. We provide a global convergence property of Algorithm 1 as follows.

Lemma 1**.**

Assume $\Omega$ is compact, then any accumulative point of the iterations from Algorithm 1, say $\{\bm{W}^{*},\bm{B}^{*(m),d},\bm{s}^{*d}_{i}\}_{m^{\prime}s,i^{\prime}s,d^{\prime}s}$, is a block-wise stationary point of the objective function in (7).

In general, multiple initializations are suggested to obtain a sound optimum in (7). The normalization on modes’ factor vectors can be performed after the last iteration. Given the estimated tensor signals $\hat{\bm{\Theta}}_{i}$’s, for binary response, we employ an $L_{1}$-penalized logistic regression model with the extracted sample features $f(\hat{\bm{\Theta}}_{i})$’s in (2) for prediction. The rank of the population-shared layers and the other tuning parameter are selected based on a grid search to minimize the prediction errors on the validation set or via a cross-validation method. More details and discussion about tuning procedure and the prediction on a new subject are provided in the supplementary material.

4 Theory

In this section, we develop the theoretical framework for the proposed model regarding both tensor signal recovery and prediction modeling. First, we introduce conditions to ensure identifiability of the proposed multilayer tensor modeling. Furthermore, we provide an accurate error bound for the recovered signal tensor and show that the estimated signal tensor converges to the true one. Finally, we present the theoretical results of the supervised learning stage to evaluate prediction performance. All proofs are provided in supplementary materials Section A.

We first address the tensor modeling identifiability issue before establishing the statistical property, which is critical for tensor representation. In particular, we provide sufficient conditions to achieve identifiable layers in the proposed multilayer tensor model, which can be verified easily in practice. For ease of notation, the following discussion focuses on the single modality model while the presented conditions can be easily extended to the multimodality model.

In the proposed framework, potential unidentifiability could occur in the multi-layer CP decomposition of tensor predictors in (5), which can be attributed to three aspects. The first two indeterminacies arise from scaling and permutation, and the last aspect is the non-uniqueness of the CP decomposition for a tensor, that is, the possibility that more than one combination of population layers and individualized layers can lead to the underlying true image tensor. In the proposed model, the scaling indeterminacy, which refers to possible rescaling over different modes’ factor vectors for each layer, is eliminated by imposing a unit-norm constraint on parameterization, that is, $\|\bar{\bm{B}}^{d}_{\cdot r}\|_{2}=1$ and $\|\bar{\bm{s}}^{d}_{i}\|_{2}=1$. Moreover, the permutation indeterminacy refers to the arbitrary reordering of the population bases. To address this point, we could align the population bases according to a descending order of the first element of mode-1 factor vectors, that is, $\bar{\bm{B}}^{1}_{11}\geq\bar{\bm{B}}^{1}_{12}\geq\cdots\geq\bar{\bm{B}}^{1}_{1R}.$

After controlling the scaling and the permutation, in general, the CP decomposition of a tensor could still be non-unique, due to the possibility of multiple combinations of rank-one tensors in decomposition. Although various identifiability conditions have been presented for a conventional CP decomposition (Coppi and Bolasco, 1988; Sidiropoulos and Bro, 2000; Kolda and Bader, 2009), they are infeasible here as all of them require checking each individual tensor $\bm{\Theta}_{i}$ ($1\leq i\leq N$) separately, which is not effective in practice, especially when the sample size $N$ is increasing.

In the following, we provide a much weaker sufficient condition based on the integrated higher-order tensor without imposing any additional constraints on parameterization. We first introduce the concept of a $k$-rank of a matrix according to (Coppi and Bolasco, 1988). Specifically, the $k$-rank of a matrix $\bm{A}$, denoted as $\mathcal{K}_{\bm{A}}$ is defined as \mathcal{K}_{\bm{A}}=\max\{k:\;\text{any k$columns of$\bm{A} are linearly independent}\}. Let $\bm{\Theta}_{[1:n]}$ denote an integrated $(D+1)$-way tensor combining $n$ individual tensors. Without loss of generality, we assume $\bm{\Theta}_{[1:n]}=[\bm{\Theta}_{1}\;\cdots\;\bm{\Theta}_{n}]$. There is a $(R+n)$-rank representation for the integrated tensor:

[TABLE]

where $\bm{w}^{[1:n]}_{r}=(w_{1r},\ldots,w_{nr})^{\prime}$ and $\tilde{\bm{u}}^{[1:n]}_{i}=(\underbrace{0,\ldots,0}_{1,\ldots,i-1},\underbrace{1}_{i},\underbrace{0,\ldots,0}_{i+1,\ldots,n})^{\prime}$. For $1\leq d\leq D$, let $\tilde{\bm{B}}^{d}_{[1:n]}=[\bm{b}^{d}_{1}\;\cdots\;\bm{b}^{d}_{R}\;\bm{s}^{d}_{1}\;\cdots\;\bm{s}^{d}_{n}]$ denote the mode-d factor matrix for integrated tensor $\bm{\Theta}_{[1:n]}$. We have the following proposition providing a sufficient condition for the identifiability of the multi-layer decomposition in (6).

Proposition 1**.**

If there exist $n$ individual tensors ($2\leq n\leq N$) in $\{\bm{\Theta}_{i},1\leq i\leq N\}$, such that,

[TABLE]

holds for the integrated high-order tenor $\bm{\Theta}_{[1:n]}$, then the multi-layer decomposition in (5) is unique given the unit-norm and the ordering constraints on the factor vectors.

In the proof of Proposition 1, we show that any rank-1 CP decomposition for a $D$-way ($D\geq 2$) tensor is unique up to scaling indeterminacy, which implies the identifiability of the individualized layer given the population bases. The above condition is relatively weak compared to those in Coppi and Bolasco (1988), Sidiropoulos and Bro (2000), and Kolda and Bader (2009), and it is easy to satisfy as it applies for an arbitrary $n$. For example, if $n=2$ and the factor matrices are of full rank, then the condition in Proposition 1 holds as long as $D\geq 2$.

Next we establish the theoretical properties for the recovered tensor signal $\hat{\bm{\Theta}}$ based on the observed sample tensor covariates $\bm{\mathcal{X}}^{(m)}_{i}$’s. We denote $\bm{\gamma}=\left(\mbox{Vec}(\bm{W})^{\prime},\{\mbox{Vec}(\bm{B}^{(m),d})\}^{\prime}_{m,d},\{\bm{s}_{i}^{d}\}^{\prime}_{i,d}\right)^{\prime}$ as the vector of all latent variable parameters. It is straightforward that $\dim(\bm{\gamma})=MR(\sum_{d=1}^{D}p_{d})+N(R+\sum_{d=1}^{D}p_{d})$. Let $\mathcal{I}=\{\bm{\mathcal{X}}_{i}^{(m)}\}_{i,m}$ denote a collection of all observed $D$-way tensors and $|\mathcal{I}|$ be the cardinality measure of $\mathcal{I}$. That is, $|\mathcal{I}|$ denotes the number of all singlemodality images, which equals $NM$ if there is no missing.

Furthermore, let $\bm{\Theta}^{(m)}=(\theta^{(m)}_{i,j_{1}\cdots j_{D}})$ and we have $\theta^{(m)}_{i,j_{1}\cdots j_{D}}=\sum\limits_{r=1}^{R}w_{ir}b_{j_{1}r}^{(m),1}\cdots b_{j_{D}r}^{(m),D}+s_{i,j_{1}}^{1}\cdots s_{i,j_{D}}^{D}$ according to the proposed multilayer model in (6). We assume that $\mathbf{E}[X^{(m)}_{i,j_{1}\ldots j_{D}}]=\theta^{(m)}_{i,j_{1}\ldots j_{D}}$, where $X^{(m)}_{i,j_{1}\ldots j_{D}}$ denotes an element of the sample image tensor $\bm{\mathcal{X}}_{i}^{(m)}$, for example, a pixel in the image, and $x^{(m)}_{i,j_{1}\ldots j_{D}}$ denotes an observed value. Then, the number of all sample elements’ observations is $|\mathcal{I}|\prod_{d=1}^{D}p_{d}$, which increases as the number of images increases. In the following, for the $(j_{1},\ldots,j_{D})$th element of image $\bm{\mathcal{X}}_{i}^{(m)}$, we define the loss function

[TABLE]

Consequently, we assume that the overall objective function is an additive form of the loss function and the penalty function, that is,

[TABLE]

where $\lambda_{|\mathcal{I}|}$ is a penalization coefficient. Suppose that $\mathcal{S}$ is the parameter space of $\mathbf{\Theta}$, and that

[TABLE]

In practice, each pixel of a tensor image can only range from white to black and is usually normalized. Hence, it is sensible to assume that $\|\bm{\Theta}\|_{\infty}\leq C_{0}$, $\|\bm{\gamma}\|_{\infty}\leq C_{1}$ and $\|\bm{\mathcal{X}}\|_{\infty}\leq C_{2}$ for large constants $C_{0}\geq 0$, $C_{1}\geq 0$ and $C_{2}\geq 0$. Then we define the vector parameter space $\mathcal{S}_{\bm{\Theta}}=\{\bm{\Theta}:\|\bm{\Theta}\|_{\infty}\leq C_{0}\}$ and $\mathcal{S}_{\bm{\gamma}}=\{\bm{\gamma}:\|\bm{\gamma}\|_{\infty}\leq C_{1}\}.$

For each $X^{(m)}_{i,j_{1}\cdots j_{D}}$, let $l_{\Delta}(\mathbf{\Theta}|X^{(m)}_{i,j_{1}\cdots j_{D}})=l(\mathbf{\Theta},X^{(m)}_{i,j_{1}\cdots j_{D}})-l(\mathbf{\Theta}_{0},X^{(m)}_{i,j_{1}\cdots j_{D}})$ be the loss difference, where $\mathbf{\Theta}_{0}$ corresponds to the unique true parameter. We first define:

[TABLE]

which is the expected loss difference. Since $\mathbf{\Theta}_{0}$ is the unique true parameter, we have $K(\mathbf{\Theta},\mathbf{\Theta}_{0})\geq 0$ for all $\mathbf{\Theta}\in\mathcal{S}_{\bm{\Theta}}$ and $K=0$ if and only if $\mathbf{\Theta}=\mathbf{\Theta}_{0}$. Therefore, we define the distance between $\mathbf{\Theta}$ and $\mathbf{\Theta}_{0}$ as $\rho(\mathbf{\Theta},\mathbf{\Theta}_{0})=K^{1/2}(\mathbf{\Theta},\mathbf{\Theta}_{0})$, and also define the variance of the loss difference as follows:

[TABLE]

Under the $L_{2}$-loss, it is expected that $K(\mathbf{\Theta},\mathbf{\Theta}_{0})=\frac{1}{NMp_{1}\cdots p_{D}}\|\mathbf{\Theta}-\mathbf{\Theta}_{0}\|_{F}^{2}$, and that $V(\mathbf{\Theta},\mathbf{\Theta}_{0})=\frac{4\sigma^{2}}{NMp_{1}\cdots p_{D}}\|\mathbf{\Theta}-\mathbf{\Theta}_{0}\|_{F}^{2}$, where $\sigma^{2}$ is assumed to be the same variance of each element of the tensor.

Theorem 1**.**

Suppose $\hat{\mathbf{\Theta}}$ is the sample estimator satisfying (11), then we have

[TABLE]

for $\tau_{|\mathcal{I}|}=\max(\varepsilon_{|\mathcal{I}|},\lambda^{1/2}_{|\mathcal{I}|})$, where $c_{1}>0$ is a constant and $\varepsilon_{|\mathcal{I}|}\sim{|\mathcal{I}|}^{-1/2}$. The best possible rate is achieved at $\tau_{|\mathcal{I}|}=|\mathcal{I}|^{-1/2}$ when $\lambda_{|\mathcal{I}|}\sim\varepsilon^{2}_{|\mathcal{I}|}$.

Theorem 1 indicates that, with an appropriate rate of the penalty term shrinking to zero, the recovered tensor signal in (11) converges to the true one as the number of images goes to infinity. In other words, assuming the identifiability conditions are satisfied, Theorem 1 ensures that the extracted information from the sample images captures the true underlying information as the sample size increases. In addition, the convergence rate is established under the $L_{2}$ distance, which is a special case of the Kullback-Leibler divergence.

Remark 1: In practice, the estimator specified in (11) is a global minimizer and may not be attainable due to the non-convex nature of the proposed optimization problem. One possible solution is to follow (1.1) of Shen (1998) and relax (11) to allow an approximate global minimizer, in the sense that the IMTL estimator approaches the global minimizer asymptotically, which is less restrictive on the requirement of the global optimum. Another solution is to develop theoretical properties directly for the IMTL estimator, which is viable because of the block-wise convexity of the proposed criterion function (Li and Zhang, 2017; Wang et al., 2017b). However, this may require good initial values in order to establish the consistency property.

Next, we investigate the theoretical property of the predictive model with the extracted sample information. In the proposed framework, we consider a generalized linear model (McCullagh and Nelder, 2018) assuming that the response variable $Y_{i}$ is associated to the true image signal ($\bm{\Theta}_{i}$) via $\mathbf{E}(Y_{i}|\bm{\Theta_{i}})=\mu(\eta_{i})\quad\text{and}\quad\eta_{i}=f(\bm{\Theta_{i}})^{T}\bm{\beta},$ where $f:\mathbb{R}^{p_{1}\times\ldots\times p_{D}\times M}\rightarrow\mathbb{R}^{p_{f}}$ is the feature-extraction mapping following (6), and $\mu(\cdot)$ is a link function. Let $F(\bm{\Theta}_{|\mathcal{I}|})=(f(\bm{\Theta}_{1})^{T},\ldots,f(\bm{\Theta}_{N})^{T})^{T}$ denote the extracted feature covariates based on the true imaging signals, and let $\mathcal{L}\left(y_{i},f(\bm{\Theta}_{i}),\bm{\beta}\right)=\mathcal{L}(y_{i},\eta_{i}\{f(\bm{\Theta}_{i}),\bm{\beta}\})$ denote a general loss function analogues to that in (2). Consequently, a general criterion function for supervised learning with true imaging information can be presented as

[TABLE]

where $p^{*}(\bm{\beta})$ is a penalty term imposed on feature covariates effect $\bm{\beta}$, $\lambda_{\bm{\beta}}$ is the associated tuning parameter, and $\bm{y}=(y_{1},\ldots,y_{N})^{\prime}$ is the sample response.

Note that the model in (12) essentially refers to a “true” model with the underlying true imaging information. Since only sample images with noises are observed, we actually use the recovered information $\hat{\bm{\Theta}}_{|\mathcal{I}|}$ and the corresponding sample features $F(\bm{\mathcal{X}}_{|\mathcal{I}|})=F(\hat{\bm{\Theta}}_{|\mathcal{I}|})$ in the supervised learning model (12), leading to a sample objective function. Define $\Omega_{\bm{\beta}}$ and $\Omega_{\eta}$ as the parameter spaces for $\bm{\beta}$ and $\eta_{i}$ ($1\leq i\leq N$), respectively. In the following, we introduce a regularity condition on the loss function $\mathcal{L}(\cdot)$:

(C1): For $\eta_{i},\tilde{\eta}_{i}\in\Omega_{\eta}$ and all $y_{i}$’s, there exists a function $K_{i}(\cdot)$ such that

[TABLE]

and $\frac{1}{N}\sum_{1}^{N}\mathbf{E}[K_{i}(y_{i})^{2}]\leq K_{0}$ holds for a positive constant $K_{0}$.

Lemma 2**.**

Suppose condition (C1) holds. Under regularity condition (R1) provided in supplementary material (A.5), and if $\frac{1}{N^{1/2}}\|F(\bm{\mathcal{X}}_{|\mathcal{I}|})-F(\bm{\Theta}_{|\mathcal{I}|})\|_{F}\rightarrow 0$, then, for any given $\bm{\beta}\in\Omega_{\bm{\beta}}$, we have

[TABLE]

Lemma 2 indicates that the sample prediction model approaches the true model asymptotically as long as the recovered imaging signal converges to the true one. Next, we investigate the asymptotic property of the sample model estimation. Let $\bm{\beta}_{0}$ be the true value of $\bm{\beta}$ denoting the effect of the true imaging signal, let $\hat{\bm{\beta}}(\bm{\mathcal{X}}_{|\mathcal{I}|})=\arg\min\limits_{\bm{\beta}}G_{N}(\bm{\beta}|\bm{y},F(\bm{\mathcal{X}}_{|\mathcal{I}|}))$ be the sample estimator with observed image covariates, and let $\tilde{\bm{\beta}}(\bm{\Theta}_{|\mathcal{I}|}){=}\arg\min\limits_{\bm{\beta}}G_{N}(\bm{\beta}|\bm{y},F(\bm{\Theta}_{|\mathcal{I}|}))$ be the oracle estimator with true imaging signals. Furthermore, we denote $\mathcal{R}(\hat{\mathcal{L}}_{N})=\frac{1}{N}\sum_{i=1}^{N}\mathcal{L}(y_{i},f(\bm{\mathcal{X}}_{i}),\hat{\bm{\beta}}(\bm{\mathcal{X}}_{|\mathcal{I}|}))$ as the empirical risk for the sample loss function, and denote $\mathcal{R}(\tilde{\mathcal{L}}_{N})=\frac{1}{N}\sum_{i=1}^{N}\mathcal{L}(y_{i},f(\bm{\Theta}_{i}),\tilde{\bm{\beta}}(\bm{\Theta}_{|\mathcal{I}|}))$ as the empirical risk for the oracle loss function with true imaging information. The next results provide the estimation consistency of the imaging covariates’ effect as well as the empirical error.

Theorem 2**.**

Suppose the conditions in Theorem 1 and the identifiability conditions in Proposition 1 hold. Under regularity conditions (R1-R4) (supplementary material A.5) and condition (C1), if $\lambda_{\bm{\beta}}\rightarrow 0$, as $|\mathcal{I}|\rightarrow\infty$, we have

[TABLE]

and also,

[TABLE]

In addition to model estimation consistency, Theorem 2 indicates that the empirical error of the sample predictive model approaches the oracle one, where the true imaging signals are known. With an appropriately selected predictive model in (2), the proposed approach is able to achieve a consistent and efficient prediction result.

Remark 2: The exact generalization error bound of the proposed model actually depends on a particular predictive model specified in (12) as well as the response types. Consider a new observation with true signal $(Y^{*}_{j},\bm{\Theta}^{*}_{j})$, in a generalized linear model, Theorem 2 indicates that we would have a consistent mean response prediction $\hat{\mu}^{*}_{j}=\mu(f(\bm{\Theta}^{*}_{j})^{T}\hat{\bm{\beta}}(\bm{\mathcal{X}}_{|\mathcal{I}|}))\rightarrow_{p}\mathbf{E}[Y_{j}^{*}|\bm{\Theta}_{j}^{*}]$, where $\mu(\cdot)$ is a link function. Furthermore, for a continuous and bounded response, the $L_{2}$-prediction-risk $\mathbf{E}(\hat{Y}_{j}^{*}{-}Y_{j}^{*})^{2}$ converges to the theoretical lower bound $\mbox{Var}_{\bm{\Theta}_{j}^{*}}(Y^{*}_{j})$ with $\hat{Y}^{*}_{j}=\hat{\mu}^{*}_{j}$; while for a binary response $Y^{*}_{j}$, let $\hat{\pi}(\bm{\Theta}_{j}^{*})=\mathbf{1}[\hat{\eta}(\bm{\Theta}_{j}^{*})>1/2]$ be the trained classification rule, where $\hat{\eta}(\bm{\Theta}_{j}^{*})=\mbox{logit}(f(\bm{\Theta}^{*}_{j})^{T}\hat{\bm{\beta}})$, then we have $\hat{\pi}(\bm{\Theta}_{j}^{*})\rightarrow_{p}\pi^{0}_{j}$, the optimal Bayes’ rule.

5 Simulation Studies

In this section, we investigate the performance of the proposed approach comparing with other competing methods through a simulation study of heterogeneous weak signals, which mimics the real data structure and also frequently arise in clinical diagnosis and cancer imaging analysis. Specifically, we simulate four-modality imaging data, where multiple modalities imaged from the same individual share random-location heterogeneous signals, while each modality contains its unique background bases. Additional simulations regarding various heterogeneous signal patterns with the singlemodality data and the impact of the contrast between population-shared components and individualized components are provided in the supplementary materials Section B.

In this study, we simulate the $m$th-modal image for the $i$th subject as $\bm{\mathcal{X}}^{(m)}_{i}=\bm{\mathcal{A}}^{(m)}_{i}+\bm{\mathcal{B}}_{i}+\bm{\mathcal{N}}_{i}$, $m=1,\ldots,4$, consisting of the true signal feature image $\bm{\mathcal{B}}_{i}$, the modality-specific background image $\bm{\mathcal{A}}^{(m)}_{i}$ and the noise background image $\bm{\mathcal{N}}_{i}$. For each subject, we randomly select $s_{i}$ pixels of $\bm{\mathcal{B}}_{i}$ to be valued of $2$ (signal pixels) with the other pixels of 0. We generate the response label $y_{i}$ from a Bernoulli distribution with a probability of $0.4$ to be $1$. The number of the signal pixels $s_{i}$ is generated from a Poisson distribution with means $\mu_{C}=25$ and $\mu_{N}=5$ for the cancer subject’s image (given $y_{i}=1$) and the normal subject’s image (given $y_{i}=0$), respectively. The noise elements of $\bm{\mathcal{N}}_{i}$ are generated from $N(0,0.2^{2})$.

Moreover, the first modality background $\bm{\mathcal{A}}^{(1)}_{i}$ is a random noise matrix with elements generated from $N(0,0.1^{2})$; the second modality has a uniform background with $\bm{\mathcal{A}}^{(2)}_{i}=w^{(2)}_{i}\mathbf{1}_{D}\mathbf{1}_{D}^{T}$, where $\mathbf{1}_{D}$ is a $D\times 1$ vector of 1’s and $w^{(2)}_{i}$ is generated from $N(0,0.1^{2})$; and both the third and fourth modality imaging have low-rank structures with $\bm{\mathcal{A}}^{(m)}_{i}=\sum_{r=1}^{5}w^{(m)}_{ir}\bm{a}^{(m),1}_{r}\circ\bm{a}^{(m),2}_{r}$ ($m=3,4$), where $w^{(3)}_{ir}$ and $w^{(4)}_{ir}$ are generated from $N(0,0.1^{2})$, $\bm{a}^{(m),1}_{r}$’s and $\bm{a}^{(m),2}_{r}$’s are generated from $N(\bm{0},\mathbf{I}_{D})$ and $N(\bm{0},0.5\mathbf{I}_{D})$, respectively, and $\mathbf{I}_{D}$ is the $D$-dimensional identity matrix. We set the sample size as $60$ and $100$, equally for the training set, the validation set and the testing set, and the marginal imaging dimension as $D=64$. Figure 5 illustrates sample normal and cancerous images.

To evaluate the prediction performance, we calculate the prediction accuracy rate ($\mathbf{P}[\hat{y}=y]$), the sensitivity ($\mathbf{P}[\hat{y}=1|y=1]$) and the specificity ($\mathbf{P}[\hat{y}=0|y=0]$) on the testing test. We compare the proposed IMTL model with the higher-order CP decomposition method (HOCPD) in (4), the marginal principal component analysis (MPCA) method by mimicking Caffo et al. (2010), the vectorizing $L_{1}$-penalized logistic regression model (VPL), the tensor regression model (TR, Zhou et al., 2013) and the convolutional neural network (CNN).

The VPL is applied on a $16,384$-dimensional vector predictor by vectorizing all four modalities and then fits an $L_{1}$-penalized logistic regression model for the binary response, which is implemented by the R package “glmnet” (Friedman et al., 2010). The MPCA first extracts feature from each individual modality separately and then fits a logistic model with all modal features. The TR and the HOCPD are applied on the integrated multimodality image (third-order tensor predictor). The TR method is implemented by Zhou (2013)’s Matlab toolbox “TensorReg”. In the following numerical studies and real data analysis, we employ an $L_{1}$-penalized logistic regression model for the IMTL, the HOCPD and the MPCA methods at the predictive stage using the extracted features. The relevant tuning parameters associated with each model are selected through minimizing the prediction error rates ($\mathbf{P}[\hat{y}\neq y]$) on the validation set, respectively.

In particular, the CNN method is implemented by Matlab toolbox “matcovnet” (Vedaldi and Lenc, 2015), which inputs four modalities as four channels. Since the tuning of the CNN is crucial, we tune the CNN to minimize the classification error rates on the validation set over the number of layers, the number of filters on each convolution layer, the filter size, depth and the stride size, the pooling window size, the pooling stride size, the pooling method, the activation function, the number of epochs and the learning rates. The detailed CNN architectures and the tuning process in our numerical studies are provided in the supplementary material Section C. In addition, we also implement the TensorFlow for the CNN and the results are consistent with the Matlab’s results.

Table 1 provides the prediction results based on 100 replications. The proposed method (IMTL) outperforms other methods with the highest prediction accuracy ($0.86$ and $0.94$) and sensitivity ($0.69$ and $0.88$) for both sample sizes ($60$ and $100$), respectively. Moreover, the IMTL, the HOPCD methods and the CNN, have significant advantages over the VPL and the MPCA methods which assume the independence between the four modalities. This indicates that integrating different modalities’ information at feature extraction enhances the prediction power. Furthermore, the proposed IMTL method achieves more than $16\%$ improvement in prediction accuracy than the HOCPD, indicating that the proposed method is more effective in utilizing correlation information among different modalities with additional individualized layers. An increased training data size improves the CNN’s prediction power, however, the IMTL still outperforms the CNN with a $10\%$ higher overall accuracy and a $13.4\%$ higher sensitivity at a sample size of 100.

Although the CNN is a powerful tool for image classification, it requires a large number of training samples to entails multiple hidden layers and involves a large number of parameters (Keshari et al., 2018; Wagner et al., 2013; Abbasi-Asl and Yu, 2017a, b). In practice, the sufficient sample size is not always attained, especially for cancer imaging data. In addition, the heterogeneity is another great challenge occurred in cancer imaging, while the CNN is less capable of capturing heterogeneous signals with relatively weak strength.

We conduct additional simulations to investigate the performance of the CNN compared with the proposed IMTL in various situations. First, we compare the performance of the CNN and the IMTL, with different training sample sizes. Figure 3 shows that the IMTL consistently outperforms the CNN, especially when the training sample size is limited, while the CNN gradually “catches up” as the sample size increases. This confirms that the CNN requires a large size of training samples in order to perform well. The image size in this simulation is only $64\times 64$. According to the trend in Figure 3, for real data with an image size of hundreds to thousands, the CNN could require quite a substantial number of training samples to guarantee a good performance.

Second, we investigate the performance of the CNN and the proposed IMTL with varying heterogeneous-signal strengths. To mimic real data, we utilize the concentration intensity of the TMVs to measure the strength and the heterogeneity of signals, which is controlled by the mean value of the number of the TMVs ($\mu_{C}$) in cancer images, while fixing the normal images’ mean ($\mu_{N}=5$). Therefore, the larger the $\mu_{C}$ is, the stronger the signal contrast between the cancerous images and the normal images is. On the other hand, if $\mu_{C}$ is larger, then the heterogeneity is less since the TMVs tend to fill the entire cancer image, and thus the heterogeneity of the randomly-located signals is reduced. We set the training sample size as 60. Figure 4 shows the prediction power of the CNN and the IMTL with a varying $\mu_{C}$. It is clear that the IMTL outperforms the CNN consistently, while the CNN improves as the signal strength becomes stronger and the heterogeneity reduces. This simulation confirms that the proposed IMTL is more powerful than the CNN in situations where the heterogeneity is high and the signal strength is weak.

6 Real Data: Multiphoton Imaging Data for Breast Cancer

We apply the proposed method to the multimodality optical imaging data for breast cancer produced by Boppart Lab (Tu et al., 2016) at University of Illinois Urbana-Champaign. To better visualize the tissue’s biological structure at cellular and molecular levels, Tu et al. (2016)’s multiphoton microscope generates multimodal images through emission of different numbers of photons. Specifically, there are two-photon auto-fluorescence (2PAF), three-photon auto-fluorescence (3PAF), second-harmonic generation (SHG) and third-harmonic generation (THG). Two-photon-fluorescence microscopy is commonly used to visualize tissue morphology and physiology at a cellular level, and three-photon-fluorescence with longer wavelength can reach deeper levels of the tissue and thus provide higher imaging resolution (Chu et al., 2005; Horton et al., 2013). This new technique is able to capture the important tumor-associated microvesicles (TMVs) which are not easily identified by conventional imaging tools such as histology imaging.

Figure 6 illustrates the four modalities for a human subject’s normal breast tissue and a human subject’s cancerous breast tissue. In contrast to the normal tissue, multiple modalities clearly indicate a large number of TMVs which spread out in the microenvironment on the cancerous imaging, particularly in the 2PAF image and the 3PAF image. The normal imaging also presents microvesicles in the 3PAF image, but they are more sporadic with much less intensity.

Prior knowledge in cancer detection shows that informative TMVs are frequently observed in the microenvironment between certain biological organizations such as at the lipid boundary area and around the stroma (Tu et al., 2016). Therefore, we study segmented imaging from three normal human individuals and two cancerous individuals. Specifically, we generate sample images through segmentation, which leads to $107$ normal subjects and $53$ cancerous subjects, each subject having four modalities of images, with each modality of $100\times 100$ pixels. To fit the predictive models and evaluate the diagnosis performance, we randomly split the total sample into a training set, a validation set and a testing set with $60$, $40$ and $60$ subjects, respectively. The prediction results are evaluated by prediction accuracy, sensitivity and specificity on the testing set based on $30$ replications.

We compare the proposed IMTL method to the five methods described in Section 5. Table 2 provides the averaged prediction results on the testing set, which indicates that the proposed IMTL method outperforms the other methods significantly in terms of achieving the highest overall prediction accuracy and sensitivity. Specifically, the prediction accuracy of the proposed method improves by 34%, 39%, 19%, 13% and 7% compared to the VPL, TR, MPCA, HOCPD and CNN approaches, respectively. The boxplot in Figure 7 suggests that the proposed method also has the smallest standard error in overall prediction accuracy. Note that the outstanding predictive power of the proposed method is mainly due to its highest sensitivity score, which is also of primary interest for cancer diagnosis. In practice, the proportion of the potential cancer patients is rather small compared to the general population, and correctly detecting cancer at earlier stages of cancer development is crucial and critical.

For this data, the VPL method and the tensor regression model perform inadequately compared to other methods due to the random location and the weak signal strength of the TMVs. Moreover, without incorporating individual-wise heterogeneity, the MPCA and the HOCPD are inefficient with low sensitivity. In general, the CNN provides an acceptable prediction result, however, due to limited sample size and heterogeneous imaging features, the CNN is not as effective in identifying cancerous subjects with only a 73% sensitivity rate compared to the IMTL’s 90% sensitivity rate. In addition, Figure 8 provides an illustration of the prediction performance using a single modality only compared to utilizing all modalities. Figure 8 clearly shows that integrating multiple modalities improves the predictive power, especially in detecting cancerous subjects.

7 Discussion

In this article, we propose an individualized multilayer tensor learning model incorporating multimodality imaging tensor covariates to predict targeted responses. In the proposed model, we extract both individual-specific and population-wise important features simultaneously from higher-order tensor covariates through different layers, and then fit a prediction model with the extracted features. We illustrate the proposed method through numerical studies and human breast cancer imaging application on both singlemodality and multimodality data.

A major contribution of the proposed method is that we achieve heterogeneous tensor decomposition through utilizing an individualized layer in addition to population-shared modality-specific structure. Our method is motivated by a multimodality imaging study for breast cancer diagnosis, where the biomarker TMVs’ are distributed with great heterogeneity and each modality has its unique background features. Most existing methods assuming homogeneous structure on signals’ features are either infeasible or inefficient in our situation. In contrast, the proposed different layers are capable of capturing individual-specific spatial features through integrating different modalities’ imaging information for the same individual. Both numerical studies and theoretical results demonstrate that the proposed method can achieve higher diagnostic accuracy.

In the proposed method, multilayer tensor decomposition in the first stage is not connected to the response variable directly, which may not guarantee an optimal feature extraction for supervised learning in the second stage. For future research, it is worth developing a supervised feature extraction scheme which can be more powerful in predicting outcomes. We point out three potential directions. The first one is to optimize both the classification loss and the image tensor reconstruction loss, where the latter can be treated as a part of regularization. The second possible direction is to incorporate sufficient dimension reduction techniques to search the “best” bases layers conditional on the response variable. The third direction is to link the feature extraction stage and the prediction stage in a hierarchical form mimicking the CNN framework, and update the first stage of feature extraction adaptively based on the prediction loss.

**Supplementary Materials

**

The online supplement contains all technical proofs, additional numerical results and computation details.

**Acknowledgments

**

The authors would like to acknowledge support for this project from the National Science Foundation grants DMS-1415308, DMS-1613190, DMS-1821198. The authors are grateful to reviewers, the Associate Editor and Editor for their insightful comments and suggestions which have improved the manuscript significantly.