{'fraction_non_alphanumeric': 0.028902232309592896, 'fraction_numerical': 0.008532766397563708, 'mean_word_length': 5.074694589877836, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MULTIPAR: Supervised Irregular Tensor Factorization with Multi-task Learning

Tensor factorization has received increasing interest due to its intrinsic ability to capture latent factors in multi-dimensional data with many applications such as recommender systems and Electronic Health Records (EHR) mining. PARAFAC2 and its variants have been proposed to address irregular tensors where one of the tensor modes is not aligned, e.g., different users in recommender systems or patients in EHRs may have different length of records. PARAFAC2 has been successfully applied on EHRs for extracting meaningful medical concepts (phenotypes). Despite recent advancements, current models' predictability and interpretability are not satisfactory, which limits its utility for downstream analysis. In this paper, we propose MULTIPAR: a supervised irregular tensor factorization with multi-task learning. MULTIPAR is flexible to incorporate both static (e.g. in-hospital mortality prediction) and continuous or dynamic (e.g. the need for ventilation) tasks. By supervising the tensor factorization with downstream prediction tasks and leveraging information from multiple related predictive tasks, MULTIPAR can yield not only more meaningful phenotypes but also better predictive performance for downstream tasks. We conduct extensive experiments on two real-world temporal EHR datasets to demonstrate that MULTIPAR is scalable and achieves better tensor fit with more meaningful subgroups and stronger predictive performance compared to existing state-of-the-art methods.


Introduction
Tensor factorization has received increasing interest due to its intrinsic ability to capture the multi-dimensional structure in the data. It has a wide range of applications including social network analysis [1,2], health data mining [3,4,5,6,7,8,9], recommender systems [10], and signal processing [11]. Canonical Polyadic (CP) [12,13], Tucker [14], and tensor singular value decomposition (SVD) [15,16] are popular regular tensor factorization methods, where each mode of the tensor has fixed size. However, in real-world cases, different people in recommender systems or patients in health data may have different length of records, which can not be handled by regular tensor factorization methods. PARAFAC2 [17,18] has been proposed for factorizing irregular tensors, where one of the mode size is not fixed. We introduce our motivating application for Electronic Health Records (EHRs) mining below for irregular tensor factorization.
EHRs are patient-centered records collected from a variety of institutions, hospitals, pharmaceutical companies across a long period of time. EHR mining can significantly improve the ability to diagnose diseases and reduce or even prevent medical errors, thereby improving patient outcomes [19]. However, directly using the raw, massive, high-dimensional, patient Figure 1: PARAFAC2 tensor factorization and longitudinal EHRs is challenging. For example, a single disease can consist of several heterogeneous subgroups yet be coded with the same diagnosis category, e.g., hypertension can be divided into hypertension resolver, hypertension, and prehypertension subgroups. Thus, researchers and healthcare practitioners seek to identify phenotypes, or disease subgroups, to better understand differences in biological mechanisms and treatment responses, which can lead to more effective and precise treatment.
Tensor factorization has been widely used to capture the multi-dimensional structure in EHRs for computational phenotyping [20,21,3,4,5,6,7,8,9]. Compared to traditional clustering-based approaches, tensor factorization can mine concise and potentially more interpretable latent information between multiple attributes (e.g., diagnosis and medications) in addition to clustering patients into subgroups. One key characteristic of EHR data is that different patients have different visit lengths, varying disease states, and varying time gaps between consecutive visits. Hence PARAFAC2 has been applied to extract pheonotypes from EHR data.
EHR Example. Table 1 shows an EHR database with a set of patient records. Each record consists of a sequence of visits. Each visit contains measurements for a set of medical features (e.g., laboratory results, medications, etc.). Circles denote missing observations. The number of visits can be different across patients.
Each patient record can be captured using a binary, numeric, or count matrix X k , where each matrix value represents the measurement associated with a particular feature for a particular visit. The entire data can be represented as an irregular tensor as shown in Figure 1 where each slice X k represents the information of patient k with I k visits and J medical features. Figure 1 illustrates the computational phenotyping process using PARAFAC2. Each slice of the irregular tensor X will be factorized by PARAFAC2 to three factor matrices. U k ∈ R I k ×R captures temporal evolution of the R phenotypes for patient k. V ∈ R J×R contains the R phenotypes. Each row of V matrix represents one latent and potentially interpretable phenotype. Each medical feature is represented with a weight indicating its contribution to the phenotype in each row. S k ∈ R R×R is a diagonal matrix with the importance membership of patient k in each one of the R phenotypes. The right side table in Figure 1 shows three example phenotypes represented in the V matrix with the top 4 highest weighted medical features in each phenotype (the example is from the MIMIC-EXTRACT dataset [22] and the associated weight is shown in parenthesis). Each phenotype represents a set of related medical features which can suggest meaningful subgroups.
Various works have proposed improvements to the basic PARAFAC2 model. SPARTan [23] modified the MTTKRP calculation order of PARAFAC2 to handle large and sparse data. COPA [9] further introduced various constraints to improve the interpretability of the factor matrices for more meaningful phenotype extraction. REPAIR [7] and LogPar [8] added a low-rankness constraint to improve PARAFAC2 robustness to missing data. Despite these advances, current PARAFAC2 models are completely unsupervised and only attempt to learn the latent factors to best recover the original observations. Some works have considered using the latent factors as features for downstream prediction  Figure 2: Overview of MULTIPAR on MIMIC-EXTRACT dataset tasks (e.g., in-hospital mortality or hospital readmission prediction using extracted phenotypes), and achieved limited performance gain than using the raw data as features. This is because the tensor factorization does not take advantage of the downstream labels. The extracted factors, while interpretable, may not be the most representative or discriminating for downstream prediction tasks. In addition, current work [9,8,7,23,24] using tensor factorization for predictive tasks only consider a single task (e.g., in-hospital mortality prediction) and ignore useful information from other prediction tasks.
To address these limitations, we propose MULTIPAR: a supervised irregular tensor factorization with multi-task learning for both phenotype extraction and predictive learning, as shown in Figure 2. MULTIPAR jointly optimizes the tensor factorization and downstream prediction together, so that the factorization can be "supervised" or informed by the predictive tasks. In addition, we use a multi-task framework to leverage information from multiple predictive tasks. It provides flexibility to incorporate both one-time or static (e.g. in-hospital mortality prediction) and continuously changing or dynamic (e.g. the need for ventilation) outcomes. To achieve this, the temporal features from U matrix are used for dynamic prediction and the features from S matrix are used for static prediction, as shown in Figure 2.
Our main hypothesis is that such a supervised multi-task framework can yield not only more meaningful phenotypes but also better predictive accuracy than performing tensor factorization independently followed by predictive learning using the phenotypes extracted from the tensor. In addition, by sharing the representation across tasks, the learned phenotypes can generalize better for each task. Our empirical studies on two large publicly available EHR datasets with representative predictive tasks (both static and dynamic) and different models (e.g. logistic regression and recurrent neural networks) verified this hypothesis.
In summary, we list our main contributions below: 1. We propose a supervised framework for PARAFAC2 tensor factorization and downstream prediction tasks such that the factorization can be "supervised" or informed by the predictive tasks. 2. We use a multi-task framework to leverage information from multiple predictive tasks and provide flexibility to incorporate both static and dynamic tasks and different models (e.g. logistic regression and recurrent neural networks). 3. We introduce a novel unified and dynamic weight selection method for weighing the tensor factorization and predictive tasks during the optimization process, where the tensor factorization is considered as one task, to achieve an overall optimized result. 4. We evaluate MULTIPAR's tensor reconstruction quality, predictability, scalability, and interpretability on two real-world temporal EHR datasets through a set of experiments, which verify MULTIPAR can identify more meaningful subgroups and yield stronger predictive performance compared to existing state-of-the-art approaches.

Irregular Tensor Factorization
PARAFAC2 model is the state-of-the-art tensor factorization framework for irregular tensor, i.e., tensors that do not align along one of its modes. The PARAFAC2 model decomposes each frontal slice of the irregular tensor X as U k S k V ⊤ , where U k ∈ R I k ×R , S k ∈ R R×R is diagonal and V ∈ R J×R . R is the target rank. Uniqueness is an important property in tensor factorization models that ensures the solution is not an arbitrarily rotated version of the PARAFAC2 model

Scalability Predictability Robustness Interpretability
Original PARAFAC2 x x x x SPARTan [23] x x x COPA [9] x x REPAIR [7] x LogPar [8] x x MULTIPAR x Table 2: Comparison of existing PARAFAC2-based models actual latent factors. In order to enforce uniqueness, Harshman [13] imposed the constraint U T U k = Φ, ∀k. This is equivalent to each U k being decomposed as U k = Q K H, where Q k ∈ R I k ×R , Q ⊤ k Q k = I, and H ∈ R R×R . Given the above modeling, the standard algorithm to fit PARAFAC2 solves the following optimization problem: Given a tensor representing the EHR data as in Figure 1 where each slice X k represents the information of patient k with I k visits and J medical features, PARAFAC2 decomposes the irregular tensor X into the factorization matrices which have the following interpretations: • U k ∈ R I k ×R represents the temporal trajectory of I k clinical visits in each one of the R phenotypes.
• V ∈ R J×R represents the relationship between medical features and phenotypes.
• S k ∈ R R×R represents the relationship between patients and phenotypes. Each column in S represents one phenotype, and if a patient has the highest weight in a specific phenotype, it means the patient is mostly associated with or exhibits a particular phenotype.
SPARTan [23] was developed to decompose large sparse datasets, and COPA [9] extended SPARTan to enhance the interpretability by adding more constraints, e.g., smoothness on U k and sparsity on V. REPAIR [7] and LogPar [8] added low-rankness constraints to improve the robustness of PARAFAC2 model to handle missing values. However, no current work has considered using downstream tasks to supervise and improve the predictability of PARAFAC2 model as Table 2 shows. We note that while we focus on introducing the predictability of MULTIPAR in this paper, the framework can be easily extended to incorporate robustness by adding low-rankness constraint on factor matrices.

Supervised and Multi-task learning Framework
Supervision can enhance traditional unsupervised tasks such as clustering for a wide range of applications, e.g., graph learning [25], pattern classification [26]. Wang et al. proposed a supervised feature extraction framework using discriminative clustering to improve graph learning model's clustering accuracy [25]. Liu te al. proposed a supervised minimum similarity projection framework using lowest correlation representation to improve pattern model's classification accuracy [26].
Over the past years, multi-task learning (MTL) [27,28] has attracted much attention in the artificial intelligence and machine learning communities. Traditional machine learning frameworks solve a single learning task each time, which ignores commonalities and differences across different tasks. MTL aims to learn multiple related tasks jointly so that the knowledge contained in one task can be leveraged by other tasks, with the hope of improving generalization performance by learning a shared representation [29,30]. MTL has been used successfully across all applications, from natural language processing [31,32] and speech recognition [33,34] to computer vision [35] and drug discovery [36,37]. However, no current work has considered improving predictability of tensor factorization using MTL.

Proposed Method
In this section, we present the MULTIPAR model in the context of EHR pheynotyping and its optimization. The general framework is applicable to any irregular tensor factorization and predictive learning tasks such as recommender systems.

Problem Formulation
We formalize the objective function for the MULTIPAR model in Definition 2. The PARAFAC2 loss for X ensures the reconstructed tensor closely approximates the original tensor. The static outcomes loss and dynamic outcomes loss are separate prediction tasks. Static outcome prediction tasks have a one-time or static labels, and dynamic outcome prediction tasks have a continuously changing or temporal dynamic labels for each time stamp. An approximate uniqueness constraint ensures tensor factorization uniqueness. For EHRs phenotype discovery, various constraints can be imposed on the factorization matrices to yield meaningful and high-interpretability phenotypes. The MULTI-PAR model accommodates such interpretability-purposed constraints in eq. (1) including: non-negativity for c 1 (S k ), sparsity for c 2 (V). We explain each of the loss components and constraints in detail below.
PARAFAC2 loss. The PARAFAC2 tensor factorization loss can ensure the reconstructed tensor closely approximate the original tensor. To accommodate different data types, the PARAFAC2 loss can be any smooth loss function, e.g., Least square loss, Poisson loss [38] and Rayleigh Loss [38].
Static outcomes loss. Previous PARAFAC2 models separate the PARAFAC2 training process and downstream prediction process. For example, in-hospital mortality prediction accuracy may be used as the metric to measure the predictability of the phenotypes extracted by the model. In the MULTIPAR model, we optimize the downstream prediction tasks and tensor factorization together by adding the prediction losses of the prediction tasks to the objective function. If the prediction task has one label per patient, we denote it as a static outcome prediction task. For illustrative purposes, we use a logistic regression model on the S matrix to predict static outcome tasks, and add the cross-entropy loss to the objective function. In fact, any differentiable loss function (e.g., square loss, exponential loss) can be incorporated in the objective function.
Dynamic outcomes loss. Different from static outcomes, dynamic outcomes have labels at each timestamp. For example, predicting whether a patient will be on a ventilator at a given future time can be used to measure the model's predictability. For illustrative purposes, we use the long short-term memory (LSTM) model on the U matrix to predict each patient's dynamic outcome labels, and add the loss of the LSTM model to the objective function. Similar to the static outcome loss, other models (e.g., gated recurrent units, vanilla recurrent neural networks) and their associated loss functions can be incorporated in the objective function.
Approximate uniqueness constraint. The optimization of the original PARAFAC2 model adopts the alternating direction method of multipliers (AO-ADMM) [39], which can not make full use of the parallel computation feature of GPUs. To adopt mainstream deep learning frameworks like Pytorch and Tensorflow, we use a stochastic gradient descent (SGD) based optimization approach. The uniqueness constraint in the original PARAFAC2 model is Q ⊤ k Q k = I. Similar to LogPar [8], to optimize Q k we relax the uniqueness constraint to Q ⊤ k Q k − I 2 F . Sparsity on V. The V matrix captures the association between a medical feature and a particular phenotype. In order to improve interpretability, we introduce a sparsity constraint on the V matrix. l 0 and l 1 norms are two popular sparsity regularization techniques. The l 0 regularization norm, also relaxed by hard thresholding, will cap the number of nonzero values in a matrix. The l 1 regularization norm, also relaxed by soft thresholding, will shrink matrix values towards zero. As hard thresholding is a non-convex optimization problem which can not be optimized by the SGD framework, we adopt the soft thresholding, which is convex and can be migrated into the SGD optimization framework.
Non-negativity on S. The diagonal matrix S indicates the importance membership of patient k in each one of the R phenotypes. Since only non-negative membership values makes sense, we zero out the negative values in S. SDW: Smooth dynamic weight selection. Numerous deep learning applications benefit from MTL with multiple regression and classification objectives. Yet the performance of MTL is strongly dependent on the relative weighting between each task's loss. Our objective function consists of several losses from the tensor factorization and the predictive tasks. Each of this loss is associated with a weight. While Definition 2 shows ρ 1 as the weight for tensor loss, ρ 2 and ρ 3 as the accumulative weights for static and dynamic tasks respectively, here we use W eight n (t) to denote the weight for each individual task n in epoch t. A key challenge is how to tune these weights for different tasks.
While the DWA weight selection [40] was proposed to dynamically change the task weights at each epoch by considering the rate of change of the loss over the epoch, the noisy nature of SGD weights can cause drastic fluctuations in the task weights between epochs. This can cause oscillating behavior between the various tasks and impedes convergence of the algorithm. Therefore, we propose a novel smooth dynamic weight selection method to choose the weight for each task. We first calculate the relative descending rate of each task loss and denote it as ω n (t − 1). t here represents an epoch index: We then calculate the weight for each task using the following equation: Intuitively, each task weight is dynamically updated based on a "smoothed" descending rate of the loss, the higher the rate (i.e. the more the task contributes to the optimization objective), the higher the weight for the task in next epoch. We use C to control the softness distribution between different tasks. If C is large enough, the weight for each task will be uniform. Different from [41], we introduce m, the weight update window size. The task weights are updated based on an average of the descending rate of the loss over m past epochs from t − m to t − 1 (instead of using one previous iteration). The main rationale for this smoothing is to reduce the SGD update uncertainty and training data selection randomness. Finally, a softmax operator, which is multiplied by the number of tasks N , ensures the sum of the weight equals N . For t = 1, we initialize all the weights to 1.

Optimization
To solve the optimization problem in Eq. (1), MULTIPAR follows an alternative optimization strategy where we optimize one variable individually with all other variables fixed. According to the subproblem differentiable, we group the variables into two groups: pure smooth subproblems which can be directly solved by SGD and proximal mappingbased smooth subproblems [42]. Proximal map is a key building block for optimizing nonsmooth regularized objective functions, e.g., the · 1 ℓ 1 -norm regularization function for inducing sparsity and · * nuclear norm regularization for inducing low-rankness. In the following, we omit the iteration number for brevity in notation.

Pure Smooth Subproblems Updates.
For the pure smooth subproblems, we use SGD to update the variables, which include the following three parts: Update of U k . The subproblem of U k takes the form as follows arg min Update of Q k . The subproblem of Q k takes the form as follows arg min Update of H. The subproblem of H takes the form as follows

Proximal Mapping-base Smooth Subproblems Updates
For the nonsmooth subproblems, we propose a proximal mapping-based [42] update, which include the following two parts.
Update of S k . The subproblem of S k takes the form as follows arg min We use projected SGD to update S k , where each step takes the following form where G[S k ] denotes the stochastic gradient of the smooth part Update of V. The subproblem of V takes the form as follows We use soft-thresholding operator to update V, where each step takes the following form: where λ is the step-size and G[V] denotes the stochastic gradient of the smooth part K k=1 The complete algorithm. The optimization procedure is summarized in Algorithm 1.

Complexity and Convergence Analysis
The following theorem summarizes the computational complexity of Algorithm 1. Proof. MULTIPAR's per-iteration complexity breaks down as follows: Line 2 costs O(R(R + J + K)); Line 3 costs O(min{R 2 I, RI 2 }), where I denotes the maximum among {I k }; Line 4,5,6 cost O(R 2 (R + J + K)). As a result, the per-iteration complexity is O(3R 2 JK). Proof. (Sketch of proof) MULTIPAR is a multi-block SGD algorithm for nonconvex optimization in nature, which has been extensively studied in optimization literature. The above theorem follows the convergence results established in [43].

Dataset
We use two real-world datasets to evaluate MULTIPAR in terms of its reconstruction quality, predictive performance, interpretability, and scalability. eICU 1 [44]: The eICU Collaborative Research Database is a freely available multi-center database for critical care research. It contains variables used to calculate the Acute Physiology Score (APS) III for patients. We select 202 diagnosis codes that have the highest frequency, as in [45]. The resulting number of unique ICU visits is 145426. The maximum number of observations for a patient is 215. We select three static outcome prediction tasks, including intubated prediction, ventilation prediction, and dialysis prediction. The ventilation prediction here is a static prediction task indicating whether a patient needs to be ventilated at the time of the worst respiratory rate, we will use "vent-res" as the name for this task.
MIMIC-EXTRACT 2 [22]: MIMIC-Extract is an open-source pipeline for transforming raw EHR data in MIMIC-III into data frames that are directly usable in common machine learning pipelines. We use the vitals labs mean table, which contains 34,472 patients with 104 features (Vital lab codes). The maximum number of observations for a patient is 240. We further normalize the data to [0,1]. We select three static outcome prediction tasks, including in-hospital mortality prediction, readmission prediction, ICU mortality prediction, and one dynamic outcome prediction task, which is ventilation prediction for every visit.

Evaluation Metrics
In order to test the tensor reconstruction quality of MULTIPAR model, we adopt the F IT ∈ (−∞, 1] score [46] as the quality measure (the higher the better): The original tensor, denoted as {X k }, serves as the ground truth. U k , S, V are factor matrices after the MULTIPAR tensor factorization. We evaluate the derived phenotypes' predictability power using the PR-AUC score of the prediction tasks. We split the data with a proportion of 8:2 as training and test sets and use PR-AUC score to evaluate the predictive power.

Methods for Comparison
We compare MULTIPAR with three baseline methods: SPARTan, COPA, and singlePAR. SPARTan and COPA are two state-of-the-art irregular tensor factorization methods. We also compare against a supervised single task PARAFAC2, which is a single-task version of MULTIPAR.
• SPARTan [23] -scalable PARAFAC2: A tensor factorization method for fitting large and sparse irregular tensor data. It only considers the tensor reconstruction loss.
• COPA [9] -scalable PARAFAC2 with additional regularizations: An irregular tensor factorization method that introduces various constraints/regularizations to improve the interpretability of the factor matrices. For both SPARTan and COPA, the extracted phenotypes are used for training the models for the downstream predictive tasks.
• SinglePAR -supervised single task PARAFAC2: The supervised irregular tensor factorization with single prediction task (single task version of MULTIPAR). The weight of tensor factorization and prediction tasks are also tuned using SDW.

Implementation Details
The multi-task dynamic weight selection of MULTIPAR has two hyper-parameters that need to be tuned. C controls the softness distribution between different tasks, and m is the weight update window size. In order to find the best C, we vary C from 0.2 to 2, and compare the prediction tasks' PR-AUC scores on different data set under different ranks.

Experiment Result
Tensor reconstruction quality analysis. For the following experiments on tensor reconstruction quality, we run each method for 5 different random initializations and report the average F IT . In addition, we evaluate model completion performance under different target ranks, R, from 10 to 60, and run 200 epochs.
First, we compare MULTIPAR's FIT with the baseline models on two datasets shown in Figure 4. MULTIPAR optimizes all prediction tasks and tensor factorization together. SPARTan and COPA first finish the tensor factorization, and then predict the downstream prediction tasks. SinglePAR optimizes the evaluated task and tensor factorization together. As Figure 4a and 4b shows, MULTIPAR outperforms all baseline methods on all datasets. In particular, MULTIPAR achieves a FIT score of 0.97 and 0.71 on MIMIC-EXTRACT and eICU respectively, a 13% and 40% relative improvement when compared to the best baseline model SinglePAR, which shows the strong tensor reconstruction ability of MULTIPAR, thanks to the "supervision" of the multiple predictive tasks. COPA performs better than SPARTan because it introduces various regularizations on the factor matrices, which can slightly improve the tensor reconstruction ability.
SinglePAR performs better than SPARTan and COPA on most of the ranks but is left behind COPA on large ranks. SinglePAR jointly optimizes prediction task and tensor factorization together. We can see that certain tasks benefit the tensor FIT while others may guide the tensor factorization into a suboptimal direction and degrade the tensor reconstruction quality. Although MULTIPAR model is supervised, thanks to the MTL, it can use all of the available outcomes across the different tasks to learn generalized representations of the data that are useful for tensor reconstruction.
Predictability analysis. A logistic regression model is trained on the patient importance membership matrix S k for static outcome prediction tasks and an LSTM model is trained on the temporal evolution matrix U k for dynamic outcome prediction task. LSTM is a variant of the recurrent neural network (RNN) that mitigates the gradient vanishing problem in traditional RNNs.  In the MIMIC-EXTRACT dataset, only the ventilation prediction task is a dynamic task, and all the tasks in the eICU dataset are static tasks. In order to illustrate the benefit of using the latent factors as features for a downstream prediction model, we also include an LSTM model trained using the original EHR data. The reason why we choose LSTM model is because the original EHR data contains different length patients' visit data, and LSTM model can handle the varying size input temporal data. The input to the LSTM model is an irregular tensor which contains k different patients, and each patient information X k consists of I K visits and J medical features. The output is the prediction label for the different patients and different visit for the dynamic task.
We evaluate the prediction accuracy as a function of the tensor factorization rank. As shown in Figures 5 and 6, MUL-TIPAR outperforms the other methods. In Figure 5, when the rank is 10, MULTIPAR outperforms the best baseline methods SinglePAR by 17%, 18%, 20% and 22% for each of the tasks respectively. This demonstrates MULTIPAR's strong generalization ability across multiple prediction tasks by leveraging the shared information between different tasks as well as the strong predictive power of the extracted phenotypes. Moreover, SinglePAR always outperforms COPA, SPARTan, and LSTM, which shows that the supervised learning framework can improve predictability. The Figure also illustrates the important role PARAFAC2 plays as the non-tensor based LSTM model performs the worst because it lacks the ability to filter out noise in the raw EHR.
Scalability analysis. Adding MTL on the PARAFAC2 framework can raise potential scalability issues on large datasets. Therefore, we evaluated the computational time of MULTIPAR compared with the other baseline methods using different data sizes and different feature sizes. We use two Titan RTX GPUs, each GPU has 24 GB of RAM, and train 50 epochs of both methods. Figure 7 shows the total training time. COPA, SPARTan, and SinglePAR shows linear scalability as the number of patients and features grows. Although MULTIPAR adds MTL, it still exhibits linear scalability similar to SinglePAR. While MTL adds some additional training time, it is not significantly more as the maximum added time is 8 minutes. Moreover, the experiment result is consistent with our computational complexity analysis, the per-iteration computational complexity is linear with respect to number of patients J and feature size K.
Interpretability analysis. Finally, we did an interpretability analysis of MULTIPAR on the MIMIC-EXTRACT dataset. We first illustrate the phenotypes discovered by MULTIPAR in Table 4. We set rank to 4, and use the V matrix to select the most important vital signs in each phenotype based on the weight. V matrix represents the membership of medical features in each phenotype, and the "weight" column in Table 4 is the weight in the V matrix. We then use the S matrix to find the patient subgroup of each phenotype, and calculate the average value of the vital signs shown in the "Average value" column. It is important to note that there is no post-processing in these extracted phenotypes.       A critical care expert reviewed and endorsed the presented phenotypes which suggest collective characteristics such as normal vital signs, abnormal renal and liver function, normal blood counts and serum electroytes, and abnormal vital signs.
The phenotypes discovered by the supervised single task model SinglePAR strongly overlap with each other as shown in Tables 5, 6, 7, and 8. Since we are incorporating in-hospital mortality prediction task in Table 5, most of the phenotypes discovered by SinglePAR are abnormal in vital signs. COPA discovered phenotypes shown in Table 9 contain more information compared to SinglePAR, which makes sense because a supervised model may guide the tensor factorization to a specific direction geared toward the task and cause information loss. However, MULTIPAR does not have information loss compared to COPA, it even provides a new phenotype (phenotype 2: abnormal in renal and liver function) which is not discovered by COPA. This verifies the benefit of MTL in MULTIPAR, which can leverage information from multiple tasks to avoid local optimum.
We then test MULTIPAR's ability to find meaningful subgroup temporal trajectories, which can help clinical care experts make precise prescriptions and treatments for specific subgroup of patients. We select the patients with the number of observations equal to 22 for visualization purposes. We select the four phenotypes for the temperature feature and systolic blood pressure feature, then use the S matrix to find the patient subgroup for each phenotype, and print the average value trajectory.
From Figure 8, we can see that the four patient subgroups (clusters) exhibit very different temporal trajectories in the temperature. Our clinical expert interpreted that the green line suggests a hyperthermic slow resolver patient subgroup which exhibits a slow decreasing trend as time increases, the red line suggests a hyperthermic fast resolver patient subgroup, which exhibits a fast decreasing trend as time increases, the dark blue line suggests a normothermic patient subgroup and the light blue line is a hypothermic patient subgroup. For the systolic blood pressure trajectory shown in Figure 8b, the green subgroup has high, increasing blood pressure, the red subgroup has high, decreasing blood pressure, the dark blue and light blue subgroups have consistently normal and low-normal blood pressure, respectively.

Conclusion
PARAFAC2 irregular tensor factorization has been studied to successfully extract meaningful medical concepts (phenotypes). Despite recent advancements, the predictability and interpretability is not satisfactory. We proposed MULTI-PAR: a supervised irregular tensor factorization with multi-task learning to jointly optimize the tensor factorization and multiple downstream prediction tasks. It is built on three major contributions: a supervised framework for PARAFAC2 tensor factorization and downstream prediction tasks; a new multi-task learning method combining tensor factorization and multiple prediction models; and a novel weight selection method for supervised multi-task optimization. We conducted extensive experiments to demonstrate that MULTIPAR can extract more meaningful phenotypes with higher predictability compared to state-of-the-art methods in a scalable way. In the future, we plan to incorporate lowrankness constraints for robustness against missing data, and incorporate more sophisticated regularization constraints to capture the complex and non-linear temporal relationships.