{'fraction_non_alphanumeric': 0.03920004668988191, 'fraction_numerical': 0.016419275139583293, 'mean_word_length': 4.563805606667389, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coupled Support Tensor Machine Classification for Multimodal Neuroimaging Data

Multimodal data arise in various applications where information about the same phenomenon is acquired from multiple sensors and across different imaging modalities. Learning from multimodal data is of great interest in machine learning and statistics research as this offers the possibility of capturing complementary information among modalities. Multimodal modeling helps to explain the interdependence between heterogeneous data sources, discovers new insights that may not be available from a single modality, and improves decision-making. Recently, coupled matrix-tensor factorization has been introduced for multimodal data fusion to jointly estimate latent factors and identify complex interdependence among the latent factors. However, most of the prior work on coupled matrix-tensor factors focuses on unsupervised learning and there is little work on supervised learning using the jointly estimated latent factors. This paper considers the multimodal tensor data classification problem. A Coupled Support Tensor Machine (C-STM) built upon the latent factors jointly estimated from the Advanced Coupled Matrix Tensor Factorization (ACMTF) is proposed. C-STM combines individual and shared latent factors with multiple kernels and estimates a maximal-margin classifier for coupled matrix tensor data. The classification risk of C-STM is shown to converge to the optimal Bayes risk, making it a statistically consistent rule. C-STM is validated through simulation studies as well as a simultaneous EEG-fMRI analysis. The empirical evidence shows that C-STM can utilize information from multiple sources and provide a better classification performance than traditional single-mode classifiers.


Introduction
Advances in clinical neuroimaging and computational bioinformatics have dramatically increased our understanding of various brain functions using multiple modalities such as Magnetic Resonance Imaging (MRI), functional Magnetic Resonance Imaging (fMRI), electroencephalogram (EEG), and Positron Emission Tomography (PET). Their strong connections to the patients' biological status and disease pathology suggest the great potential of their predictive power in disease diagnostics. Numerous studies using vector-and tensor-based statistical models illustrate how to utilize these imaging data both at the voxel-and Region-of-Interest (ROI) level and develop efficient biomarkers that predict disease status. For example, Anderson et al. (2011) proposes a classification model using functional connectivity MRI for autism disease and reaches 89% diagnostic accuracy for subjects under 20. Schindlbeck and Eidelberg (2018) utilizes network models and brain imaging data to develop novel biomarkers for Parkinson's disease. Many works in Alzheimer's disease research such as Morris et al. (2009);Liu et al. (2014); Khazaee et al. (2016); Gavidia-Bovadilla et al. (2017); Long et al. (2017); Ding et al. (2019); Li et al. (2020) use EEG, MRI and PET imaging data to predict patient's cognition and detect early-stage Alzheimer's diseases. Although these studies have provided impressive results, utilizing imaging data from a single modality such as individual MRI sequences are known to have limited predictive capacity, especially in the early phases of the disease. For instance, Li et al. (2020) uses brain MRI volumes from regions of interest to identify patients in early-stage Alzheimer's disease. They use one-year MRI data from Alzheimer's Disease Neuroimaging Initiative (ADNI) and obtain 77% prediction accuracy. Although such a performance is favorable compared to other existing approaches, the diagnostic accuracy is relatively low due to the limited information from MRI data. In recent years, it has been common to acquire multiple neuroimaging modalities in clinical studies such as simultaneous EEG-fMRI or MRI and fMRI. Even though each modality measures different biological signals, they are interdependent and mutually informative. Learning from multimodal neuroimaging data may help integrate information from multiple sources and facilitate biomarker developments in clinical studies. It also raises the need for novel supervised learning techniques for multimodal data in statistical learning literature.
Existing statistical approaches to multimodal data science are dominated by unsupervised learning methods. These methods analyze multimodal neuroimaging data by performing joint matrix decomposition and extracting common information across different modalities. During optimization, the decomposed factors bridging two or more modalities are estimated to interpret the connections between different modalities. Examples of these methods include matrixbased joint Independent Component Analysis (jICA) Calhoun et al. (2006); Groves et al. (2011); Lei et al. (2012); Liu et al. (2009) ;Sui et al. (2011); Acar et al. (2019) which assume bilinear correlations between factors in different modalities. However, these matrix-vector based models cannot preserve the multilinear nature of original data and the spatiotemporal correlations across modes as most neuroimaging modalities are naturally in tensor format. Recently, various coupled matrix-tensor decomposition methods have been introduced to address this issue Acar et al. (2014Acar et al. ( , 2017bAcar et al. ( , 2019; Chatzichristos et al. (2018Chatzichristos et al. ( , 2020; Karahan et al. (2015); Mosayebi and Hossein-Zadeh (2020). These methods impose different soft or hard multilinear constraints between factors from different modalities providing more flexibility in data modeling.
Current supervised learning approaches for multimodal data mostly concatenate data modalities as extra features without exploring their interdependence. For example, Zhou et al. (2013); Li et al. (2018) build generalized regression models by appending tensor and vector predictors linearly for image prediction and classification. Pan et al. (2018) develops a discriminant analysis by including tensor and vector predictors in a linear fashion.  proposes an integrative factor regression for multimodal neuroimaging data assuming that data from different modalities can be decomposed into latent factors. More recently, Gahrooei et al. (2021) proposed multiple tensor-on-tensor regression for multimodal data, which combines tensoron-tensor regression from Lock (2018) with traditional additive linear model. Another type of integration utilizes kernel tricks and combines information from multimodal data with multiple kernels. Gönen and Alpaydın (2011) provides a survey on various multiple kernel learning techniques for multimodal data fusion and classification with support vector machines. Combining kernels linearly or non-linearly instead of original data in different modalities provides more flexibility in information integration. Bach (2008) proposed a multiple kernel regression model with group lasso penalty, which integrates information by multiple kernels and selects the most predictive data modalities.
Despite these accomplishments, the current approaches have several shortcomings. First, they mainly focus on exploring the interdependence between multimodal imaging data, ignoring the representative and discriminative power of the learned components. Thus, the methods cannot further bridge the imaging data to the patients' biological status, which is not helpful in biomarker development. Second, the supervised techniques such as integrate information primarily by data or feature concatenation without explicitly considering the possible correlations between different modalities. This lack of consideration of interdependence may cause issues like overfitting and parameter identifiability. Third, even though methods from ; Gahrooei et al. (2021) have considered latent structures for multimodal data, these models are designed primarily for linear regression and are not directly applicable to classification problems. Fourth, the aforementioned multimodal analysis methods are mainly vector based methods, which cannot handle large-size multi-dimensional data encountered in contemporary data science. As discussed in Bi et al. (2020), tensors provide a powerful tool for analyzing multi-dimensional data in statistics. As a result, developing a novel multimodal tensor-based statistical framework for supervised learning can be of great interest. Finally, although many empirical studies demonstrate the success of using multimodal data, there is a lack of mathematical and statistical clarity to the extent of generalizability and associated uncertainties. The absence of a solid statistical framework for multimodal data analysis makes it impossible to interpret the generalization ability of a certain statistical model.
In this paper, we propose a two-stage Coupled Support Tensor Machine (C-STM) for multimodal tensor-based neuroimaging data classification. The proposed model addresses the current issues in multimodal data science and provides a sound statistical framework to interpret the interdependence between modalities and quantify the model consistency and generalization ability. The major contributions of this work are: (1) Individual and common latent factors are extracted from multimodal tensor data, for each sample or subject, using Advanced Coupled Matrix Tensor Factorization (ACMTF) (Acar et al., 2014(Acar et al., , 2017a. The extracted components are then utilized in a statistical framework. Most of the work on ACMTF do not work on each subject separately and the extracted factors are utilized for a signal analysis rather than a subsequent statistical learning framework. Specifically, the work on supervised approaches with CMTF is limited. (2) Building a novel Coupled Support Tensor Machine with both the coupled and noncoupled tensor CP factors for classification. In this regard, multiple kernel learning approaches are adopted to integrate components from multi-modal data.
(3) For the validation of our work, we provide both theoretical and empirical evidence. We provide theoretical results such as classification consistency for statistical guarantee. A thorough numerical study has been conducted, including a simulation study and experiments on real data to illustrate the usefulness of the proposed methodology.
A Matlab package is also provided in the supplemental material, including all functions for C-STM classification. The source codes are available at our Github repository 2 .

Related Work
In this section, we review some background and prior work on tensor decompositions and multiple kernel learning.

Notations
In this work, we denote numbers and scalars by letters such as , , . Vectors are denoted by boldface lowercase letters, e.g. , . Matrices are denoted by boldface capital letters like , .
Multi-dimensional tensors are denoted by boldface Euler script letters such as X X X, Y Y Y. The order of a tensor is the number of dimensions of the data hypercube, also known as ways or modes. For example, a scalar can be regarded as a zeroth-order tensor, a vector is a first-order tensor, and a matrix is a second-order tensor. Let X X X ∈ R 1 × 2 ×···× be a tensor of order , where 1 , 2 ,..., denotes the ( 1 , 2 , . . . , )th element of the tensor. Vectors obtained by fixing all indices of the tensor except the one that corresponds to th mode are called mode-fibers and denoted as where the mode-fibers of the tensor X X X are the columns of X X X ( ) and the remaining modes are organized accordingly along the rows.

Canonical/Polyadic (CP) Decomposition
Let X X X ∈ R 1 × 2 ×...× be a tensor with modes. Rank-Canonical/Polyadic (CP) decomposition of X X X is defined as: where ( ) ∈ R × , = 1, .., are defined as factor matrices whose columns are ( ) and "•" represents the vector outer product. The second equality in (1) is called Kruskal tensor (see Kruskal (1977)), which is a convenient representation for CP tensors. We denote a Kruskal tensor by = (1) , ..., ( ) or = ; (1) , ..., ( ) where ∈ R is a vector whose entries are the weights of rank one tensor components. In the special case of matrices, corresponds to singular values of a matrix. In general, it is assumed that the rank is small so that equation (1) is also called low-rank approximation for a tensor X X X. Such an approximation can be estimated by an Alternating Least Square (ALS) approach (see Kolda and Bader (2009)).
Although there are other tensor decomposition structures such as Tucker Decomposition or Tensor Train Decomposition, the advantage of CP is that the factors extracted by a CP decomposition are unique up to permutation. The uniqueness of the factors make CP decomposition more interpretable.

Coupled Matrix Tensor Factorization
Motivated by the fact that joint analysis of data from multiple sources can potentially unveil complex data structures and provide more information, Coupled Matrix Tensor Factorization (CMTF) (Acar et al., 2011) was proposed for multimodal data fusion. CMTF estimates the underlying latent factors for both tensor and matrix data simultaneously by taking the coupling between tensor and matrix data into account. This feature makes CMTF a promising model in analyzing heterogeneous data, which generally have different structures and modalities.
During latent factor estimation, CMTF solves an objective function that approximates a CP decomposition for the tensor modality and a singular value decomposition for the second modality with the assumption that the factors from one mode of each modality are the same. Given X X X 1 ∈ R 1 × 2 ×...× and 2 ∈ R 1 × 2 , without loss of generality assume that the factors X X X 1 2 3 1 2 Figure 1.: Illustration of Coupled Tensor Matrix Model from the first mode of the tensor X X X 1 span the column space of the matrix 2 . CMTF then tries to estimate all factors by minimizing: where ( ) are the factor matrices for modality and mode . The factor matrices (1) 1 = (1) 2 are the coupled factors between tensor and matrix data. An illustration of this coupling is given in Figure 1. These factor matrices can also be represented in Kruskal form, . By minimizing the objective function ( 1 , 2 ), CMTF estimates latent factors for the tensor and matrix data jointly which allows it to utilize information from both modalities. Acar et al. (2011) uses a gradient descent algorithm to optimize the objective function (2). Although this model is formulated for the joint decomposition of a th order tensor and a matrix, extensions to two or more tensors with couplings across multiple modes are possible.
In real data, couplings across different modalities might include shared or modality-specific (individual) components. Shared components correspond to those columns of the factor matrices that contribute to the decomposition of both modalities, while individual components carry information unique to the corresponding modality. Although CMTF provides a successful framework for joint data analysis, it often fails to obtain a unique estimation for shared or individual components. As a result, any further statistical analysis and learning from CMTF estimation will suffer from the uncertainty in latent factors. To address this issue, Acar et al. (2017a) proposed Advanced Coupled Matrix Tensor Factorization (ACMTF) by introducing a sparsity penalty to the weights of latent factors in the objective function (2), and restricting the norm of the columns of the factors to be unity to provide uniqueness up to a permutation. This modification provides a more precise estimation for latent factors compared to CMTF (Acar et al. (2017a,b)). In our framework, we utilize ACMTF to extract the latent factors which are in turn used to build a classifier for multimodal data.

CP-STM for Tensor Classification
CP-STM has been previously studied by Tao et al. (2005); He et al. (2014He et al. ( , 2017 and uses CP tensor to construct STM types of model. Assume there is a collection of data = {(X X X 1 , 1 ), (X X X 2 , 2 ), ..., (X X X , )}, where X X X ∈ X X X ⊂ R 1 × 2 ×...× are -way tensors. X X X is a compact tensor space, which is a subspace of R 1 × 2 ×...× .
∈ {1, −1} are binary labels. CP-STM assumes the tensor predictors are in CP format, and can be classified by the function which minimizes the objective function By using tensor kernel function where X X X 1 =

=1
(1) The STM classifier can be written as where X X X is a new -way rank-tensor, = [ 1 , ..., ] is the coefficient vector, is a diagonal matrix whose diagonal elements are 1 , .., and (X X X) = [ (X X X 1 , X X X), ..., (X X X , X X X)] is a column vector, whose values are kernel values computed between training and test data. We denote the collection of functions in the form of (5) with H H H , which is a functional space also known as Reproducing Kernel Hilbert Space (RKHS). The optimal classifier CP-STM ∈ H H H can be estimated by plugging function (5) into objective function (3), and minimize it with Hinge or Squared Hinge loss. The coefficients of the optimal CP-STM model is denoted by * . The classification model is statistically consistent if the tensor kernel function satisfying the universal approximating property, which is shown by Li and Maiti (2019).

Multiple Kernel Learning
Multiple kernel learning (MKL) creates new kernels using a linear or non-linear combination of single kernels to measure inner products between data. Statistical learning algorithms such as support vector machine and kernel regression can then utilize the new combined kernels instead of single kernels to obtain better learning results and avoid the potential bias from kernel selection (Gönen and Alpaydın (2011)). A more important and related reason for using MKL is that different kernels can take inputs from various data representations possibly from different sources or modalities. Thus, combining kernels and using MKL is one possible way of integrating multiple information sources.
Given a collection of kernel functions { 1 (·, ·), ... (·, ·)}, a new kernel function can be constructed by where is a linear or non-linear function and is a vector whose elements are the weights for the kernel combination. Linear combination methods are the most popular in multiple kernel learning, where the kernel function is parameterized as X X X 1 2 3 1 2 +..+ Individual Factors (Tensor Modality) The weight parameters can be simply assumed to be the same (unweighted) (Pavlidis et al. (2001); Ben-Hur and Noble (2005)), or be determined by looking at some performance measures for each kernel or data representation (Tanabe et al. (2008); Qiu and Lane (2008)). There are few more advanced approaches such as optimization-based, Bayesian approaches, and boosting approaches that can also be adopted (Lanckriet et al. (2004) 2002)). In this work, we only consider linear combination (7), and select the weight parameters in a heuristic data-driven way to construct our C-STM model.

Let
= {(X X X 1,1 , 1,2 , 1 ), ..., (X X X ,1 , ,2 , )} be training data, where each sample ∈ {1, . . . , } has two data modalities ,1 , ,2 and a corresponding binary label ∈ {1, −1}. In this work, following Acar et al. (2011), we assume that the first data modality is a third-order tensor, X X X ,1 ∈ R 1 × 2 × 3 , and the other is a matrix, ,2 ∈ R 4 × 3 . The third mode of X X X ,1 and the second mode of X X X ,2 are assumed to be coupled for each , i.e. the factor matrix is assumed to be fully or partially shared across these modes. Utilizing this coupling, one can extract factors that better represent the underlying structure of the data, and preserve and utilize the discriminative power of the factors from both modalities. Our approach, C-STM, consists of two stages: Multimodal tensor factorization, i.e., ACMTF, and coupled support tensor machine as illustrated in Figure 2. In this section, we present both stages and the corresponding procedures.

Multimodal Tensor Factorization
In this work, the first aim is to perform a joint factorization across two modalities for each training sample, . Let ,1 = ; (1) ,1 , (2) ,1 , (3) ,1 denote the Kruskal tensor of X X X ,1 , and ,2 = ; (1) ,2 , (2) ,2 denote the singular value decomposition of ,2 . The weights of the columns of each factor matrix ( ) , , where is the index for modality and denotes the mode, are denoted by and and the norms of these columns are constrained to be 1 to avoid redundancy. The objective function of ACMTF (Acar et al., 2014(Acar et al., , 2017a is then given by: (3) ,1 = (2) ,2 (1) where Σ Σ Σ is a diagonal matrix whose elements are the singular values of the matrix ,2 and ( ) , , ∈ R denotes the columns of the factor matrices for the object X X X , . The objective function in (8) includes penalties for the number of non-zero weights in both tensor and matrix decomposition. Thus, the model identifies the shared and individual components. These factors are then considered as different data representations for multimodal data, and used to predict the labels in C-STM classifier.

Coupled Support Tensor Machine (C-STM)
C-STM uses the idea of multiple kernel learning and considers the coupled and uncoupled factors from ACMTF decomposition as various data representations. As a result, we use three different kernel functions to measure their similarity, i.e., inner products. One can think of these three kernels inducing three different feature maps transforming multimodal factors into different feature spaces. In each feature space, the corresponding kernel measures the similarity between factors in this specific data modality. The similarities of multimodal factors are then integrated by combining the kernel measures through a linear combination. This procedure is illustrated in Figure 2. In particular, kernel 1 is a tensor kernel (equation (4)) since the first individual factors are tensor CP factors. The kernel function for C-STM is defined as for two pairs of decomposed tensor matrix factors ( ,1 , ,2 ) and ( ,1 , ,2 ). (3) * ,1, is the average of the estimated shared factors 1 2 [ (3) ,1, + (2) ,2, ]. 1 , 2 , and 3 are the three weights parameters combining three kernel functions, which can be tuned by cross-validation.
With the kernel function in (9), C-STM model tries to estimate a bivariate decision function from a collection of functions H H H such that where L (X X X , ) = max 0, 1 − (X X X ) · is Hinge loss. H H H is defined as the collection of all functions in the form of due to the well-known representer theorem (Argyriou et al. (2009)) for any pair of testing data (X X X 1 , 2 ) and for ∈ R . For all possible values of , equation (11) defines the data collection H H H . is a diagonal matrix whose diagonal elements are labels from the training data . (X X X 1 , 2 ) is a × 1 vector whose -th element is (X X X ,1 , ,2 ), (X X X 1 , 2 ) . The optimal C-STM decision function, denoted by = * (X X X 1 , 2 ), can be estimated by solving the quadratic programming problem where is the kernel matrix constructed by function (9). Problem (12) is the dual problem of (10), and its optimal solution * also minimizes the objective function (10) when plugging functions in the form of (11). For a new pair of test points (X X X 1 , 2 ), the class label is predicted as Sign (X X X 1 , 2 ) .

Model Estimation
In this section, we first present the estimation procedure for coupled tensor matrix decomposition (8), and then combine it with the classification procedure to summarize the algorithm for C-STM.
To satisfy the constraints in the objective function (8), the function ( ,1 , ,2 ) is converted to a differentiable and unconstrained form given by: where ℓ 1 norm penalties in (8) are replaced with differentiable approximations; and are Lagrange multipliers and > 0 is a very small number. This unconstrained optimization problem can be solved by nonlinear conjugate gradient descent (Acar et al. (2011(Acar et al. ( , 2014; Mosayebi and Hossein-Zadeh (2020)).
is a normalized matrix whose columns have unit ℓ 2 norms. By combining all of the partial derivatives, the partial derivative of the objective function is given by: which is a 2 + 5 dimensional vector. As mentioned in Acar et al. (2014), a nonlinear conjugate gradient method with Hestenes-Stiefel updates is used to optimize (13). The procedure is described in Algorithm 1.
Once the factors for all data pairs in the training set are extracted, we can create the kernel matrix using the kernel function in (9). By solving the quadratic programming problem Algorithm 1 ACMTF Decomposition 1: procedure ACMTF 2: Input: Multimodal data (X X X 1 , 2 ), r, , S (Upper limit for the number of iterations) 3: 11: 12: [ , ] = [ , ] 10: Solve the quadratic programming problem (12) and find the optimal * .

11:
Output: * (12), we can obtain the optimal decision function . This two-stage procedure for C-STM estimation is summarized in Algorithm 2.

Theory
In this section, we provide some preliminary theoretical results to validate the C-STM model. The first proposition provides a sketch of proof of convergence for the coupled matrix tensor decomposition.
Proposition 5.1. Suppose for every pair of multimodal data X X X ,1 and ,2 , the optimal latent factor estimate is the optimal solution for the objective function (8), which is denoted by * ,1 , * ,2 . The conjugate gradient descent algorithm 1 converges to stable estimates for tensor and matrix components * ,1 , * ,2 where: where is the number of iterations and (., .) is a distance measure between Kruskal tensor sets such as ℓ 2 distance between factors with an appropriate selection of permutations.
This proposition is a direct result of the convergence property of nonlinear conjugate gradient descent algorithm with line search (Nocedal and Wright (2006); Hager and Zhang (2006); Zoutendĳk (1970); Wolfe (1969Wolfe ( , 1971). The convergence rate of non-linear conjugate gradient descent is linear. For detailed convergence properties of nonlinear conjugate gradient descent with Hestenes-Stiefel updates on non-convex objectives, the readers are referred to Nocedal and Wright (2006); Powell (1984). The next result discusses the statistical property of C-STM. Let's assume the risk of a decision function, , is The expectation is taken over the joint distribution defined on X X X × Y Y Y, which is a data domain. The function 1 1 1{·} is an indicator function measuring the loss of classification function . It is also known as the "zero-one" loss since its value is zero when the decision function provides correct prediction and is one otherwise. If there is a * : X X X → Y Y Y from the collection of all measurable functions such that * = arg min R ( ), its risk is called the Bayes risk for the classification problem with data from X X X × Y Y Y. We denote the Bayes risk as R * = R ( * ). With different training sets , we can estimate a sequence of decision functions under the same training procedure. This sequence of decision function { } is called a decision rule. Obviously, C-STM is a decision rule if different training sets are provided. A decision rule is statistically consistent if R ( ) converges to the Bayes risk R * as the size of training data increases (see eg. Devroye et al. (2013)). The consistency property is desirable for classification rules, because a consistent rule guarantees to reconstruct the whole data distribution with more training data / observations. The reconstruction here means the Bayes risk of the classification problem will be eventually the same as the risk of estimated classifier with sufficient training data, and thus will be known. Our next result shows that C-STM is a statistically consistent decision rule.

AS.4
The hyper-parameter in the regularization term = satisfies: This proposition is an extension of our previous result for the statistical consistency of CP-STM. The proof of this proposition is provided in Appendix A.
For each simulation case, we generate 50 pairs of tensor and matrix from both class, collecting 100 pairs of observations in total. We then perform a random training and testing set separation by randomly choosing 20 samples as the testing set, and use the remaining data as the training set. The random selection of testing set is conducted in a stratified sampling manner such that the proportion of samples from each class remains the same in both training and testing sets. For all models, we report the model prediction accuracy, the proportion of correct predictions over total predictions, on the testing set as the performance metric. The random training and testing set separation is repeated for 50 times and the average prediction accuracy of these 50 repetitions for all the cases are reported in Figure 3. Additionally, the standard deviations are illustrated by the error bars in the figure. The results of CP-STM, CMDA, and DGTDA with tensor data are denoted by CPSTM1, CMDA1, and DGTDA1 respectively in the figure. The results using matrix data are denoted by CPSTM2, CMDA2, and DGTDA2.
From the Figure 3, we can conclude that our C-STM has a more favorable performance in this multimodal classification problem comparing with other competitors. Its accuracy rates are significantly larger than other methods in most cases. Particularly, we can see that the accuracy rates of C-STM (pink) are increasing from case 1 to case 5, while the accuracy rates of CP-STM using tensor data remain the same. This is because the difference between class mean vectors for the first tensor factor does not change from case 1 to case 5. However, the gap between class mean vectors in matrix factor increases. Due to this fact, both C-STM and CP-STM (yellow) which utilize matrix data are getting better performance from case 1 to case 5. More importantly, C-STM always outperforms CP-STM with matrix data as it enjoys the extra class information from multimodalities. In case 6 and case 7 where class information are in single data modalities, the advantage of C-STM is not as significant as the previous cases, though its performance are slightly better than CP-STM. This indicates C-STM can provide robust classification results even though extra data modalities do not provide any other class information. In case 8 where the class information is from the share factors, C-STM turns out to better recover the share factors and provides significantly better classification accuracy. Through this simulation, we showed that C-STM has a clear advantage of using multimodal data in classification problems, and are robust to redundant data modalities.

Trial Classification for Simultaneous EEG-fMRI Data
In this section, we present the application of the proposed method on simultaneous EEG-fMRI data. The simultaneous electroencephalography (EEG) with functional magnetic resonance imaging (fMRI) is one of the most popular non-invasive multimodal brain imaging techniques to study human brain function. EEG records electrical activity from the scalp resulting from ionic current within the neurons of the brain. Its millisecond temporal resolution makes it possible to record event-related potentials that occur in response to visual, auditory and sensory stimuli (Teplan et al. (2002); Abreu et al. (2018)). While EEG provides high temporal resolution, its spatial resolution is limited by the number of electrodes placed on the scalp and thus provides less spatial resolution compared to other neuroimaging modalities such as magnetic resonance imaging (MRI) and Positron Emission Tomography (PET). As a result, it has been commonplace to record EEG data in conjunction with a high spatial resolution modality. As another powerful tool in studying human brain function, blood oxygenation level dependent (BOLD) functional magnetic resonance imaging (fMRI) provides signals with much higher spatial resolution to reflect hemodynamic changes in blood oxygenation level at all voxels related to neuronal activities (Ogawa et al. (1990); Belliveau et al. (1991); Kwong et al. (1992); Filippi et al. (2005)). Recording simultaneous EEG and fMRI can provide high resolution information at both the spatial and temporal dimensions at the same time. Thus, developing novel machine learning techniques to utilize such multimodal data is of great significance. In this application, we apply our C-STM model to a binary trial classification problem on a simultaneous EEG-fMRI data.
The data is obtained from the study Walz et al. (2013). In this study, there are seventeen individuals (six females, average age 27.7) participated in three runs each of analogous visual and auditory oddball paradigms. The 375 (125 per run) total stimuli per task were presented for 200 ms each with a 2-3 s uniformly distributed variable inter-trial interval. A trial is defined as a time window in which subjects receive stimuli and make responses. In the visual task, a large red circle on isoluminant gray backgrounds was considered as the target stimuli, and a small green circle were the standard stimuli. For the auditory task, the standard and oddball stimuli were, respectively, 390 Hz pure tones and broadband sounds which sound like "laser guns". During the experiment, the stimuli were presented to all subjects, and their EEG and fMRI data are collected simultaneously and continuously. We obtain the EEG and fMRI data from OpenNeuro website (https://openneuro.org/datasets/ ds000116/versions/00003). We utilize both EEG and fMRI in this data set with our C-STM model to class stimulus types in all the trials. Through our numerical study, we want to demonstrate the fact our C-STM model enjoys the advantage of data multimodality and provides more accurate class predictions. The data from Subject 4 are dropped since its fMRI data are corrupted. Due to the fact that the number of trials from each subject are different, we further provide a table in appendix B to show the number of trials for each subject.
We pre-process both the EEG and fMRI data with Statistical Parametric Mapping (SPM 12) (Ashburner et al. (2014)) and Matlab. The EEG data is collected by a custom built MR-compatible EEG system with 49 channels. Walz et al. (2013) provides a version of rereferenced EEG data with 34 channels which are used in our experiment. This version of EEG data are sampled at 1,000 Hz, and are downsampled to 200 Hz at the beginning of pre-processing. We then remove both low-frequency and high-frequency noise in the data using SPM filter functions. As the last step of EEG pre-processing, we define trials from Brain Imaging Data Structure (BIDS) files Niso et al. (2018) and extract EEG data epochs recorded within the trial-related time windows. The time window for each trial is considered to go from 100 ms before the stimulus onset until 500 ms after the stimulus. For each trial, we construct a three-mode tensor corresponding to the EEG data for all subjects where the modes represent channel × time × subject. We denote it as X X X ,1 ∈ R 34×121×16 . The fMRI data is collected by 3T Philips Achieva MR Scanner with 170 volumes (TR = 2s) per session. Each 3D volume contains 32 slices. The voxel size in the image is 3 x 3 x 4 mm. For each subject, we realign all the fMRI volumes from multiple sessions to the mean volume, and co-register the participant's T1 weighted anatomical scan to the mean fMRI volume. Next, we normalize all the fMRI volumes to match the MNI brain template (Lancaster et al. (2007)) by creating segments from co-registered T1 weighted scan, and keep the voxel size as 3 x 3 x 4 mm. All normalized fMRI volumes are then smoothed by 3D Gaussian kernels with full width at half maximum (FWHM) parameter being 8 × 8 × 8. After the pre-processing, we further perform a regular statistical analysis (Lindquist et al. (2008); Worsley et al. (2002)) to extract fMRI volumes from visual and auditory stimulus related voxels. Such data are also known as Region of Interest (ROI) data. We describe the ROI voxel identification and data extraction in Appendix B. We extract fMRI volmes from 178 voxels (in Figure 4a) for auditory oddball tasks, and 112 voxels for auditory tasks. As a result, fMRI data are modeled by matrices whose rows and columns stand for voxels and subjects: X X X ,2 ∈ R 16×178 for auditory task data, and X X X ,2 ∈ R 16×112 for visual task data. There is no time mode in fMRI data because the trial duration is less than the repetition time of fMRI (time for obtaining a single 3D volume fMRI). For each trial, there is only one 3D scan of fMRI collected from a single subject. The ROI data then becomes a vector for this subject in the trial as we extract volumes from the regions of interest.
To classify trials with oddball and standard stimulus, we collect 140 multimodal data samples (X X X ,1 , X X X ,2 ) from auditory tasks, and 100 samples from visual tasks. For both types of tasks, the numbers of oddball and standard trials are equal. We consider the trials with oddball stimulus as the positive class, and the trials with standard stimulus as the negative class. Like the procedures in our simulation study, we select 20% of data as testing set, and use the remaining 80% for model estimation and validation. The classification accuracy, precision (positive predictive rate), sensitivity (true positive rate), and specificity (true negative rate) of classifiers are calculated using the test set at each experiment. The experiment is repeated multiple times, and the average accuracy, precision, sensitivity, and specificity, and their standard deviations (in subscripts) are reported in Table 2. The single mode classifiers CPSTM, CMDA, and DGTDA are also applied on either EEG or fMRI data as a comparison. The single mode classifiers applied on EEG data are denoted by appending the number "1" after their names, and those applied on fMRI data are denoted by appending the number "2". The area under the curve (AUC) for all the classifiers are also reported in Table 2.
The results in Table 2 show that the trial classification accuracy for C-STM using multimodal data is better than any classifier based on single modality with a significant improvement in  terms of average accuracy rates and average AUC values. This improvement is observed for classification of both auditory and visual tasks. This observation agrees to the conclusion from our simulation study. Similar to our simulation study, the tensor discriminant analysis does not work as well as CP-STM and C-STM. In addition, it is obvious that the performance of tensor discriminant analysis using fMRI data are better than using EEG data. This is within our expectation, since the regions we extracted from fMRI data are identified by group level fMRI statistical analysis (see Appendix B). The data in these regions have already shown significant differences between different trials in the traditional study, and thus are easy to classify. On the other hand, there is no prior analysis and feature extraction procedure applied on EEG data, leaving a low signal to noise ratio in EEG data. However, C-STM still can take advantage of using EEG data and further increase the classification accuracy, highlighting its robustness and potential in processing noisy mulitmodal tensor data.

Conclusion
In this work, we have proposed a novel coupled support tensor machine classifier for multimodal data by combining the advanced coupled matrix tensor factorization (ACMTF) and support tensor machine (STM). The most distinctive feature of this classifier is its ability to integrate features across different modalities and structures. The proposed approach can simultaneously take matrix-and tensor-shaped data for classification and can be easily extended to inputs with more than two modality. The coupled tensor matrix decomposition unveils the intrinsic correlation structure between data across different modalities, making it possible to integrate information from multiple sources efficiently. Such a decomposition also makes the whole method robust and applicable to large-scale noisy data with missing values. The newly designed kernel functions in C-STM provide feature-level information fusion, combining discriminant information from different modalities. Moreover, the kernel formulation makes it possible to utilize the most discriminative features from each modality by tuning the weight parameters in the function. Our theoretical results demonstrate that the C-STM decision rule is statistically consistent.
The most important theoretical extension of our current approach would be the development of excess risk for C-STM. In particular, we are looking for an explicit expression for the excess risk in terms of data factors from multiple modalities to quantify the contribution of every single modality in minimizing the excess risk. By doing so, we are able to interpret the importance of each data modality in classification tasks. In addition, quantifying the uncertainty of tensor and matrix factors estimation and their impact on the excess risk will build the foundation to the next level.
Future work will focus on learning the weight parameters in the kernel function via optimization. As Gönen and Alpaydın (2011) introduced, the weights in the kernel function can be further estimated by including a group lasso penalty in the objective function. Such a weight estimation procedure can identify the most significant data components and reduce the burden of parameter selection. In addition, the proposed framework can be extended to multimodal tensors with more than two modalities, and for regression problems.
In conclusion, we believe C-STM offers many encouraging possibilities for multimodal data integration and analysis. Its capability of handling multimodal tensor inputs will make it appropriate in many advanced data applications in neuroscience and medical research. We anticipate that this method will play an important role in a variety of applications.

Appendix A. Proof of Theorem 5.2
Proof. To prove the proposition 5.2, we introduce few more notations here. Let L be the loss function satisfying the condition AS.2. We denote the classification risk for an arbitrary decision function, , as The expectation is taken over the joint distribution of X X X × Y Y Y. Notice that this risk notation, R L ( ), is different from our notation R ( ) in section 5 since we use the Lipschitz continuous loss L instead of the "zero-one" loss to measure the classification error. L is also called surrogate loss for classification problems. Examples of such surrogate loss functions include Hinge loss and Squared Hinge loss. Comparison of these loss functions and their statistical properties can be found in Zhang et al. (2004). If we denote the Bayes risk under the surrogate loss L as R * L , i.e. R * L = min R L ( ) for all measurable function , then the result from Zhang et al. (2004) says R L ( ) → R * L indicates R ( ) → R * for any decision rule { }. This conclusion holds as long as the surrogate loss is "self-calibrated" (see Steinwart and Christmann (2008)). Since we use Hinge loss in our problem, and Hinge loss is known to be Lipschitz and self-calibrated, our assumption AS.2 holds in our discussion. Thus, we only need to show R L ( ) → R * L for the proof of our proposition 5.2. Given the tuning parameter satisfying condition AS.4, we denote where H H H is the reproducing kernel Hilbert space (RKHS) generated by the kernel function (9). As we mentioned in the section 3, H H H is also know as the collection of functions which are in the form of equation (11). Now we further assume = arg min ∈H H H · || || 2 + R L ( ) Then is the optimal decision function from H H H such that it minimizes the expected risk. Comparing with , we can understand that is the version of when the size of training data is as large as possible. If we denote R L, ( ) = 1 =1 L ( (X X X ), ), then R L, ( ) is a sample estimate of R L ( ). With , we can show that through triangular inequality. Since the Bayes risk under loss function L is defined as R * = min :X X X→Y Y Y R ( ) over all functions defined on X X X, we can immediate show that This is the result of using condition AS.1 and AS.2 in the proposition 5.2. is in the RKHS and thus bounded by some constant depending on . * * * is also continuous on compact subspace X X X (because all the tensor components considered are bounded in condition AS. → 0, making → ∞. Further notice that we introduce since it is independent of . As a result, its norm, even though is bounded by , is a constant and is not changing with respect to . By Rademacher complexity, the following inequality holds with probability at least 1 − , where 0 < < 1 is not the optimal in training data √︂ log 2/ 2 Rademacher Complexity again R L ( ) + || || 2 + 4 ( ) √ + 2( 0 + ( ) ) √︂ log 2/ 2 Let = 1 2 , and large such that for all > , The inequality exists because || || is a constant with respect to , and all other terms are converging to zero. Thus R L ( ) R L ( ) + 2 R L ( ) + with probability 1 − 1 2 . We conclude that for any arbitrary . This establishes the weak consistency of CP-STM. For strong consistency, we consider for each (2019)), R L ( ) → R L ( ) almost surely. The proof is finished.

Appendix B. Data Pre-processing for Section 7
We provide further details about our EEG-fMRI data pre-processing and fMRI data extraction in this section. Most of the processing steps are referred from Henson et al. (2019).

B.1. fMRI Data
The fMRI data processing includes three major steps, which are pre-processing, regions of interests (ROI) identification, and data extraction. We describe all these steps here. All the steps are performed by SPM 12 in Matlab. There are five steps in the image pre-processing part including realignment, co-registration, segment, normalization, and smoothing.
• Realignment: It is a procedure to align all the 3D BOLD volumes recorded along the time to remove artifacts caused by head motions, and also to estimate head position. For each task, there are three sessions of fMRI scans, providing 510 scans in total for each subject. These scans are realigned within subject to the average of these 510 scans. (average across time) In SPM, we create three independent sessions to load all the fMRI runs, and choose not to reslice all the images at this step. The reslicing will be done in normalization step. Avoiding extra reslicing can avoid introducing new artifacts. The mean scan is created in this step for co-registration. • Co-registration: Since all the fMRI scans are aligned to the mean scan, we have to transform the T1 weighted anatomical scan to match their orientation. Reason for doing this is that all the data will finally be transformed to a standardized space. Estimating such a transformation with T1 weighted scan can provide a high accuracy, since anatomical scans have higher resolutions. Matching the orientation of T1 weighted scan with all the fMRI scans makes it possible to apply the transformation estimated from T1 scan directly on fMRI data. In this step, we let the mean fMRI scan to be stationary, and move T1 anatomical scan to match it. A resliced T1 weighted scan is created in this step. • Segment: This step estimate a deformation transformation mapping data into MNI 152 template space Lancaster et al. (2007); Brett et al. (2001). A forward deformation field is created in this step. • Normalization: In this step, the forward deformation is applied to all realigend fMRI scans, transforming all the data into MNI template space. The voxel size is set to be 3 × 3 × 4 mm, which is the same as the original images. This pre-processing procedure is applied to auditory and visual fMRI scans separately and independently. For each task, the processed fMRI are used to for statistcal analysis introduced in Lindquist et al. (2008); Worsley et al. (2002). These models are basic linear mixed effect model with auto-regression covariance structure. Since these models are standard and are out of the scope of this dissertation, we do not introduce them in this part. For the first level (subject level) analysis, we use the model to estimate two contrast images: standard stimulus over baseline and oddball stimulus over baseline. These two are difference of average BOLD signals during stimulus time and that during no stimulus (baseline) time. They can be understand as the estimateˆin a regression model = + . These contrasts are then pooled together in the group-level analysis. For each voxel, the group-level analysis performs a T-test to compare the BOLD signals in standard contrasts and oddball contrasts. For voxels whose test results is significant, SPM highlighted them as the regions of interest (ROI). The ROI of auditory and visual tasks are presented in the figure B1 and figure B2 with P-values. Figure B1.: Auditory fMRI Group Level Analysis The coordinates of these activate voxels are also provided in the statistical analysis results. To extract ROI data, we can use "spm_get_data" function in SPM 12. Since we are classifying trials, we only take one fMRI scan for each trial. This is because the trial duration (0.6 sec) is less than the repetition time (2 sec) of fMRI data. For each trial, we take the -th fMRI scan Figure B2.: Visual fMRI Group Level Analysis where "k = round(onset / TR) + 1". This option is also inspired by SPM codes.

B.2. EEG Data
The pre-processing of EEG data is relatively easy comparing to fMRI, since all EEG are already converted to MAT file and are re-referenced. Thus, we only need to resample it using Matlab Signal Processing toolbox to a lower sampling rates, which is 200 Hz in our case. Then, we use function "ft_preproc_lowpassfilter" and "ft_preproc_highpassfilter" from SPM 12 toolbox to filter the data. Finally, we split EEG into epochs which starts 100 ms before the onset and ends 500 ms after the onset. According to Henson et al. (2019), such a duration is long enough to capture the event-related potential for EEG data.

C.1. Multimodal Tensor Factorization
The proposed model requires the selection of three different parameters, namely, , , and rank . To select these parameters, we closely follow best practices outlined in previous work on CMTF (Acar et al., 2011), ACMTF (Acar et al., 2014) and CCMTF (Mosayebi and Hossein-Zadeh, 2020). First of all, one of these parameter can be set to 1 as a pivot, and following previous work, we set = 1. The selection of rank is directly related to the selection of . As enforces sparsity over the singular values, it directly minimizes the rank. With sufficiently large , we can estimate the low-rank part through optimization. For the selection of in real data, we set = 5 following the work of Mosayebi and Hossein-Zadeh (2020) where it was shown through CORCONDIA tests (Kamstrup-Nielsen et al., 2013) that = 3 is sufficiently large for oddball data. In the case of the simulation study, = 5 is again sufficiently large as the data were generated from rank = 3 factors. Finally, based on our empirical results and the results presented in (Mosayebi and Hossein-Zadeh, 2020) we set = 0.001 using k-fold cross validation.

C.2. C-STM
The parameters in C-STM include kernel weights 1 , 2 , 3 and regularization parameter in the optimization. The weight parameters are selected based on the procedure in Gönen and Alpaydın (2011), and is selected using cross-validation.