{'fraction_non_alphanumeric': 0.047059180024847734, 'fraction_numerical': 0.01333292930517257, 'mean_word_length': 4.444976076555024, 'pattern_counts': {'":': 0, '<': 0, '<?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': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SUSTain: Scalable Unsupervised Scoring for Tensors and its Application to Phenotyping

This paper presents a new method, which we call SUSTain, that extends real-valued matrix and tensor factorizations to data where values are integers. Such data are common when the values correspond to event counts or ordinal measures. The conventional approach is to treat integer data as real, and then apply real-valued factorizations. However, doing so fails to preserve important characteristics of the original data, thereby making it hard to interpret the results. Instead, our approach extracts factor values from integer datasets as scores that are constrained to take values from a small integer set. These scores are easy to interpret: a score of zero indicates no feature contribution and higher scores indicate distinct levels of feature importance. At its core, SUSTain relies on: a) a problem partitioning into integer-constrained subproblems, so that they can be optimally solved in an efficient manner; and b) organizing the order of the subproblems’ solution, to promote reuse of shared intermediate results. We propose two variants, SUSTainM and SUSTainT, to handle both matrix and tensor inputs, respectively. We evaluate SUSTain against several state-of-the-art baselines on both synthetic and real Electronic Health Record (EHR) datasets. Comparing to those baselines, SUSTain shows either significantly better fit or orders of magnitude speedups that achieve a comparable fit (up to 425× faster). We apply SUSTain to EHR datasets to extract patient phenotypes (i.e., clinically meaningful patient clusters). Furthermore, 87% of them were validated as clinically meaningful phenotypes related to heart failure by a cardiologist.

Factorization methods have been successfully used for EHR-based unsupervised phenotyping [21,22,24,41,42,50]. In many of those settings, the problem can be formulated via Nonnegative Matrix Factorization (NMF) [35] e.g., minimizing the squared Frobenius norm of the error: min X − U V T F 2 | U ≥ 0, V ≥ 0 (1) X ∈ ℤ + M × N is a non-negative integer input matrix whose X(i, j) cell reflects the event counts for the i-th (out of M) patient with respect to the j-th (out of N) features. Given an input number R of desired phenotypes, the matrix U ∈ ℝ M × R corresponds to a membership matrix of the patients with respect to the R phenotypes. And the matrix V ∈ ℝ N × R provides the phenotypes' definition: the non-zero elements of the r-th column V (:, r) reveal the potentially relevant features to the r-th phenotype.
Interpreting those factors is crucial in order to determine whether and to what extent a patient exhibits a phenotype, as well as which set of candidate features should be considered to compose each r-th phenotype so that it is clinically meaningful. However, this can be challenging if the resulting factors contain arbitrary (non-negative) real values. Real-valued factors distort the count nature of input data; thus, identifying cases and controls based on counts of relevant medical features [13,34] is no longer possible. Also, the possible ranges and relative differences of elements in real-valued factors is arbitrary, thus impeding the practitioner's assessment of their relative importance. In practice, ad hoc heuristics have been introduced with limited success: a) hard thresholding to the ranked list of factor elements, which is usually arbitrary and leads to poor model fit; b) the factor values are hidden altogether and only the elements' ranking is preserved, which omits valuable information regarding the individual elements' actual importance.

Contributions:
To tackle these challenges, we propose Scalable Unsupervised Scoring for Tensors (SUSTain), a framework extracting the factor values as scores, constrained to a small integer set. SUSTain offers a straightforward interpretation protocol: a score of zero indicates no feature contribution and higher scores indicate distinct levels of feature importance.
Our methodology relies on identifying a problem partitioning into integer-constrained subproblems so that each one of them can be solved optimally in an efficient manner; at the same time, their solution order is organized so as to promote re-use of shared intermediate results. SUSTain can handle both matrix and tensor inputs, through SUSTain M and SUSTain T methods, which we formulate in Sections 3.1 and 3.2 respectively.
SUSTain yields faster and more scalable approaches than baselines achieving comparable fit, as evaluated on both synthetic (publicly-available) and real healthcare datasets. For example, as shown in Figure 1, SUSTain M achieves the same level of accuracy as the most accurate baseline up to 425× faster. SUSTain T can handle large-scale tensor inputs for which the most accurate baseline fails and scales linearly with the number of patients.
SUSTain's interpretation protocol is particularly meaningful for unsupervised phenotyping: it is easily understood by medical experts who are used to simple and concise, scoring-based descriptions of a patient's clinical status (e.g., risk scores 1 ). While recent work derives risk scores for predictive modeling (supervised learning) [48], our application of SUSTain extracts scores for unsupervised phenotyping based on unlabeled EHR data. In Table 1, we provide a representative phenotype extracted through our method, as part of a case study we performed on phenotyping heart failure patients. The meaningfulness of the phenotype candidates extracted through this case study was confirmed by a cardiologist, who annotated 87% of them as clinically meaningful phenotypes related to heart failure. We summarize our contributions as: • Scalable unsupervised scoring: We propose SUSTain, a fast and scalable approach decomposing integer multi-aspect data into integer scores, preserving the original integer characteristics.
• SUSTain can handle matrix and tensor input: We present SUSTain for both matrix (Section 3.1) and tensor (Section 3.2) inputs, through SUSTain M and SUSTain T methods, respectively.

• Evaluation on various datasets:
We evaluate both the matrix and tensor versions on both synthetic (publicly-available) and real healthcare datasets.
• Phenotyping heart failure patients: The interpretability of the extracted scoring-based phenotypes was confirmed by a cardiologist, who annotated 87% of them as clinically meaningful.

BACKGROUND
In Table 2 we summarize the notations used throughout the paper. Let x 0 ∈ ℝ n . The euclidean projection of x 0 to a set C ⊆ ℝ n is defined as Π C x 0 = argmin x − x 0 2 2 | x ∈ C ; thus, it is the problem of determining the vector x* among all x ∈ C which is the closest to x 0 w.r.t. the Euclidean distance [5]. A matrix X is called rank-1 if it can be expressed as the outer product of 2 non-zero vectors: X = x ∘ y. The Khatri-Rao Product (KRP) [46] is the "matching column-wise" Kronecker product: for two matrices U ∈ ℝ M × R , V ∈ ℝ N × R their KRP is as follows: A tensor is a multi-dimensional array. The tensor's order denotes the number of its dimensions, also known as ways or modes (e.g., matrices are 2-order tensors). A d-order tensor X is called rank-1 if it can be expressed as the outer product of d non-zero vectors: X = a 1 ∘ a 2 ∘ ⋯ ∘ a d . A fiber is a vector extracted from a tensor by fixing all modes but one. For example, a matrix column is a mode-1 fiber. A slice is a matrix extracted from a tensor by fixing all modes but two. Matricization, also called reshaping or unfolding, logically reorganizes tensors into other forms without changing the values themselves. The mode-n matricization of a d-order tensor X ∈ ℝ I 1 × I 2 × ⋯ × I d is denoted by X (n) ∈ ℝ I n × I 1 I 2 …I n − 1 I n + 1 …I d and arranges the mode-n fibers of the tensor as columns of the resulting matrix. The Matricized-Tensor Times Khatri-Rao Product (MTTKRP) [3] w.r.t. mode-n is the matrix multiplication X (n) A ⊙ ( − n) , where A ⊙ ( − n) corresponds to the Khatri-Rao product of all the modes except the n-th. MTTKRP is the bottleneck operation in many sparse tensor algorithms.

THE SUSTAIN FRAMEWORK
First we present SUSTain for matrix input. Then, we describe how SUSTain can be extended for general high-order tensor input. Finally, we provide our interpretation protocol of SUSTain for unsupervised phenotyping.

SUSTain for matrix input
Model: For an integer input matrix X ∈ ℤ + M × N and a certain target rank R, the problem can be defined as: where ℤ τ = 0, 1, …, τ is the set of nonnegative integers up to τ, ℤ + = 1, 2, …, ∞ is the set of positive integers and Λ is a diagonal R-by-R matrix. The above problem can be also normalization as in the corresponding real-valued models (e.g., NMF [25,35]), due to the integer constraints. Finally, note that the integer set ℤ τ can easily vary for different factor matrices and even allow negative integers; this can also happen for the input matrix X. The formulation in Problem (2) favors simplicity of presentation and matches the need of phenotyping applications.

Fitting Algorithm:
We employ an alternating updating scheme to tackle the non-convex optimization Problem (2). Our scheme leads to optimal solutions to each one of the subproblems in an efficient manner, while organizing the order of updates so as to promote reuse of already computed intermediate results.
We follow the intuition behind the Hierarchical Alternating Least Squares (HALS) framework, which enables isolating and solving for each k-th rank-1 component separately.
where R k corresponds to the "residual matrix" and is considered fixed when solving for the k-th rank-1 component. The objective can be written as [30]: We set: If U (:, k) T R k V (:, k) > 0 then the minimum value of J for λ(k) ∈ ℤ + is obtained at max 1, round λ k * where round() rounds to the nearest integer. If U (:, k) T R k V (:, k) ≤ 0, then the minimum objective value for λ(k) ∈ ℤ + is attained at λ(k) = 1. Combining these two cases, the optimal λ(k) ∈ ℤ + is given by: In practice, R k may be large (M × N) and dense, even if the input is sparse (as happens in our main motivating application); thus its explicit materialization should be avoided [16,27]. Expanding the above expression gives: Next, solving Problem To solve the above, we apply the Optimal Scaling Lemma [8] for the integer constraint. This Lemma states that for any set C of constraints imposed on b, it holds that: is the unconstrained solution to the above problem. This means that the optimal solution of the constrained problem is simply the projection of the unconstrained solution onto the constraint set C. Thus, the optimal solution of Problem (6) is: Since ℤ τ N is the Cartesian product of subsets of the real line, i.e., ℤ τ thus project each scalar coordinate individually. For a real-valued scalar α, projecting onto ℤ τ gives [47]: Finally, expanding R k in Expression (7), combining with (8), (9) and setting: b V (: , k) + X T U : , k − V Λ U T U : , k U T U k, k λ(k) (10) gives the optimal solution for Problem (6): where min(), max(), round() are taken element-wise.
Having derived the updates for λ(k), V (:, k), in Relations (5) and (11)  To exploit that, we choose to successively update λ(k) and V(:, k) during the same iteration and iterate ∀k ∈ {1, … , R}. As a result of the proposed update order, the only non-negligible additional operation in order to compute both λ(k) andV (:, k) is to re-compute t ≔ V Λ U T U] :, k after having updated λ(k).
Re-computing t can be further optimized by observing that only the contribution of the k-th component t k ≔ V (:, k) * λ(k) [U T U] k, k has to be adjusted. Thus, we can store t and t k , compute λ′(k), and then adjust t as: Updating U (:, k) can be executed in symmetric fashion to V (:, k). In Algorithm 1, we present our main procedure to update both the factor matrices U and V and λ values in an alternating fashion. In Algorithm 2, we provide the definition of SUSTain_Update_Factor which updates a single factor (denoted as F) and the vector λ.

SUSTain M
Require: X ∈ ℝ M × N , target rank R and upper bound τ

9: end for
Computational Complexity: The asymptotic cost of executing Algorithm 2 is 2R 2 I flops (i.e., floating-point operations), ∀R > 5. This step costs 2R 2 N when updating V and 2R 2 M when updating U. In Algorithm 1, assuming the input X is sparse, the cost of each one of X V and X T U is 2 nnz (X) R flops. Also, computing V T V and U T U cost 2R 2 N and 2R 2 M flops respectively. Thus, the total cost is: 4R (nnz(X) + (M + N)R) flops.

SUSTain for tensor input
Model: For a tensor X ∈ ℝ I 1 × I 2 × …I d of order d and a certain target rank R, the problem can be defined as where n = {1, … , d}, ℤ τ = 0, 1, …, τ is the set of nonnegative integers up to τ and ℤ + = 1, 2, …, ∞ is the set of positive integers. Our model is an extension of SUSTain M presented in Section 3.1 for high-order tensors. It can be viewed as a constrained version of the CP tensor model [9,18,20].
Fitting Algorithm: Similarly to the matrix case, we set: Thus, Problem (12) becomes: We matricize the above expression w.r.t. mode-n and utilize the fact that the mode-n matricization of a rank-1 tensor b 1 We set A ⊙ ( − n) ≔ A (d) ⊙ ⋯ ⊙ A (n + 1) ⊙ A (n − 1) ⊙ ⋯ ⊙ A (1) as the Khatri-Rao Product of all the factor matrices except the n-th and as the Hadamard product of the Gram matrices of all the factor matrices except the n-th.
Then, Objective (14) becomes Perros et al. Page 9 Solving the above for λ(k) can be handled equivalently to the corresponding matrix case (Relation (4)). Thus, the optimal solution for λ(k) ∈ ℤ + is: By exploiting that [29]: and expanding: (18) where M (n) (:, k) is the Matricized-Tensor Times Khatri-Rao Product (MTTKRP) [3] operation w.r.t. mode n, we get the optimal solution for λ(k) ∈ ℤ + as: Next, we transpose the Objective (16) and solving for A (n) (:, k) can be handled as in the matrix case (Relation (7)) through the(Optimal Scaling Lemma [8]. Thus, the optimal A (n) (: , k) ∈ ℤ τ I n is given by: Finally, combining Equation (18) into the above gives: Note the direct correspondence of the above formulations for λ(k), A (n) (:, k) with the core update Algorithm 2 we used for the matrix case. If we set F ← A (n) , M ← M (n) , C ← C (−n) then we can simply use Algorithm 2 to update a single factor A (n) and the λ values. Also, we can exploit the development of existing scalable software libraries computing the bottleneck MTTKRP kernel for sparse data efficiently [3]. In Algorithm 3, we summarize the operations of our methodology for tensor input.
Computational Complexity: Updating the n-th mode in Algorithm 3 requires: 3 R nnz(X) flops to compute the MTTKRP using state-of-the-art libraries for sparse tensors [3], 2 R 2 I n flops to compute A (n) T A n and (d − 1) R 2 flops to update C (−n) as in Equation

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript (15). As discussed in the matrix case, the dominant cost of Algorithm 2 is 2 R 2 I n flops.

Interpretation for phenotyping
Given the EHRs of a certain cohort, we form a patient-by-diagnoses matrix X, whose X(i, j) cell is the number of encounters of patient i where encounter diagnosis j was recorded. In that case, the patient membership vector U (i, :) of SUSTain M provides the distinct levels of frequency of each one of the R phenotypes throughout the medical history of the i-th patient. Likewise, each column V (:, r) indicates the frequency levels of each medical feature w.r.t. the r-th phenotype. Table 1 summarizes a phenotype example that accounts for the largest share of heart failure patients. Finally, due to the integer box (i.e., {0, … , τ}) constraints employed on the factor matrices, we can interpret the integer λ(r) values as scaling up the input encounter counts for the r-th phenotype. Thus, phenotypes with higher λ(r) values are expected to describe more persistent medical conditions, with higher number of associated encounters.
The above interpretation can be extended to the tensor case. Consider a tensor X whose X i, j, k cell defines the count of encounters of patient i where medication k was ordered for the patient with diagnosis j as the order indication. Factorizing this tensor using SUSTain T yields a patient factor A (1) which can be interpreted similarly to the U factor in the matrix case. Also, the factor matrices A (2) , A (3) corresponding to diagnosis and medication or procedure phenotypes can be interpreted similarly to the V factor in the matrix case. The same applies to the λ(r) model values.
Sutter: This dataset corresponds to EHRs from Sutter Palo Alto Medical Foundation (PAMF) Clinics. The patients are 50 to 80 years old adults chosen for a heart failure study [11]. To form a patient-by-diagnosis matrix input, we extracted the number of encounter records with a specific diagnosis for each patient. To form a patient-diagnosis-medication tensor input, we used the medication orders, reflecting the ordered medications and the indicated diagnosis. We adopt standard medical concept groupers to group the available ICD-9 diagnosis codes [45] into Clinical Classification Software (CCS) [1] diagnostic categories (level 4). We also group the normalized drug names (i.e., combining all branded names and the generic name for a medication) based on unique therapeutic subclasses using the Anatomical Therapeutic Chemical Classification System [39].

CMS: We used a publicly-available CMS Linkable 2008-2010 Medicare Data
Entrepreneurs' Synthetic Public Use File (DE-SynPUF) 2 that contains three years of claim records synthesized (i.e., to protect privacy) from 5% of the 2008 Medicare population. CMS creates twenty 5% subsamples of the claims data. We used the carrier claims data available from DE-SynPUF for the patients belonging to Samples 1 & 2. We increase the number of samples (i.e., number of patients) considered for the experiments related to assessing scalability. We used the diagnostic code information to build the input matrix and the diagnoses and procedures recorded to build the input tensor. In particular, we group the available ICD-9 diagnosis codes [45] into CCS [1] diagnostic categories (level 4) and use the CCS flat code grouper [1] to transform the CPT procedure codes available into procedure categories.

4.1.2
Baselines.-Below, we describe our efforts to design competitive baseline methods producing the target models in Problems 2 and 12, for the matrix and the tensor cases respectively.
Round: This baseline rounds the factor matrices from nonnegative matrix/tensor factorization. In the matrix case, we used the implementation of Nonnegative Matrix Factorization (NMF) [27,28] and projected all the entries of the resulting factor matrices to ℤ τ . We also set λ to an all-ones vector as NMF typically does not have the diagonal matrix Λ. A typical issue of naively rounding NMF solutions is that values that are lower than 0.5 are rounded to 0, so a potentially large part of model information can be lost.
In the tensor case, we used the CP-ALS algorithm as in the Tensor Toolbox [4], adjusted to impose non-negativity constraints [28] on the factor matrices. Also, in contrast to the NMF case, CP-ALS produces a λ vector of nonnegative real values. In order to alleviate the effect of zeroing out values less than 0.5 we compute the cube root of the λ vector element-wise and form a vector λ . Then, we absorb this scaling in the factor matrices by multiplying A (n) diag(λ ), ∀n = 1, …, d where d is the input tensor's order. Finally, we set λ to an allones vector and project all the entries of the resulting factor matrices to ℤ τ .

Scale-and-round:
We design a more sophisticated scale-and-rounding heuristic which scales the factor matrices of the real-valued solutions before performing the rounding. This step further alleviates the problem of zeroing out values less than 0.5.
In the tensor case, we absorb the scaling of the λ output of the real-valued solution into the factor matrices as in "Round", and extend the Scale-and-round matrix approach accordingly.

AILS: Alternating Integer Least Squares approach:
We used the Integer Least Squares (ILS) with box constraints approach which is proposed in [6,49]. This approach was recently unified within an Integer Matrix Factorization framework [15]. We exploit the redundancy among ILS problems targeting the same factor matrix, so that the QR factorization in the reduction phase is only computed once. Note that solving general ILS problems is NP-hard [15], which is reflected in the runtime of this method in the experiments. We enabled the extraction of the integer λ values through an ILS by noticing that vectorizing the original problem as can be transformed to [7]: min (V ⊙ U)λ − vec(X) F 2 | λ ∈ ℤ + R which gives the ILS to solve for. Note that we attempted to extend this approach for tensor input; however, the materialization of the Khatri-Rao product of all the factor matrices failed due to out of memory problems even for the smallest target rank for both of the datasets used. To illustrate the magnitude of this issue, the size needed for the Khatri-Rao product of all factor matrices for Sutter data and R = 5 is: 248347 * 552 * 555 * 5 * 8 bytes ≈ 3 Terabytes. Regarding the accuracy-time trade-off evaluation, we initialize with several schemes and for each method we choose the one providing the highest fit. The schemes are the following: a) round heuristic, b) scale-and-round heuristic, c) random: random initialization with integers within the required range and λ set to all-ones vector, d) random & sampling: random initialization of the patients factor and sampling from the input data to populate the rest of the factors. In the matrix case, we initialize each j-th column of V by random sampling of input patient vectors and scaling them to lie on ℤ τ if needed. In the tensor case, for each sampled slice X i, : , : , we populate each j-th component of A (2) , A (3) by sampling the row and column of X i, : , : with the maximum sum. Note that when we measure execution time for each approach, we do take into account the time spent for its initialization.

Evaluation metrics.-We
In the scalability evaluation, we initialize each method with the random & sampling scheme (d) described above; this provided better starting points than using pure random initialization. For this experiment, we ignore the initialization time, since we want to focus on the methods' scalability behavior.

Implementation details.-
We used MatlabR2017b for our implementations, along with functionalities for sparse tensors from the "Tensor Toolbox" [4] and for nonnegative matrix factorization from the "nonnegfac-matlab" [27] toolbox. The ILS solver we use for the AILS baseline is included in the state-of-the-art MILES software [10].
The zero-lock problem refers to the case when a single column is zeroed out, thus zeroing out an entire rank-1 component of the solution. To avoid that in our scheme, we add the smallest perturbation possible (+1) to a randomly-chosen coordinate of the vector zeroed out.
In both SUSTain and AILS, we break the iterations when the successive difference of the objective drops below 1e − 4. Finally, the parameter τ is set to 5 driven by discussions with medical experts and similarity to many medical scoring systems.

4.1.6
Hardware.-We conducted our experiments on a server running Ubuntu 14.04 with 1TB of RAM and four Intel E5-4620 v4 CPU's with a maximum clock frequency of 2.10GHz. Each of the processors contains 10 cores with 2 threads each.

Matrix case experiments
Accuracy-Time trade-off: In Figure 1, we showcase the accuracy-time trade-off regarding the Sutter PAMF dataset. SUSTain M is at least 60× faster (R = 40) than the most accurate baseline (AILS). For R = 5, SUSTain M achieves 425× speedup over AILS: as compared to the ≈ 22 minutes spent by AILS, our approach executes in ≈ 3 seconds for the same level of accuracy. Even for R = {10, 20} SUSTain M achieves 98× and 110× faster computations than AILS. At the same time, SUSTain M achieves up to 16% higher fit than the scale-and-round heuristic, operating on comparable running times. Note that for R = 5, our approach is even faster (and more accurate) than the scale-and-round baseline as well, since initializing with random factors provided a better final fit than initializing with the scale-and-round result. We also remark that the naive round heuristic achieves a fit of zero, which is a by-product of zeroing out the majority of the model factor elements.
In Figure 2, we provide the results of the same experiment regarding the CMS dataset. For the same level of accuracy, SUSTain M is at least an order of magnitude faster than AILS, and up to 38× faster for R = 20. It also achieves up to 14% higher fit over the scale-androunding heuristic for comparable execution time. Table 4, for fixed R = 10, we measure a single iteration's time for increasing subsets of CMS patients. The NMF execution time is considered for Round and Scale-and-round heuristics, since their post-processing cost is negligible. SUSTain M can execute very fast (a single iteration in ≈ 3 seconds) even for ≈ 985 thousand patients.

Tensor case experiments
Accuracy-Time trade-off: In Figure 3, we provide the fit-time trade-off for varying target rank of our input tensor datasets. As discussed in Section 4.1.2, the extension of AILS approach to tensors cannot scale for any dataset or target rank considered. Overall, SUSTain T achieves up to 9% and 12% increase in fit over the scale-and-round heuristic w.r.t. the Sutter PAMF and CMS datasets respectively. Note that the fit of the scale-andround approach decreases for successively increasing target rank values (e.g., transitioning from R = 20 to R = 40 for CMS data). This indicates that heuristic approaches which simply post-process real-valued solutions may not fully exploit the available target rank.
Scaling for larger number of patients: In Table 5, we report the time spent for one iteration of increasingly larger subset of patients considered from the CMS data, with fixed target rank (R = 10). The time measured for the heuristic approaches corresponds to the execution time of CP-ALS, since the post-processing cost is negligible. We observe that SUSTain T achieves linear scale-up w.r.t. increasing number of patients. We also remark that the dominant cost in both SUSTain T and the CP-ALS is the MTTKRP computation, which explains the comparable running time.

Case study on Phenotyping HF patients
Cardiovascular disease (CVD) is the leading cause of death worldwide and heart failure (HF) is a dominant cause of morbidity and mortality. HF is traditionally characterized by reduced ejection fraction (HFrEF) and preserved ejection fraction (HFpEF). But, recent evidence suggests that HF is more heterogeneous than is reflected by ejection fraction. We used SUSTain to explore this heterogeneity in an incident HF cohort.

Cohort and data selection:
We select only the HF case patients from the Sutter PAMF dataset. For each incident HF case, we extracted data in the 12-months before and the 12months after the initial HF diagnosis date, which resulted in 70, 531 clinical encounters. We used all the data modalities available, i.e., medication orders and indications and encounter diagnoses. The size of the resulting (patient-by-diagnosis-by-medication) tensor is 3, 497 × 396 × 367; the tensor contains a total of 92, 662 non-zero elements.
Choosing the number of phenotypes: We use the stability-driven criterion introduced in [51]. The intuition behind this criterion is in promoting a target rank for which several runs with different initial points return reproducible factors. We choose the diagnosis factor matrix as the factor under assessment. Let D 1 and D 2 be the diagnosis factor matrix for 2 different runs with the same target rank. Then, the cross-correlation matrix C ∈ ℝ R × R is computed between the columns of D 1 , D 2 and the dissimilarity between them is computed as [51]: Note that when D 1 can be transformed to D 2 by column permutation, then diss(D 1 , D 2 ) = 0. If B is the number of repetitions for each target rank, then the following relation computes the average dissimilarity over B(B − 1)/2 pairs of resulting factors: We used the "staNMF" toolbox 3 to compute the above score for each target rank on the range {5, … , 20}. The input to SUSTain T were B = 20 initial points of the round heuristic. R = 15 phenotypes were selected based on the above criterion. For the target rank chosen, we pick the solution yielding the highest fit.

SUSTain provides concise and accurate solutions:
We observed that besides preserving the input data properties and providing a natural interpretation for medical experts, SUSTain implicitly imposes sparse factors. To assess the factors' conciseness, we compare their fit with the achieved fit of the real-valued model (NN CP-ALS), which is post-processed to achieve factor sparsity (as would be done by a practitioner). For each of the feature factors (diagnosis and medication) of the real-valued model, we only consider the top-k elements for each column (i.e., most important elements of each phenotype). For the patient factor, we consider the top-k elements for each row (i.e., most important phenotypes for each patient). In each case, the value of k is chosen so that the sparsity level is close to the one achieved by SUSTain T . We provide the results in Table 6, where we notice that for the same level of sparsity, SUSTain T achieves ≈ 8.6% increase in fit. Thus, the integer factors of SUSTain T decompose the input more accurately for the same level of sparsity than the real-valued counterpart.
Phenotype discovery: In Table 1 and Table 7, we provide representative phenotypes extracted through our method. A subset of annotations provided by the cardiologist are as follows: hyperlipidemia (the one in Table 1), HF with reduced LVEF (HFrEF), hypertension (HTN), HTN which is more difficult to control, persistent and chronic atrial fibrillation, depression, diabetes, comorbidities of aging, prior pulmonary embolism. Overall, 13 out of 15 phenotype candidates were annotated as clinically meaningful phenotypes related to heart failure.

RELATED WORK Discrete factorization-based approaches:
Dong et al. [15] proposed an Integer Matrix Factorization framework via solving Integer Least Squares subproblems. As we experimentally evaluated, this approach is orders of magnitude slower than SUSTain while achieving the same level of accuracy. Kolda and O'Leary [30,40] proposed a Semidiscrete Matrix Decomposition into factors containing ternary values ({−1, 0, 1}). Despite its demonstrated success for compression purposes, a direct application of this approach would introduce negative values into the factors, thus hurting interpretability for nonnegative input. Finally, several prior works target binary factorization (e.g., [32, 33, 36-38, 44, 53]). In contrast to strictly binary factors, SUSTain captures the quantity embedded in the input data, which reveals important information (e.g., relative phenotype prevalence and associated feature frequencies).

HALS fitting algorithms:
Our fitting algorithms follow the intuition of Hierarchical Alternating Least Squares (HALS) framework [12] (aka rank-one residue iteration [23]), which enables formulating the solution for each k-th rank-1 component separately. However, plain HALS does not tackle the challenges involved with either imposing integer constraints or solving for the vector λ.

CONCLUSIONS
The accuracy and scalability of SUSTain on "native" integer data derives from two key insights. One is expected: just rounding or applying related transformations to real-valued solutions is inherently limited. The second may be more surprising: while discrete constraints might appear to make the problem more challenging, in fact, a careful organization of the problem into subparts can mitigate that complexity. In our case, we identify a problem partitioning of integer-constrained subproblems that leads to an optimal and efficient solution; and, we also define the order of alternating updates so as to enable reuse of shared intermediate results. Consequently, SUSTain outperforms several baselines on both synthetic (publicly-available) and real EHR data, showing either a better fit or orders-of-magnitude speedups at a comparable fit.
Moving forward, there are many other sources of integer values in real-world data. These include, for instance, ordinal values. Thus, whereas this paper targets event counts, extensions for other cases is a ripe target for future work.
Lastly, to enable reproducibility of our work, we open-source our implementations and make them publicly available at: https://github.com/kperros/SUSTain. for the Sutter and the CMS tensor input. SUSTain T achieves up to 9% and 12% higher fit respectively over scale-and-rounding heuristics.