In classical classification models, such as logistic regression, a conditional class-generating probability distribution is modeled as a simple function of the observed features with unknown parameters to be trained. However, the raw observed features may be high-dimensional, and most of them might be uninformative and hard to interpret (e.g., pixel values of an image). Therefore, it would be desirable to extract more informative and interpretable low-dimensional features prior to the classification task. For instance, the multi-layer perceptron or deep neural networks (DNN) in general [8,9] use additional feature extraction layers prior to the logistic regression layer. This allows the model itself to learn the most effective (supervised) feature extraction mechanism and the association of the extracted features with class labels simultaneously.
Matrix factorization (MF) is a classical unsupervised feature extraction framework, which learns latent structures of complex datasets and is regularly applied in the analysis of text and images [13,32,44]. Various matrix factorization models such as singular value decomposition (SVD), principal component analysis (PCA), and nonnegative matrix factorization (NMF) provide fundamental tools for unsupervised feature extraction tasks [17,55,2,25]. Extensive research has been conducted to adapt matrix factorization models to perform classification tasks by supervising the matrix factorization process using additional class labels. Note that matrix factorization and classification are not necessarily aligned objectives, so some degree of trade-off is necessary when seeking to achieve both goals simultaneously. Supervised matrix factorization (SMF) provides systematic approaches for such multi-objective tasks. Our goal is to use SMF to learn low-rank latent factors that offer interpretable, data-reconstructive, and class-discriminative features, addressing challenges posed by high-dimensional data. The general framework of SMF was introduced in [33]. A similar SMF-type framework of discriminative K-SVD was proposed for face recognition [61]. A stochastic formulation of SMF was proposed in [30]. SMF has also found numerous applications in various other problem domains, including speech and emotion recognition [16], music genre classification [62], concurrent brain network inference [62], structure-aware clustering [59], and object recognition [27]. Recently, supervised variants of NMF, as well as PCA, were proposed in [5,26,47]. See also the survey work of [15] on SMF.
Various SMF-type models have been proposed in the past two decades. We divide them into two categories depending on whether the extracted low-dimensional feature or the feature extraction mechanism itself is supervised. We refer to them as feature-based and filter-based SMF, respectively. Feature-based SMF models include the classical ones by Mairal et al. (see, e.g., [33,30]) as well as the more recent model of convolutional matrix factorization by [22] for a contextual text recommendation system. Filter-based SMF models have been studied more recently in the supervised matrix factorization literature, most notably from supervised nonnegative matrix factorization [5,26] and supervised PCA [47].
In spite of vast literature on SMF, due to the high non-convexity of the associated optimization problem (see ( 4)), algorithms for SMF mostly lack rigorous convergence analysis and there has not been any algorithm that provably converges to a global minimizer of the objective at an exponential rate. We summarize our contributions below.
• We formulate a general class of SMF-type models (including both the feature-and the filter-based ones) with high-dimensional features as well as low-dimensional auxiliary features (see ( 4)). • We provide a novel framework that 'lifts' SMF as a low-rank matrix estimation problem in a combined factor space and propose an efficient algorithm that converges exponentially fast to a global minimizer of the objective with an arbitrary initialization (Theorem 3.5) . We numerically validate our theoretical results (see Fig. 2). • We theoretically compare the robustness of filter-based and feature-based SMF, establishing that the former is computationally more robust (see Theorem 3.5) while the latter is statistically more robust (see Theorem 4.1). • Applying our method to microarray datasets for cancer classification, we show that not only it is competitive against benchmark methods, but it is able to identify groups of genes including well-known cancer-associated genes (see Fig. 3).
Throughout this paper, we denote by R p the ambient space for data equipped with standard inner project h•, •i that induces the Euclidean norm k•k. We denote by {0, 1, . . . , } the space of class labels with + 1 classes. For a convex subset ⇥ in an Euclidean space, we denote ⇧ ⇥ the projection operator onto ⇥. For an integer r 1, we denote by ⇧ r the rank-r projection operator for matrices. For a matrix A = (a ij ) ij 2 R m⇥n , we denote its Frobenius, operator (2-), and supremum norm by
For each 1 i m and 1 j n, we denote A[i, :] and A[:, j] for the ith row and the jth column of A, respectively. For each integer n 1, I n denotes the n ⇥ n identity matrix. For square symmetric matrices A, B 2 R n⇥n , we denote A B if v T Av v T Bv for all unit vectors v 2 R n . For two matrices A and B, we denote [A, B] and [A k B] the matrices obtained by concatenating (stacking) them by horizontally and vertically, respectively, assuming matching dimensions.
Suppose we are given with n labeled signals (y i , x i , x 0 i ) for i = 1, . . . , n, where y i 2 {0, 1, . . . , } is the label, x i 2 R p is a high-dimensional feature of i, and x 0 i 2 R q is a low-dimensional auxiliary feature of i (p q). For a vivid context, think of x i as the X-ray image of a patient i and x 0 i denoting some biological measurements, such as gender, smoking status, and body mass index. When making predictions of y i , we use a suitable r (⌧ p) dimensional compression of the high-dimensional feature x i as well as the low-dimensional feature x 0 i as-is. We assume such compression is done by some matrix of (latent) factors W = [w 1 , . . . , w r ] 2 R p⇥r that is reconstructive in the sense that the observed signals x i can be reconstructed as (or approximated by) the linear transform of the 'atoms' w 1 , . . . , w r 2 R p for some suitable 'code' h i 2 R r . More concisely, X data = [x 1 , . . . , x n ] ⇡ WH, where H = [h 1 , . . . , h n ] 2 R r⇥n . In practice, we can choose r to be the approximate rank of data matrix X data (e.g., by finding the elbow of the scree plot).
Figure 1: Overall scheme of the proposed method for SMF-H. Now, we state our probabilistic modeling assumption. Fix parameters W 2 R p⇥r , h i 2 R r , 2 R r⇥ , and 2 R q⇥ . Let h : R ! [0, 1) be a score function (e.g., h(•) = exp(•) for multinomial logistic regression). We assume y i is a realization of a random variable whose conditional distribution is specified as [P (y i = 0 | x i , x 0 i ) , . . . , P (y i = | x i , x 0 i )] = g(a i ) := C[1, h(a i,1 ), . . . , h(a i, )],
(1) where C is the normalization constant and a i = (a i,1 , . . . , a i, ) 2 R is the activation for y i defined in two ways, depending on whether we use a 'feature-based' or 'filter-based' SMF model:
One may regard ( , ) as the 'multinomial regression coefficients' with input feature (h i , x 0 i ) or (W T x i , x 0 i ). In (2), we may regard the code h i (coming from x i ⇡ Wh i ) or the 'filtered signal' W T
x i as the r-dimensional compression of x i . Note that these two coincide if we have perfect factorization x i = Wh i and the factor matrix W are orthonormal, i.e., W T W = I r , but we do not necessarily make such an assumption.
There are some notable differences between SMF-H and SMF-W when predicting the unknown label of a test point. If we are given a test point (x test , x 0 test ), the predictive probabilities for its unknown label y test is given by (1) with activation a computed as in (2). This only involves straightforward matrix multiplications for SMF-W, which can also be viewed as a forward propagation in a multilayer perceptron [35] with W acting as the first layer weight matrix (hence named 'filter'). However, for SMF-H, one needs to solve additional optimization problems for testing. Namely, for every single test signal (x test , x 0 test ), its correct code representation h test needs to be learned by solving the following 'supervised sparse coding' problem (see [33]):
A more efficient heuristic testing method for SMF-H is by approximately computing h test by only minimizing the second term in (3).
In order to estimate the model parameters (W, H, , ) from observed training data (x i , y i ) for i = 1, . . . , n, we consider the following multi-objective optimization problem:
where
, and `(•) is the classification loss measured by the negative log-likelihood:
In ( 4), the tuning parameter ⇠ controls the trade-off between the two objectives of classification and matrix factorization. The above is a nonconvex problem involving four blocks of parameters that could have additional constraints (e.g., bounded norm). This problem entails many classical models as special cases. When ⇠ 1 so that effectively only the second term in ( 4) is being minimized with respect to W and H, it becomes the classical matrix factorization problem [31,28,29], where one seeks to find factor matrix W that can best reconstruct the feature vectors X data via the factorization X data ⇡ WH. In Figure 3, we will demonstrate that the best reconstructive factor matrix could be significantly different from the supervised factor matrix learned by SMF and may not be very effective for classification tasks.
The SMF training problem ( 4) is a nonconvex and possibly constrained optimization problem, generally with non-unique minimizers. Since it is difficult to solve exactly, approximate procedures such as block coordinate descent (BCD) (see, e.g., [57]) are often used. Such methods utilize the fact that the objective function in ( 4) is convex in each of the four (matrix) variables. Such an algorithm proceeds by iteratively optimizing for only one block while fixing the others (see [33,5,26,47]). However, convergence analysis or statistical estimation bounds of such algorithms are quite limited. Appealing to general convergence results for BCD methods (e.g., [18,58]), one can at most guarantee asymptotic convergence to the stationary points or polynomial convergence to Nash equilibria or of the objective (4), modulo carefully verifying the assumptions of these general results. We also remark that [30] provided a rigorous justification of the differentiability of a feature-based SMF model.
The main finding of our work is that the non-convexity of the SMF problem ( 4) is 'benign', in the sense that there exists an algorithm globally convergent to a global optimum at an exponential rate. We use a 'double-lifting' technique that converts the nonconvex SMF problem (4) into a low-rank factored estimation with a convex objective. This is reminiscent of the tight relation between a low-rank matrix estimation and a nonconvex factored estimation problem, which has been actively employed in a body of works in statistics and optimization [3,45,36,64,53,56,42,41,51]. Our exponentially convergent SMF algorithms are versions of low-rank projected gradient descent in the algorithm (44) that operate in the double-lifted space.
Our key idea to solve ( 4) is to transform it into a variant of the low-rank matrix estimation problem (6) and then use a Low-rank Projected Gradient Descent (LPGD) algorithm (7):
In (6), one seeks to minimize an objective f w.r.t. a paired matrix parameter Z = [✓, ] within a convex constraint set ⇥ and an additional rank constraint rank(✓) r. In (7), ⇧ r denotes applying rank-r projection on the first factor ✓ while keeping the same. Problem (6) encompasses a variety of significant problems, such as matrix regression [11,38] and matrix completion [48,37]. For these problems, algorithms of type (7) have been studied in [41,56].
To illustrate how the SMF problem (4) transforms into a low-rank matrix estimation (6), we consider a much simpler version of SMF-H. That is, instead of combining matrix factorization with multinomial logistic regression for multi-class classification problems, we combine it with linear regression. Thus, the response variable y in this discussion assumes all values in the real line. For additional simplicity, we assume there are no auxiliary features X aux . We seek to solve matrix factorization and linear regression problems simultaneously for data matrix X data 2 R p⇥n and response variable
This is a three-block optimization problem involving three factors W 2 R p⇥r , H 2 R r⇥n and 2 R r⇥1 , which is nonconvex and computationally challenging to solve exactly. Instead, consider reformulating this nonconvex problem as the following matrix factorization problem:
Indeed, we now seek to find two decoupled matrices (instead of three), one for T and W stacked vertically, and the other for H. A similar idea of matrix stacking was used in [61] for discriminative K-SVD. Proceeding one step further, another important observation we make is that it is also equivalent to finding a single matrix
⇥n of rank at most r that minimizes the function f in (9), which is convex (specifically, quadratic) in ✓: (See Fig. 1 Training).
For SMF-W, consider the following analogous linear regression model: min
where the right-hand side above is obtained by replacing H with W T X data in (9). Note that the objective function depends only on the product of the two matrices W and [ , H]. Then, we may further lift it as the low-rank matrix estimation problem by seeking a single matrix ✓ := [W , WH] 2 R p⇥(1+n) of rank at most r that solves (6) with f being the function in (10).
Motivated by the observation we made before, we rewrite SMF-H in (4) as
where
⇥n , and ⇥ is a convex subset of R (+p)⇥n ⇥ R q⇥ . We have added a L 2 -regularization term for A and with coefficient 0, which will play a crucial role in well-conditioning (11).
For solving (11), we propose to use the LGPD algorithm (7): We iterate gradient descent followed by projecting onto the convex constraint set ⇥ of the combined factor [✓, ] and then perform rank-r projection of the first factor ✓ = [A k B] via truncated SVD until convergence. Once we have a solution [✓ ? , ? ] to (11), we can use SVD of ✓ ? to obtain a solution to (4). Let ✓ ? = U⌃V T denote the SVD of ✓. Since rank(✓ ? ) r, ⌃ is an r ⇥ r diagonal matrix of singular values of ✓. Then U 2 R (+p)⇥r and V 2 R n⇥r are semi-orthonormal matrices, that is, U T U = V T V = I r . Then since ✓ ? = [( ? ) T k W ? ]H ? , we can take
V T and [( ? ) T k W ? ] = U⌃ 1/2 . We summarize this approach of solving (4) for SMF-H in Algorithm 1. Here, SVD r denotes rank-r truncated SVD and the projection operators ⇧ ⇥ and ⇧ r are defined in Subsection 1.1.
As for SMF-W, we can rewrite (4) with additional L 2 -regularizer for A = W and as
where +n) and ⇥ 2 R p⇥(+n) ⇥ R q⇥ is a convex set. Algorithm 1 for SMF-W follows similar reasoning as before with the reformulation above.
By using randomized truncated SVD for the efficient low-rank projection in Algorithm 1, the periteration complexity is O(pn min(n, p)), while that for the nonconvex algorithm is O((pr + q)n).
While the LPGD algorithm is in general more expensive per iteration than the nonconvex method, the iteration complexity is only O(log ✏ 1 ) thanks to the exponential convergence to the global optimum (will be discussed in Theorem 3.5). To our best knowledge, the nonconvex algorithm for SMF does not have any guarantee to converge to a global optimum, and the iteration complexity of the nonconvex SMF method to reach an ✏-stationary point is at best O(✏ 1 ) using standard analysis. Hence for ✏ small enough, Algorithm 1 achieves an ✏-accurate global optimum for SMF with a total computational cost comparable to a nonconvex SMF algorithm to achieve an ✏-stationary point.
Algorithm 1 Lifted PGD for SMF Input: X data 2 R p⇥n ; X 0 aux 2 R q⇥n (auxiliary features); Y label 2 {0, 1, . . . , } n Parameters: ⌧ > 0 (stepsize); N 2 N (iterations); r 1 (rank); 0 (L 2 -reg. param.)
We have discussed that one can cast the SMF problem (4) as the following 'factored estimation problem' min T,S, f (TS T , ). Note that such problems generally do not have a unique minimizer due to the 'rotation invariance'. Namely, let R be any r ⇥ r orthonormal (rotation) matrix (i.e.,
Hence the best one is to obtain parameters up to rotation that globally minimize the objective value. Our main result, Theorem 3.5, establishes that this can be achieved by Algorithm 1 at an exponential rate. First, we introduce the following technical assumptions (3.1-3.3). Assumption 3.1. (Bounded activation) The activation a 2 R defined in (2) assumes bounded norm, i.e., kak M for some constant M 2 (0, 1). Assumption 3.2. (Bounded eigenvalues of covariance matrix) Denote = [ 1 , . . . , n ] 2 R (p+q)⇥n , where
, where X aux = [x 0 1 , . . . , x 0 n ]. Then, there exist constants , + > 0 such that for all n 1,
Assumption 3.1 limits the norm of the activation a as an input for the classification model in ( 4) is bounded. This is standard in the literature (see, e.g., [36,60,23]) in order to uniformly bound the eigenvalues of the Hessian of the (multinomial) logistic regression model. Assumption 3.2 introduces uniform bounds on the eigenvalues of the covariance matrix. Assumption 3.3 introduces uniform bounds on the eigenvalues of the ⇥ observed information as well as the first derivative of the predictive probability distribution (see [10] and Sec. D in Appendix for more details). Under Assumption 3.1 and the multinomial logistic regression model h(•) = exp(•), one can derive Assumption 3.3 with a simple expression for the bounds ↵ ± , as discussed in the following remark. Remark 3.4 (Multinomial Logistic Classifier). Let `denote the negative log-likelihood function in (5), where we take the multinomial logistic model with the score function h(•) = exp(•). Denote ( ḣ1 , . . . , ḣ ) := r a `(y, a) and Ḧ(y, a) := r a r a T `(y, a). Then in this special case, we have ḣj (y, a) = g j (a) 1(y = j) and Ḧ(y, a) i,j = g i (a) (1(i = j) g j (a)) (See ( 28) and ( 30
Now, we state a special case of our first main result, specifically when the model is 'correctly specified', allowing the rank-r SMF model to effectively account for the observed data. This implies the existence of a 'low-rank stationary point' of F , as also demonstrated in [56]. However, we also handle the general case in Appendix (see Theorem D.1). Theorem 3.5. (Exponential convergence) Let Z t := [✓ t , t ] denote the iterates of Algorithm 1. Assume 3.1-3.3 hold. Let µ and L be as in (15), fix ⌧ 2 ( 1 2µ , 3 2L ), and let ⇢ := 2(1 ⌧ µ) 2 (0, 1). Suppose L/µ < 3 and let
In the statement above, we write kZk 2
Note that we may view the ratio L/µ that appears in Theorem 3.5 as the condition number of the SMF problem in (4), whereas the ratio L ⇤ /µ ⇤ for µ ⇤ := ↵ and L ⇤ := + ↵ + as the condition number for the multinomial classification problem. These two condition numbers are closely related. First, note that for any given µ ⇤ , L ⇤ and sample size n, we can always make L/µ < 3 by choosing sufficiently large ⇠ and so that Theorem 3.5 holds. However, using large L 2 -regularization parameter may perturb the original objective in (4) too much that the converged solution may not be close to the optimal solution. Hence, we may want to take as small as possible. Setting = 0 leads to
For SMF-W, if the multinomial classification problem is well-conditioned (L ⇤ /µ ⇤ < 3) and the ratio ⇠/n is in the above interval, then SMF-W enjoys exponential convergence in Theorem 3.5. However, the condition for SMF-H in ( 16) is violated, so L 2 -regularization is necessary for guaranteeing exponential convergence of SMF-H.
The proof of Theorem 3.5 involves two steps: (1) We establish a general exponential convergence result for the general LPGD algorithm (7) in Theorem C.2 in Appendix. (2) We compute the Hessian eigenvalues of the SMF objectives ( 11)-( 12) and apply the result to obtain Theorem 3.5. The proof contains two challenges: first, the low-rank projection in ( 7) is not non-expansive in general. To overcome this, we show that the iterates closely approximate certain 'auxiliary iterates' which exhibit exponential convergence towards the global optimum. Secondly, the second-order analysis is highly non-trivial since the SMF problem (4) has a total of four unknown matrix factors that are intertwined through the joint multi-class classification and matrix factorization tasks. See Appendix D for the details.
In this section, we formulate a generative model for SMF (4) and state statistical parameter estimation guarantee. Fix dimensions p q, and let n 1 be possibly growing sample size, and fix unknown true parameters B ? 2 R p⇥n , C ? 2 R q⇥n , ? 2 R q⇥ . In addition, fix A ? 2 R ⇥n for SMF-H and A ? 2 R p⇥ for SMF-W. Now suppose that class label, data, and auxiliary features are drawn i.i.d. according to the following joint distribution: 8 > > < > > :
where each " i (resp., " 0 i ) are p ⇥ 1 (resp., q ⇥ 1) vector of i.i.d. mean zero Gaussian entries with variance 2 (resp., ( 0 ) 2 ). We call the above the generative SMF model. In what follows, we will assume that the noise levels and 0 are known and focus on estimating the four-parameter matrices.
The (L 2 -regularized) negative log-likelihood of observing triples (y i , x i , x 0 i ) for i = 1, . . . , n is given as
, where c is a constant and F is as in (11) or ( 12) depending on the activation type with tuning parameter ⇠ =foot_0 2 2 . The L 2 regularizer in F can be understood as Gaussian prior for the parameters and interpreting the right-hand side above as the negative logarithm of the posterior distribution function (up to a constant) in a Bayesian framework. Note that the problem of estimating A and B are coupled due to the low-rank model assumption in (17), while the problem of estimating C is standard and separable, so it is not of our interest. The joint estimation problem for [A, B, ] is equivalent to the corresponding SMF problem (4) with tuning parameter ⇠ = (2 2 ) 1 . This and Theorem 3.5 motivate us to estimate the true parameters A ? , B ? , and ? by the output of Algorithm 1 with ⇠ = (2 2 ) 1 for O(log n) iterations. Now we give the second main result. Roughly speaking, it states that the estimated parameter Z t is within the true parameter Z ? = [A ? , B ? , ? ] within O(log n/ p n) with high probability, provided that the noise variance 2 is small enough and the SMF objective ( 11)-( 12) is well-conditioned. ). Then following holds with probability at least 1 1 n : For all t 1 and n 1,
We remark that Theorem 4.1 implies that SMF-H is statistically more robust than SMF-W. Namely, in order to have an arbitrarily accurate estimate with high probability, one needs to have µ p n log n. Combining with the expression in (15) and the well-conditioning assumption L/µ < 3, one needs to require ⇠ = ⌦(n), hence small noise variance 2 = O(1/n) for SMF-W. However, for SMF-H, this is guaranteed whenever 2 = o(1/( p n log n)) and ⇡ µ.
5 Simulation and Numerical Validation a Semi-synthetic MNIST data (" = 2) b Job postings data c LPGD vs.BCD MNIST Job postings SMF SMF SMF SMF SMF SMF SMF SMF We numerically verify Theorem 3.5 on a semi-synthetic dataset generated by using MNIST image dataset [24] (p = 28 2 = 784, q = 0, n = 500, = 1) and a text dataset named 'Real / Fake Job Posting Prediction' [1] (p = 2840, q = 72, n = 17880, = 1). Details about these datasets are in Sec. G in Appendix. 1 We used Algorithm 1 with rank r = 2 for MNIST and r = 20 for job postings datasets. For all experiments, = 2 and stepsize ⌧ = 0.01 were used.
We validate the theoretical exponential convergence results of our LPGD algorithm (Algorithm 1) in Figure 2. Note that the convexity and smoothness parameters µ and L in Theorem 3.5 are difficult to compute exactly. In practice, cross-validation of hyperparameters is usually employed. For ⇠ 2 {0.1, 1, 5, 10, 20} in Figure 2, we indeed observe exponential decay of training loss as dictated by our theoretical results for Algorithm 1. We also observe that the exponential rate of decay in training loss increases as ⇠ increases. According to Theorem 3.5, the contraction coefficient is ⇢ = (1 ⌧ µ), which decreases in ⇠ since µ increases in ⇠ (see (15)). The decay for large ⇠ 2 {10, 20} seems even superexponential. Furthermore, 2c shows that our LPGD algorithm converges significantly faster than BCD for training both SMF-H and SMF-W models at ⇠ 2 {5, 10} (other values of ⇠ omitted).
We apply the proposed methods to two datasets from the Curated Microarray Database (CuMiDa) [14]. CuMiDa provides well-preprocessed microarray data for various cancer types for various machine-learning approaches. One consists of 54,676 gene expressions from 51 subjects with binary labels indicating pancreatic cancer; Another we use has 35,982 gene expressions from 289 subjects with binary labels indicating breast cancer. The primary purpose of the analysis is to classify cancer patients solely based on their gene expression. We compare the accuracies of the proposed methods -SMF-W and SMF-H with a binary logistic classifier trained using Algorithm 1 -against the following benchmark algorithms: SMF-W and SMF-H trained using BCD; 1-dimensional sevenlayer Convolutional Neural Networks (CNN); three-layer Feed-Forward Neural Networks (FFNN); Naive Bayes (NB); Support Vector Machine (SVM); Random Forest (RF); Logistic Regression with Matrix Factorization by truncated SVD (MF-LR). For the last benchmark method MF-LR, we use rank-r SVD to factorize X train ⇡ U⌃V T and take W = U⌃ and H = V T . For testing, we use W T
X test as input to logistic regression for both filter and feature methods since We normalize gene expression for stable matrix factorization and interpretability of regression coefficients. We split each data into 50% of the training set and 50% of the test set and repeat the comparison procedure 5 times. A scree plot is used to determine the rank r. Other parameters are chosen through 5-fold cross-validation (⇠ 2 {0.1, 1, 10} and 2 {0.1, 1, 10}), and the algorithms are repeated in 1,000 iterations or until convergence. As can be seen in the table in Figure 3a, the proposed methods show the best performance for both types of cancers.
An important advantage of SMF methods is that they provide interpretable results in the form of 'supervised factors'. Each supervised factor consists of a latent factor and the associated regression coefficient. That is, once we train the SMF model (for = 1) and learn factor matrix W = [w 1 , . . . , w r ] 2 R p⇥r and vector of regression coefficients = [ 1 , . . . , r ] 2 R 1⇥r , each column w j of W describes a latent factor and the corresponding regression coefficient j tells us how w j is associated with class labels. The pairs (w j , j ), which form supervised latent factors, provide insights into how the trained SMF model perceives the classification task. See Fig. 1 for illustration.
In the context of microarray analysis for cancer research, each w j corresponds to a weighted group of genes (which we call a 'principal gene group') and j represents the strength of its association with cancer. SMF learns supervised gene groups (Fig. 3a, c) with significantly higher classification accuracy than the unsupervised gene groups (Fig. 3b, d). In Fig. 3a, c, both gene groups (consisting of p genes) have positive regression coefficients, so they are positively associated with the log odds of the predictive probability of the corresponding cancer. Remarkably, our method detected the well-known oncogene BRCA1 of breast cancer and other various genes (in Fig. 3e) that are known to be prognostic markers of breast/pancreatic cancer (see Human Protein Atlas [49]) in these groups of extreme coefficients (top five). The high classification accuracy suggests that the identified supervised principal gene groups may be associated with the occurrence of breast/pancreatic cancer.
We propose an exponentially convergent algorithm for nonconvex SMF training using new lifting techniques. Our analysis demonstrates strong convergence and estimation guarantee. We compare the robustness of filter-based and feature-based SMF, finding that the former is computationally more robust while the latter is statistically more robust. The algorithm's exponential convergence is numerically verified. In cancer classification using microarray data analysis, our algorithm successfully identifies discriminative gene groups for various cancers and shows potential for identifying important gene groups as protein complexes or pathways in biomedical research. Our analysis framework can be extended to more complex classification models, such as combining a feed-forward deep neural network with a matrix factorization objective. While our convergence analysis holds in certain parameter regimes, we discuss them in detail. We have tested our method and convergence analysis on various real-world datasets but recommend further numerical verification on a wider range of datasets.
We provide our implementation of Algorithm 1 in our code repository https://github.com/ljw9510/ SMF/tree/main.