Exponentially Convergent Algorithms for Supervised Matrix Factorization
Authors: Joowon Lee
Hanbaek Lyu
Weixin Yao
Venue: NeurIPS 2023 poster (2023)
Keywords: Supervised matrix factorization
multi-objective optimization
global convergence
linear convergence
statistical estimation
TL;DR: Supervised matrix factorization, a multi-objective non-convex problem, is lifted to a low-rank matrix estimation problem and is solved by using projected low-rank gradient descent algorithms. Exponential convergence to global optima is established.

Abstract: Supervised matrix factorization (SMF) is a classical machine learning method that simultaneously seeks feature extraction and classification tasks, which are not necessarily a priori aligned objectives. 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. Training SMF model involves solving a nonconvex and possibly constrained optimization with at least three blocks of parameters. Known algorithms are either heuristic or provide weak convergence guarantees for special cases. In this paper, 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 provably converges exponentially fast to a global minimizer of the objective with arbitrary initialization under mild assumptions. Our framework applies to a wide range of SMF-type problems for multi-class classification with auxiliary features. To showcase an application, we demonstrate that our algorithm successfully identified well-known cancer-associated gene groups for various cancers.

Exponentially Convergent Algorithms for Supervised Matrix Factorization
Joowon Lee Department of Statistics
Hanbaek Lyu Department of Mathematics University of Wisconsin Madison, WI, USA hlyu@math. wisc.edu
University of Wisconsin Madison, WI, USA jlee2256dvisc. edu
Weixin Yao Department of Statistics
University of California, Riverside, CA, USA weixinyducr.c .ddu edu
Abstract
Supervised matrix factorization (SMF) is a classical machine learning method that simultancously seeks feature extraction and classification tasks, which are not necessarily a priori aligned objectives. Our goal to use SMF to learn low-rank latent factors that offer interpretable, data-reconstruuctvv,, and clas--disrrminnaiiee features, addressing challenges posed by high-dimensional data. Training SMF model involves solving nonconvex and possibly constrained optimization with at least three blocks of parameters. Known algorithms are either heuristic or provide weak convergence guarantees for special cases. In this paper, provide a novel framework that 'lifts' SMF as low-rank matrix estimation problem in a combined factor space and propose efficient algorithm that provably converges exponentially fast to a global minimizer of the objective with arbitrary initialization under mild assumptions. Our framework applies to wide range of SMF-type problems for multi-class classification with auxiliary features. To showcase an application, demonstrate that algorithm successfully identified well-known cancer-associated gene groups for various cancers.
Introduction
classical classification models, such as logistic regression, a conditional class- geass-ggnnrainng probability distribution is modeled as simple function of the observed features with unknown parameters to be trained. However, the raw observed features may be nigh-dimensional, and most of them might be uninformative and hard interpret (e.g., pixel values of image). Therefore, it would be desirable to extract more informative and interpretable oow-dimensional features prior to the classification task. For instance, the multi-layer perceptron deep neural networks (DNN) in general [8, 9] use additional feature extraction layers prior to the logistic regression layer. This allows the model itself learn the most effective (supervised) feature extraction mechanism and the association of the extracted features with class labels simultaneously.
Matrix factorization (MF) is classical unsupervised feature extraction framework, which learns latent structures of complex datasets and is regularly applied in the analysis text and images [13, 32, 44]. Various matrix factorization models such singular value decomposition (SVD), principal component analysis (PCA), and nonnegative matrix factorization (NMF) provide fundamental tools for unsupervised feature extraction tasks [17, 2, 25]. Extensive research has been conducted adapt matrix factorization models to perform classification tasks by supervising the matrix factorization
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

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, nata-reconstructivve, and faas-discriminative features, addressing challenges posed by high-dimensionall data. The general framework of SMF was introduced [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-awe 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, c.g., [33, 30]) as well as the more recent model of convolutional matrix factorization by [22] for 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].
Contributions 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 global minimizer of the objective 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 novel framework that 'lifts' SMF as low-rank matrix estimation problem in 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 thcoretically compare the robustness of filter-based and centure-based SMF, establishing that the former is computationally more robust (see Theorem 3.5) while the latter is statistically more robust (see Theorem 4.
Applying our method to microarray datasets for cancer classification, show that not only it is competitive against benchmark methods, but it able to identify groups of genes including well-known cancer-associated genes (see Fig. 3).
1.1 Notations
Throughout this paper, we denote by R2 the ambient space for data equipped with standard inner project (-, >) that induces the Euclidean norm ||-||- We denote by 1,. A8 the space of class labels with i + classes. For a convex subset in an Euclidean space, we denote IT the projection operator onto O For an integer 2 1, we denote by II, the rank-r projection operator for matrices. xx, For a matrix A = (ds))ss E m we denote its Frobenius, operator (2-), and supremum norm by || = E., dzy. || == suPxER", ||x||= ||A]| || == max;,; la:yl. respectively. For each 1 IsiS S m and S S n, we denote A [ and A [::.j] for the ith row and the jih column of A, identity matrix. For square symmetric matrices A, B E Raxx we denote A < B if v" Av S ''' Bv for all unit vectors v e IR". For two matrices A and B, we denote [A, B] and [A B] the matrices obtained by concatenating (stacking) them by horizontally and vertically, respectively, assuming matching dimensions.
respectively. For each integer n 2 1, I, denotes the
1.2 Model setup
Suppose we are given with n labeled signals (yi, x,, x() for i = 1.... n, where yi E (0, 1,..,,,. x} the label, x; E TRP is a high-dimensional feature of and x E T29 a low-dimensional auxiliary feature of i (p > q). For a vivid context, think of x; as the X-ray image of patient i and x, denoting

---

some biological measurements, such as gender, smoking status, and body mass index. When making predictions of V4, we use a suitable r (< p) dimensional compression of the high-dimensional feature X; well as the low-dimensional feature x as-is. We assume such compression is done by some matrix of (latent) factors W = [wi,..., w,] E TRPXT that is reconstructive in the sense that the observed signals x; can be reconstructed as (or approximated by) the linear transform of the 'atoms' WI....; Wr E _º for some suitable 'code" h, E R". More concisely, Xdata = 1x.....,X. Xn] Rs WH, where H .... hn] |R*x", In practice, we can choose to be the approximate rank of data matrix Xaam (e.g., by finding the elbow of the scree plot).
Model
Interpretable Output
B
Group-Wte association
Logistic Regression
AJ
Sux
01 2.855
0 5,515
Rag.
Rep. Coeff.
860220
0SN151
PIGH1
A1
SLCOSA1
2LO0151
Nark Facor
X
W
H
Training
ATP1A2
14577
Code
TM48-1
KGT014
DPY1
Input
Dictionary
min
W
8/M015
CA17
SLC2542
EP104
Nermal
Canerr
Nomal
Label (e.g. cancer nomai)
Gone
asseciation
X
High-dim feature (e.a, gene, image)
rank(Z) S
Y1 1,653
Coveriste
72 0.164
Augmented lowrrenk
Low-dim feature (e.g.. ago, sex)
Age
Sex
Figure 1: Overall scheme of the proposed method for SMF-H.
Now, we state our probabilistic modeling assumption. Fix parameters W E Rexr, h, E R",
3 TR' and 7 RRa*=_ Let h I [0, oo) be score function (e.g., h(.) exp(.) for multinomial logistic regression). We assume y is a realization a random variable whose conditional distribution is specified
IP (y: x,, xx,,), = x., x)) = g(a:) = C(1,h((,,,) ((a.))........ M(s,..))
where C is the normalization constant and a, = (a.. ac,x) R is the activation for y defined in two ways, depending on whether we use a 'feature-bsedd or 'filter-based' SMF model:
533Th: + y" 7"x a, 1.3 x; xx* for filter-based (SMF-W).
for feature-based (SMF-H),
(2)
One may regard ()) as the 'multinomial regression coefficients' with input feature (h; or (W" x,, x(). In (2), we may regard the code h, (coming from x; & Whe) or the 'filtered signal' WTx, as the r-dimensional compression of x,. Note that these two coincide we have perfect factorization x; Wh; and the factor matrix W orthonormal, i.e., W W I, 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. we are given a test point (Xrest, X the predictive probabilities for its unknown label ytest is given by (1) with activation computed in (2). This only involves straightforward matrix multiplications for SMF-W, which can also be viewed as forward propagation 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(ax)) its correct code representation hes needs to be learned by solving the following 'supervised sparse coding' problem (see [33]):
min min Ely, 3 h + 7 X(x)) + Xiest WhI IE(.. h
(3)
more efficient heuristic testing method for SMF-H is by approximately computing has by only minimizing the second term in (3).
In order to estimate the model parameters (W, H, B, 7) from observed training data (x4, v.) for i L.... n2, we consider the following multi-objective optimization problem:
min W. W.
Ely,, a;) + + Xatta WH 2
(4)
3

---

where Xaaaa Xi...... x, E Rux", a, is as in (2), and E(.) is the classification loss measured by
the negative log-likelihood:
Ely. a) = log 1 Σ h(ac)
(v=c) log h(az).
(5)
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 ç > 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 find factor matrix W that can best reconstruct the feature vectors X via the factorization Xdata 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.
1.3 Related works
The SMF training problem (4) is 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, 47J). However, convergence analysis or statistical estimation bounds of such algorithms are quite limited. Appealing 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 ceature-based SMF model.
The main finding of our work is that the non convexxty of the SMF problem (4) is 'benign', in the sense that there exists an algorithm globally convergent to a global optimun at an exponential rate. We use 'double-lifting' technique that converts the nonconvex SMF problem (4) into a low-rank factored estimation with convex objective. This 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, 51J. Our exponentially convergent SMF algorithms are versions of low-rank projected gradient descent in the algorithm (44) that operate in the double-lifted space.
2 Methods
2.1 Sketch of key idea
Our key idea to solve (4) is to transform it into variant of the low-rank matrix estima- tion problem (6) and then use Low-rank Projected Gradient Descent (LPGD) algorithm (7):
min 2=,0,, (0)
f (Z)
(6)
Z II, IT (Z.- (Z-11-V (7Z-1-)7((4--)))
0 fixed. (7)
In (6), one seeks to minimize objective / w.r.t. paired matrix parameter Z = [0,7] within convex constraint set O and an additional rank constraint rank(0) S r. In (7), 11, denotes applying rank-r projection on the first factor 0 while keeping 7 the same. Problem (6) encompasses variety of significant problems, such as matrix regression [1 [11, 38] and matrix completion [48, 37J. For these problems, algorithms of type (7) have been studied in [41, 56].
To illustrate how the SMF problem (4) transforms into 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 X We seek to solve matrix factorization and linear regression problems simultancously for data matrix X E gxxx and response variable

---

Y Rix,
min |Y W.H.B
33' Eataa WHI
(8)
This three-block optimization problem involving three factors W RRx H E R and B grxi which is nonconvex and computationally challenging to solve exactly. Instead, consider reformulating this nonconvex problem as the following matrix factorization problem:
B (
3 LVEW
2
min f W,
H - VEX Jara
H
(9)
Indeed, we now seek to find two decoupled matrices (instead of three), one for B" and W stacked vertically, and the other for H. A similar idea of matrix stacking was used [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 0 = [ H WH] E (z++))) of rank at most r that minimizes the function f in (9), which is convex (specifically, quadratic) in O: (See Fig. Training).
For SMF-W. W, consider the following analogous linear regression model: min f (W[3, H)) = Y 3² WT X.mmal W.H,
Ell Xamta WH
(10)
where the right-hand side above is obtained by replacing with W in (9). Note that the objective function depends only on the product of the two matrices W and [3. HJ. Then, we may further lift as the low-rank matrix estimation problem by seeking a single matrix [W33, WH] E [Px(1+n) of rank most that solves (6) with being the function in (10).
2.2 Algorithm
Motivated by the observation we made before, we rewrite SMF-H in (4) as
"
min [0,r|E Sr
(0, 7) = Σ Eyn Al; + + Xdata B + A (IAI) + l1>l1))),
(11)
where A = 3 H. B- WH, 0 |A BJ (*++))xn and θ is a convex subset of (R(++)) x R We have added L2 -e-ruullarization term for A and yy with coefficient À 2 0, which will play 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 [0, 71 and then perform rank-r projection of the first factor 0 [A || B) via truncated SVD until convergence. Once we have solution [0*. >Y1 to l), we can use SVD of 0" to obtain a solution to (4). Let 0* = UEVT denote the SVD of 0. Since rank(O") S r. Σ an r X r diagonal matrix of singular values of A. Then U E TR. x)+ and V E R= are oothoorthonormal matrices, that is, UTU VTV I,. Then since 0* [(=))]] W*' we can take H' //²2VT and [(6*)* W 1/
We summarize this approach solving (4) for SMF-H in Algorithm 1. Here. SVD, denotes rank-r truncated SVD and the projection operators II and II, are defined in Subsection 1.1.
As for SMF-W, we can rewrite (4) with additional regulauurrr for A W3 and y as
min (0.y]E
F F0,, 7) = Σ Cly,, ATx x xx) + Xl B
(1A lr>l)), (12)
Sr
where 0 = A, B) = W|B, H] E <RPx(*+*) and θ E (ex((=n) x Rax* is convex set. Algorithm 1 for SMF-W 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 per- iteration complexity is (pn min(n, (n, p)). while that for the nonconvex algorithm is O(ppr + g)n). While the LPGD algorithm is in general more expensive per iteration than the nonconvex method, the iteration complexity only (log e') 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 e-stationary point is at best O(e-1) using standard analysis. Hence for small enough, Algorithm I achieves an e-aceurate global optimum for SMF with a total computational cost comparable to nonconvex SMF algorithm achieve e-stationary point.
5

---

Algorithm 1 Lifted PGD for SMF
Input: Xdata E [Rxx X R (auxiliary features); Y label E {0. kk]
Parameters: > 0 (stepsize); N TEN N (iterations); r 2 (rank); X 2 0 (L2-reg, param.)
Constraints: Convex O C (x+p)xn x RRx for SMF-H. C C [x+ x R for SMF-W W;
Initialize Wo E Rexr, Ho E R Bo R Yo E R
[0.- [ Ho Ho
WoH E ((+x))n
(D for SMF-H)
10 + Wof HoH E RPx(stn)
(o for SMF-W)
for to N do
02 11, (Ile
(0 -V (0R-11
1k- (1-
(i See Appendix A for computation)
Yk f 7k-1 -V 1k-1)
0N N UEVT
(D rank-r SVD)
Wx] +- UE HN < (E)//VVT (o SMF-H)
WN < U, BNv Hx] EVT
( SMF-W)
Output: (WN, HN, BN. YN)
3 Global convergence guarantee
We have discussed that one can cast the SMF problem (4) as the following 'factored estimation problem' minT,,sr7 (TS", 7). Note that such problems generally do not have unique minimizer due to the 'rotation invariance'. Namely, let R be any r x r orthonormal (rotation) matrix (i.e., R'R RR I,). Then (((RR) ((TR)(SR)) 7) S(TRR" -> = S(TST 7)) 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 an exponential rate. First, we introduce the following technical assumptions (3.1-3.3).
Assumption 3.1. (Bounded activation) The activation a E RR* defined in (2) assumes bounded norm, i.e., ||a| < M for some constant M E co).
Assumption 3.2. (Bounded eigenvalues of covariance matrix) Denote & = [ G E R(P+g)x* where o; = [x x
E RP+4 (so &- Xaata Xaxx), where X
Then, there exist constants 8
8 0 such that for all n 1.
6 < Amin Ininnn" "Tp")) S Amax (n 1
(m-1"p"") 5
(13)
Assumption 3.3. (Bounded stiffness and eigenvalues of observed information) The score function
h: IR + 10, oc) is twice continuously differentiable. Further, let observed information H a) := V. V.a. n Ely, a) for y and a. Then, for the constant M > 0 in Assumption 3.1, there are constants Tmax a + s.t. Ymax == | maXice5n
| VaL(vo. a and
a = inf min Amin ()), at := sup max Amax (H.,))) (14) || 15A5
|| ISa5n
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 k x x 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 simple expression for the bounds at, as discussed in the following remark.
Remark 3.4 (Multinomial Logistic Classifier), Let é denote the negative log- likellhoooo function in (5), where we take the multinomial logistic model with the score function h(-) exp(-). Denote (hi,...,hw.) = Val(y, a) and H a) V. V.V.r a)) Then in this special case, we have hy) = (a) 1(y = j) and a), and (30) in Appendix). Under Assumption 3.1, according to Lemma B.1, we can take Tmax s (1+2(xx11)e))
g:(a) (1 =) -)- (a)) (See (28)
--M 1 + Hi--- $ 2, a 14e-" and a classification, at < 1/4.
(1+e (I+e ((p For binary

---

Now define the following quantities:
min(2E, (2E, 2X nō n²)a)) 1 min, 2
max 21 + a+) max(((, 2A + a+)
for SMF-W for SMF-H
(15)
Now, we state special case of our first main result, specifically when the model 'correctly specified", allowing the rank-r SMF model to effectively account for the observed data. This implies the existence of '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 == 10, denote the iterates of Algorithm 1. Assume 3. 3.1.3.3 hold. Let μ and L be as (15), fix (2 27.), and let p == 2(1 TU) E (0, 1). Suppose Lfu < 3 and let Z = 10", ,-²" be any stationary point of F over θ s.l. rank S r. Then Z" is the unique global minimizer of F among all Z 10, 70 with rank r. Moreover, [Zz- 1 S p² [Z ( for 1 1.
101 + lyll Note that we may view the In the statement above, we write 1 10, 1= ratio T./1 that appears in Theorem 3.5 as the condition number of the SMF problem in (4), whereas the ratio I." //1" for 1³ == a and I," == a+a+ 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/μ < by choosing sufficiently large & and so that Theorem 3.5 holds. However, using large L=-regularriztion parameter À may perturb the original objective in (4) too much that the converged solution may not be close to the optimal solution. Hence, may want take À as small as possible. Setting À 0 leads to
L /
0<
3.
L'
<
3p" 2
for SMF- W
3
6
(16)
(2,
< 3
for SMF-H.
min(2ç.0)
0)
For SMF-W W, the multinomial classification problem well-conditioned /u" < 3) and the ratio e/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 Lg-regularization necessary for guarantecing 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.
4 Statistical estimation guarantee
In this section, we formulate a generative model for SMF (4) and state statistical parameter estimation guarantee. Fix dimensions p >> g, and let 2 1 be possibly growing sample size, and fix unknown true parameters B* E Rpxn C* [R ; > E g addition, fix A E gexn for SMF-H and A I for SMF-W. Now suppose that class label, data, and auxiliary features are drawn i.i.d. according to the following joint distribution:
x; = [: + x =C
y4 x ~ Multinomial(L,g((
rrank[[A* || B+] ST for SMF- H. SMF-H. [rank[[A* (Ark[[A',, B*) ST for SMF. W.
(17)
[A 1 + Txx
a, (AA)) x (xx SMF-W,
where each (resp., e) are 1 (resp., y 1) vector of i.i.d. mean zero Gaussian entries with variance 2² (resp., (o We call the above the generative SMF model. In what follows, we will assume that the noise levels and o' are known and focus estimating the four forrmear matrices.
7

---

The (Lz-regularized) negative log-likelihood of observing triples (vr, x,, x,) for i = 1,..., n is given as Cn == F(A, B, + ب XXxx C + c, where is a constant and F is as in (11) or (12) depending on the activation type with tuning parameter { = - - The L2 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 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& Ç = (2o²) This and Theorem 3.5 motivate us to estimate the true parameters B" and Y by the output of Algorithm 1 with & -= (2o²) for (loo n) iterations.
(17), while the problem of estimating
Now we give the second main result. Roughly speaking, it states that the estimated parameter Z, is within the true parameter Z = [A* Y3,""] within Ollog n/ vn) with high probability, provided that the noise variance o² is small enough and the SMF objective (11)-(12) is well-conditioned.
Theorem 4.1. (Statistical estimation for SMF) Assume the model (17) with fixed p. Suppose
Assumptions 3. 1-3.3 hold. Let p., T. be as in (15), p= 2(1 ([]) and O(1) if 7 =ZZF(Z) E
@ and O( O n)) otherwise. Let Z denote the iterates of Algorithm I with the tuning parameter & (2o2)=1 1 -=-regularizatton parameter X > O, and stepsize
3). Then
2p -p' Ioo
sollowing holds with probability at least - :. For all L 2 and n 2 1, 1| - Z
Z F V VzF(L" (
[v/n n+)) provided L/P < Furthermore, logn)gg Ollog n/ vn) if + and o² O(L/n).
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 μ > v log n. Combining with the expression in (15) and the well-conditioning assumption L/p < 3, one needs to require [n)), hence small noise variance o² = O(1/n for SMF-W. However, for SMF-H. this is guaranteed whenever o² o(1/(Vn logn)) and À = J.
5 Simulation and Numerical Validation
5017 (LPG0)
SUIT-W (LPG0I
SMF44 (LUNCD)
SUITW C.700)
[-5.00
10.00
305
10
ISINN
10"1
10
Training 3od0
10"
0.13
10
20
11
Lo-3
10-4
10-3
10-2
200
1.0
10
suupodd goT
10-4
10-1
10-3
10:1
0.5 1.5 2.0
20
10
20
Elapoet timan
Elapsed
Elapoed time
Elapsed time (5)
Elapsed time (s)
SNEWWORSE
a Semi-synthetic MNIST data
b Jobo postings data
cLPGD BCB
Figure 2: (a-b) Training loss vs. clapsed CPU time for Algorithm 1 (with binary logistic classifier) on the semi-synthetic MNIST and Job postings datasets for several values of E log scale. Average training loss over ten runs and the shades representing the standard deviation shown. (c) Comparison between LPGD (Algorithm 1) and BCD algorithms for SMF.
We numerically verify Theorem 3.5 on a semi-synthetic dataset generated by using MNIST image dataset [24] (p 28 784, q 0, 500, K 1) and text dataset named 'Real / Fake Job Posting Prediction [1] (p 2840, q = 72, n = 17880, x = 1). Details about these datasets are in for MNIST and 20 for job postings
Sec. G in Appendix. We used Algorithm with rank
datasets. For all experiments, >= = 2 and stepsize 7 = 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
We provide our implementation of Algorithm I in our code repository SMF/ htt//tre//ma/ /main.
8

---

Ć E 1,, 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 p = (1 - TH²,, which decreases ( since y increases in & (see (15)). The decay for large & E {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 W models at & E (5, 103 (other values of & omitted).
6 Application: Microarray Analysis for Cancer Classification
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-leaming approaches. One consists of 54,676 gene expressions from 51 subjects with binary labels indicating pancreatic cancer; Another we use has 35,987 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 binary logistic classifier trained using Algorithm 1 - against the following benchmark algorithms: SMF- and SMF-H trained using BCD; 1-dimensional seven- layer 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 Xarain 23 UVV2 and take W = UΣ and H = vT. For testing, we use WT X input to logistic regression for both filter and feature methods since |[ WHesal F is minimized when Hest = (WTW)- W T = W Xest with orthogonal W.
Known cancer-asodiated
Breast Cancer Datection
Pancreatic Cenoer Diteccion
genes detected by our method
Supervised Princ pal Gene Groups
(Accuracy=0.85, 0.85, [-score=00883 85)
C
Supervised Principel Gene Groups (Accuracy=0.9e 96, F-score=0.97)
Breast
Cancer
Parcrrattic
Cancar
Group-wise associafion Group-vviee asscciation
Group association Group-wise assuciation
Favorable Prognostic Marker
10.480
0.525
2.855
5.516
(PTERTAS)
GRK
PALD1
STXXPSL
(D82))
CAL2
G032
C2C62L
SEC120
CSN191
(RAGGEFLD
Cowe
GAPOEPLT
RMRM1
SANGS
ARTAP113
(SHBOL)
NF2
SLOCBA1
ATP112
KI446403
KD45
(HHDPL)
(HHITZ
OPY11LL
(KTT044)
SPR45S
UN0002398
APTR
BRCAT
EPHK4
GOLGANGP
APA
BRCAD
117AA
OPY18L
1M45F
(Oiscogene)
LIN000663
RN
ATR201
USP10
9AMD15
(CA1 2
1144
Unfevorable Prognostic Marker
Carsrrdte aaacedel oa
Narmal Cencer
Normal t-s Cancer
Metvcdi
Paaowado
Eivan
Unsupervised Principal Gene Groups
d
Unsugervised Principal Gere Groups
SMF- (wes))
0.886 (8.5)
0,915 (001)
(Acouracy=0.72 F-ecore=0.71
(Accuracy=0.773 73, F-ecore-0,
SMF
(SC3))
0.88 (310)
Grouprstee resociaaion
Groupwisee mosociatioos
Group wise association
Group Sc association
CNN
0.090 (237)
UL53 (U0h)
MB--31
I1454 (1.27)
CoCL5
RASALL
120
BUIIC
0815-1300)
(L810 (802)
HF 0.81
0.544 (0/2)
(a-b) Two selected superviseddunsuperviin principal gene groups (low-dimensional compression of genes) learned by rank- 16 SMF- W/SVD and their associated logistic regression coefficients for breast cancer detection. (c-d) Similar to a-b learned by rank-2 SMF-W/SVD for pancreatic cancer detection. (e) Blue-circled genes within each gene group of extreme coefficients coincide with known prognostic markers (for pancreatic cancer) and oncogene (for breast cancer). (f) Average classification accuracies and their standard deviations (in parenthesis) for various methods on two cancer microarray datasets over five-fold cross-validation. The fighest-performing instances are marked in bold.
Figure 3:
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 times. A scree plot is used to determine the rank r. Other parameters are chosen through 5-fold cross-validation (g E {0.1, 1, 10) and X E (0.1, 1, 10)), and the algorithms are repeated 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 i. 1) and learn factor matrix W [W1,.... wr] E TRxx and vector of regression coefficients 3 = [ 8,] E RRIx. each column w, of W describes a latent factor and the corresponding regression coefficient, 3 tells us how w, associated with class labels. The pairs By), 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 corresponds to a weighted group of genes (which we call 'principal gene group') and B, 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.
7 Conclusion and Limitations
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-basee SMF, finding that the former is computationally more robust while the latter is statistically more robust. The algorithm's 's exponential convergence is numer- ically 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 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 wider range of datasets.
References
[1] Real / fake job posting prediction. htttps://ww..in http//wwww.gglee..e//// wwt com/datasets//hiit wwtiaasstssssssiiiiieeeeee saabb real-or- f-lf-aae--arae fake-aaeeejor jaw----bboooopreddiiee prediction. Accessed: 2023-05-16.
[2] Hervé Abdi and Lynne J Williams. Principal component analysis. Wiley interdisciplinary reviews: computational statistics, 2(4):433-459, 2010.
[3] Alekh Agarwal, Sahand Negahban, and Martin J Wainwright. Fast global convergence rates of gradient methods for high-dimensional statistical recovery. Advances Neural Information Processing Systems, 23, 2010.
[4] Sanae Amani and Christos Thrampoulidis. Ucb-based algorithms for multinomial logistic regression bandits. Advances in Neural Information Processing Systems, 34, 2021.
[5] Woody Austin, Dylan Anderson, and Joydeep Ghosh. Fully supervised non-negative matrix factorization for feature extraction. In IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium, pages 5772-5775. IEEE, 2018
[6] Shivam Bansa. fake job postings dataset. https: // www. kaggle, com/ datasets/ shivamb/ (eaamb/rell-a fake- fake- patee-joooosttnppeeeee prediction, 2017.
[7] Amir Beck. First-order methods in optimization. SIAM, 2017.
[8] Christopher M Bishop et al. Neural networks for pattern recognition. Oxford university press, 1995.
10

---

[9] Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning. volume 4. Springer, 2006.
[10] Dankmar Böhning. Multinomial logistic regression algorithm. Annals of the institute of Statistical Mathematics, 44(1):197-200) 1992.
[11] Emmanuel Candes and Benjamin Recht. Exact matrix completion via convex optimization. Communications of the ACM, 55(6):11-111,, 2012.
[12] Moody T Chu, Robert E Funderlic, and Robert Plemmons. Structured low rank approximation. Linear algebra and its applications, 36::156-6110 157-172, 2003.
[13] Michael Elad and Michal Aharon. Image denoising via sparse and redundant representations over learned dictionaries. IEEE Transactions on Image processing, 15(12):373-33755 2006.
[14] Bruno César Feltes, Eduardo Bassani Chandelier, Bruno lochins Grisci, and Márcio Dom. Cumida: An extensively curated microarray database for benchmarking and testing of machine learning approaches in cancer research. Journal of Computational Biology, 26(4):376-386 2019. PMID: 30789283.
[15] Mchrdad J Gangeh, Ahmed K Farahat, Ali Ghodsi, and Mohamed S Kamel. Supervised dictionary learning and sparse representation-a review. arXiv preprint arXiv: 2015,
[16] Mehrdad J Gangeh, Pouria Fewzee, Ali Ghodsi, Mohamed S Kamel, and Fakhri Karray. Multi- view supervised dictionary learning in speech emotion recogmition. JEEE/ACM Transactions on Audio, Speech, and Language Processing, 20666 2014.
[17] Gene H Golub and Christian Reinsch. Singular value decomposition and least squares solutions. In Linear algebra, pages 134-151. Springer, 1971.
[18] Luigi Grippo and Marco Sciandrone. On the convergence of the block nonlinear gauss-seidel method under convex constraints. Operations research letters, 2673):127-136, 2000.
[19] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012.
[20] Prateck Jain, Raghu Meka, and Inderjit Dhillon. Guaranteed rank minimization via singular value projection. Advances in Neural Information Processing Systems, 23, 2010.
[21] Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi. Low-rank matrix completion using alternating minimization. Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 665-674, 2013.
[22] Donghyun Kim, Chanyoung Park, Jinoh Oh, Sungyoung Lee, and Hwanjo Yu. Convolutional matrix factorization for document context-aware recommendation. In Proceedings of the 10th ACM conference on recommender systems, pages 233-240, 2016.
[23] Guillaume Lecué and Shahar Mendelson. Sparse recovery under weak moment assumptions. Journal of the European Mathematical Society, 19(3):881-944, 2017.
[24] Yann LeCun and Corinna Cortes. MNIST handwritten digit database. 2010.
[25] Daniel D Lee and H Sebastian Seung. Learning the parts of objects by non-negative matrix factorization. Nature, 401(6755);788, 1999.
[26] Johannes Leuschner, Maximilian Schmidt, Pascal Fernsel, Delf Lachmund, Tobias Boskamp. and Peter Maass. Supervised non-negative matrix factorization methods for maldi imaging applications. Bioinformatics,
2019.
[27] Zhengming Li, Zheng Zhang, Jie Qin, Zhao Zhang, and Ling Shao. Discriminative fisher embedding dictionary learning algorithm for object recognition. IEEE transactions on neural networks and learning systems, 31(3):786-800, 2019.
[28] Julien Mairal. Optimization with first-order surrogate functions. International Conference on Machine Learning, pages 783-791, 2013.

---

[29] Julien Mairal. Stochastic majorization inpiorizatoo-mmmimizattooo algorithms for large-scale optimization. In Advances in Neural Information Processing Systems, pages
2013.
[30] Julien Mairal, Francis Bach, and Jean Ponce. Task-driven dictionary learning. IEEE transactions on pattern analysis and machine intelligence, 34(4):791-804, 2011.
[31] Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro. Online learning for matrix factorization and sparse coding. Journal of Machine Learming Research, 11(Jan):1-60,, 2010.
[32] Julien Mairal, Michael Elad, and Guillermo Sapiro. Sparse representation for color image restoration. IEEE Transactions on Image Processing, 17(1):53-66 2007.
[33] Julien Mairal, Jean Ponce, Guillermo Sapiro, Andrew Zisserman, and Francis Bach. Supervised dictionary learning. Advances in Neural Information Processing Systems, 21:1033- 2008.
[34] Raghu Meka, Prateek Jain, and Inderjit Dhillon. Guaranteed rank minimization via singular value projection. arXiv preprint 4XXv 2009.
[35] Fionn Murtagh. Multilayer perceptrons for classification and regression. Neurocomputing, 115-6) 1991.
[36] Sahand Negahban and Martin Wainwright. Estimation of (near) low-rank matrices with noise and high-dimensional scaling. The Annals of Statistics, 39(2): 2011.
[37] Sahand Negahban and Martin Wainwright. Restricted strong convexity and weighted ma- trix completion: Optimal bounds with noise. The Journal of Machine Learning Research, 13(1):1665-116) 2012.
[38] Sahand Negahban, Bin Yu, Martin J Wainwright, and Pradeep Ravikumar. A unified framework for high-dimensional analysis of m-estimators with decomposable regularizers. Advances in neural information processing systems, 22, 2009.
[39] Yu Nesterov. Gradient methods for minimizing composite functions. Mathematical program- ming, 140(1):125-161, 2013.
[40] Dohyung Park, Anastasios Kyrillidis, Srinadh Bhojanapalli, Constantine Caramanis, and Sujay Sanghavi. Provable non-convex projected gradient descent for class of constrained matrix optimization problems. stat, 1050:4, 2016.
[41] Dohyung Park, Anastasios Kyrillidis, Constantine Caramanis, and Sujay Sanghavi. Finding low-rank solutions via nonconvex matrix factorization, efficiently and provably. SIAM Journal on Imaging Sciences, 11 11(4):1165-2244 2018
[42] Dohyung Park, Anastasios Kyrillidis, Constantine Carmanis, and Sujay Sanghavi. Non-square matrix sensing without spurious local minima via the burer-monteiro approach. Artificial Intelligence and Statistics, pages 65-74. PMLR, 2017.
[43] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit- learn: Machine learning in python. the Journal of machine Learning research, 12:2825-2830, 2011.
[44] Gabriel Peyré. Sparse modeling of textures. Journal of Mathematical Imaging and Vision, 34(1):17-31, 2009.
[45] Pradeep Ravikumar, Martin J Wainwright, Garvesh Raskutti, and Bin Yu. covariance estimation by minimizing ll ll-penalized deg-retermnnan divergence. Electronic Journal of Statistics, 5:935 980, 2011.
[46] Benjamin Recht, Maryam Fazel, and Pablo A Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM review, 5203),77-501, 2010.
[47] Alexander Ritchie, Laura Balzano, Daniel Kessler, Chandra S Sripada, and Clayton Scott. Supervised pea: A multiobjective approach. arXiv preprint arXiv:2011.05309, 2020.
12

---

[48] Angelika Rohde and Alexandre B Tsybakov. Estimation of high- low-rank matrices. 2011.
[49] Evelina Sjöstedt, Wen Zhong, Linn Fagerberg, Max Karlsson, Nicholas Mitsios, Csaba Adori, Per Oksvold, Fredrik Edfors, Agnieszka Limiszewska, Feria Hikmet, al. An atlas of the protein-coding genes in the human, pig, and mouse brain. Science, 36766482) 2020.
[50] Guan Tian and Ting-Zhu Huang. Inequalities for the minimum eigenvalue of m- m-matrices. The Electronic Journal of Linear Algebra, 20:291-302, 2010.
[51] Tian Tong, Cong Ma, and Yuejie Chi. Accelerating ill-conditioned low-rank matrix estimation via scaled gradient descent. Journal of Machine Learning Research, 22(150))1-63, 2021.
[52] Joel A Tropp et al. An introduction to matrix concentration inequalities. Foundations and Trendso Machine Learning, 8(1-2):1-230 2015.
[53] Stephen Tu, Ross Boczar, Max Simchowitz, Mahdi Soltanolkotabi, and Ben Recht. Low-rank solutions of linear matrix equations via procrustes flow. In International Conference on Machine Learning, pages 964-973. PMLR, 2016.
[54] Roman Vershynin.
probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018
[55] Michael E Wall, Andreas Rechtsteiner, and Luis M Rocha. Singular value decomposition and principal component analysis. In A practical approach to microarray data analysis, pages 91-109. Springer, 2003.
[56] Lingxiao Wang, Xiao Zhang, and Quanquan Gu. unified computational and statistical framework for nonconvex low-rank matrix estimation. In Artificial Intelligence and Statistics, pages 981-990. PMLR, 2017.
[57] Stephen J Wright. Coordinate descent algorithms. Mathematical Programming, 151(1):3-34, 2015.
[58) Yangyang Xu and Wotao Yin. A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion. SIAM Journal imaging sciences, (33):1758-178) 2013.
[59] Yael Yankelevsky and Michael Elad. Structure- aware classification using supervised dictionary learning. In 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4421-4425 IEEE, 2017.
[60] Pavel Yaskov. Controlling the least eigenvalue of a random gram matrix. Linear Algebra and its Applications, 304:108-122 123, 2016.
[61] Qiang Zhang and Baoxin Li. Discriminative k-svd for dictionary learning in face recognition. In 2010 IEEE computer society conference on computer vision and pattern recognition, pages IEEE, 2010.
[62] Shijie Zhao, Junwei Han, Jinglei Lv, Xi Jiang, Xintao Hu, Yu Zhao, Bao Ge, Lei Guo, and Tianming Liu. Supervised dictionary learning for inferring concurrent brain networks. IEEE transactions on medical imaging, 34(10):0036-2045, 2015.
[63] Tuo Zhao, Zhaoran Wang, and Han Liu. nonconvex optimization framework for low rank matrix estimation. Advances in Neural Information Processing Systems, 28:559, 2015.
[64] Qinqing Zheng and John Lafferty. A convergent gradient descent algorithm for rank mini- mization and semidefinite programming from random linear measurements. arXiv preprint arXiv:150..00881 2015.
13


Reviews:
[{"review_id": "xTOr8j1Hzr", "reviewer": "Reviewer_pg2E", "summary": "This paper proposes a variant of supervised dictionary learning (SDL), provides some theoretical guarantee on finding the global minimizer of the problem with arbitrary initialization, and showcase its application in pancreatic cancer.", "questions": "See above.", "limitations": "See above.", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "Theoretical analysis on dictionary learning (DL) has a rich literature. However, theoretical analysis on convergence of supervised DL is scarce. So, this is of interest. The paper is written clear.", "weaknesses": "The paper lacks a proper literature, introduction, and formulation of the dictionary learning problem. Several statements in this paper are not properly claimed. The paper in its formulation is not rigorous. Here are some indication of this.\n\n1. Abstract: \n\t- SDL is not properly defined. There is no info on its key property which is trying to represent data with a sparse combination of columns of a dictionary, which itself is learned.\n\n\t- The global convergence with random initialization is a bold statement. Indeed, this is misleading, as this paper (I explain below) is solving a different problem than DL or SDL. I note that all theoretical analysis on DL work with the assumption that the initial estimation of the dictionary must be close to the true one [1,2,3,4] (to name a few)\n\n\t- There is no reference to the sparsity characteristic of sparse coding/dictionary learning. So, then what is the motivation behind using SDL and not another method?\n\n2. Line 24, DL is not introduced properly. DL is a method trying to learn a sparse representation that can explain the data through a dictionary in a generating fashion.\n\n3. Line 26, This is confusing. DL is a specific model with sparsity. NMF, PCA, ... are different models. However, the statement is mixing all up.\n\n4. Line 26, These citations are proper but old (there are many other recent works on DL).\n\n5. Line 41, DL extract a high-dimension sparse features (D is overcomplete) not a low-dim feature.\n\n6. Line 48, Literature on theoretical analysis on convergence and recovery for DL is missing.\n\n7. Line 50, Theoretical analysis on DL have generative perspective such that there is a sparse representation that has generated the data (or in the SDL case, is also mapped to a class label), and the goal is given data, learn the dictionary and recover that sparse representation. DL is not about global minimizer of an objective (as the solution to overcomplete dictionary may not be unique without sparsity).\n\n8. Why there is no notion of sparsity on H in (3)? Same on line 117.\n\n9. Related work is missing DL literature. Moreover, DL is a bi-convex problem.\n\n10. Line 154 is bi-convex in (beta, W) and (H). (6) objective is exactly similar to the one above it. Moreover, the objective is bi-convex not convex.\n\n11. For section 3, theoretical analysis on SDL may aim to recover the sparse representation that has generated the data, and has produced the label through a probabilistic generative model. This paper is solving another problem, hence their convergence analysis does not address the SDL problem.\n\nOverall, this paper's analysis is not on recovery of the sparse representation that has generated (explained) the data. Hence, statements are confusing in this regard. The paper provides convergence on another parameter theta.\n\n\n[1] Agarwal, A., Anandkumar, A., Jain, P., & Netrapalli, P. (2016). Learning sparsely used overcomplete dictionaries via alternating minimization. SIAM Journal on Optimization, 26(4), 2775-2799.\n\n[2] Chatterji, N. S., & Bartlett, P. L. (2017). Alternating minimization for dictionary learning: Local convergence guarantees. arXiv preprint arXiv:1711.03634.\n\n[3] Rambhatla, S., Li, X., & Haupt, J. (2018, September). NOODL: Provable Online Dictionary Learning and Sparse Coding. In International Conference on Learning Representations.\n\n[4] Arora, S., Ge, R., Ma, T., & Moitra, A. (2015, June). Simple, efficient, and neural algorithms for sparse coding. In Conference on learning theory (pp. 113-149). PMLR.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your evaluation of our work, which greatly contributes to its enhancement. We value your expertise in the field and understand your concern about the term DL being closely linked to sparse representation. While we acknowledge the prevalent association of DL with recovering overcomplete dictionaries and sparse representations, our primary focus is on *undercomplete/low rank* SDL in *high-dimensional* setting.\n\nWe would like to request the opportunity to delve deeper into our approach's motivations and implications either during the rebuttal process or through subsequent discussions. At a minimum, we are committed to clearly distinguishing our high-dimensional SDL via undercomplete DL from the conventional overcomplete DL. To facilitate reader comprehension, we are open to altering our terminology. Terms like \"supervised undercomplete dictionary learning\" or \"supervised matrix factorization\" could aptly encapsulate our approach's essence and better distinguish it from the supervised version of the standard overcomplete DL with sparse representation.\n\nIn response to your feedback:\n\n**(1)-(2)**\n\n**Motivation for Undercomplete SDL**: Our motivation for undercomplete dictionaries in SDL stems from the classification of high-dimensional data (e.g., genomic data analysis for cancer classification). We aim to learn a low-rank basis that offers interpretable, data-reconstructive, and class-discriminative features, addressing challenges posed by high-dimensional data. The overcomplete dictionary approach, with its associated sparse representation, proves computationally infeasible and complex to interpret. Our preference for undercomplete dictionaries aligns with our goals. Additionally, we find that assuming sparse representation over an undercomplete dictionary for high-dimensional data isn't always reasonable. As a result, we've excluded the sparsity regularizer for the code matrix $H$.\n\n**Related Literature on Undercomplete/Low-rank DL**: Although most DL literature emphasizes overcomplete dictionaries and sparse representation, notable works explore undercomplete or low-rank dictionaries for various purposes (examples refs below). Notably, [1] introduces undercomplete dictionaries for efficient SDL tasks on high-dimensional data, resembling our approach (without theoretical guarantees). In [2], the identifiability of sparse component analysis is examined under an undercomplete dictionary. [3] employs undercomplete dictionary learning for low-rank dictionary learning, leading to unsupervised low-rank feature extraction.\n\n[1] Mohseni-Sehdeh et al. \"A Fast Dictionary-Learning-based Classification Scheme Using Undercomplete Dictionaries.\" Signal Processing (2023): 109124.\n\n[2] Cohen et al. \"Identifiability of complete dictionary learning.\" SIAM Journal on Mathematics of Data Science 1.3 (2019): 518-536.\n\n[3] Parsa et al. \"Low-rank dictionary learning for unsupervised feature selection.\" Expert Systems with Applications 202 (2022): 117149.\n\n**(3)**\nWe'll revise this to \"matrix factorization\" to better encompass a broader spectrum of basis-learning problems.\n\n**(4)**\nWe'll enhance our reference section with more recent DL sources.\n\n**(5)**\nDL with overcomplete dictionary does extract high-dimensional features, but DL with either low-rank or undercomplete dictionary extracts low-dimensional features. We will clarify this point. \n\n**(6)**\nAcknowledging the importance of theoretical analysis on DL, we'll incorporate the references you suggested, which delve into aspects like local convergence and recovery guarantees.\n\n**(7)**\nThank you for your insightful comment. The classical references [30, 33] define the SDL problem as an optimization problem, with the possibility of undercomplete $r<p$ dictionary not excluded. Our formulation of SDL (eq. (3)) mirrors [30] (eq. (4)) without sparseness regularizer under $r<p$. Our work guarantees that one can find *some* global minimizer of the non-convex SDL objective exponentially fast from any initialization (under some conditions). This is by no means about recovering the ground-truth dictionary and representation. In fact, as we noted in the first paragraph of Sec. 3 (and also as the reviewer correctly pointed out), our SDL optimization problem does not have a unique global minimizer. However, if we transform separate matrix factors into a combined low-rank matrix (denoted $\\theta$), then uniqueness is guaranteed. Under this setting, we also obtain statistical estimation guarantee under generative SDL models, see Sec. 4. \n\n**(8)**\nSee our response for 1 and 2.\n\n**(9)**\nWhile we're constrained by space, we concur on the importance of showcasing recent DL advances. These references would fortify our work's positioning.\n\n**(10)**\nAcknowledging the non-convexity of (6), if we introduce the combined low-rank matrix factor $\\theta$ (in line 162), then the objective function in (6) is quadratic $\\theta\\mapsto \\lVert \\theta_{0} - \\theta \\rVert_{F}^{2}$ for a fixed $\\theta_{0}$. One can then minimize this convex function under a low-rank constraint on $\\theta$,  yielding solutions for $\\beta$, $W$, and $H$. These points are stated in lines 160-163. We will further clarify our discussion in lines 156-163. \n\n**(11)**\nWhile the recovery of the ground truths is a natural goal of an SDL that draws a direct analogy from the literature of overcomplete DL, it is not what we aim for in this work. Our underlying assumption is that the observed high-dimensional labeled data is generated by a linear combination of unknown low-rank features. In this setting, there are infinitely many, equally effective undercomplete dictionaries and representations (not necessarily sparse) that could have generated the observed labeled data. Hence our goal is to find an optimally effective pair of class-discriminative dictionaries and data representation by solving an optimization problem. \n\nWe sincerely hope our clarifications assuage your concerns, encouraging your reconsideration. "}, {"review_id": "Hcup6kQLLc", "reviewer": "Reviewer_zTTY", "summary": "This paper explores the optimization of supervised dictionary learning, focusing on a non-convex objective function with a matrix factorization structure. The authors propose a reformulation of the problem as a minimization task with a low-rank constraint. To solve this reformulated problem, they employ a projection gradient descent-type method. The paper provides convergence results and statistical analysis, along with practical applications to showcase the efficacy of the proposed algorithm.", "questions": "See above.", "limitations": "Yes.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Optimizing non-convex structured problems poses intriguing and challenging tasks. The problem addressed in this paper represents a classic learning paradigm in the machine learning field, making the authors' contribution well-motivated. Additionally, they demonstrate the effectiveness of their new PGD-variant through classification tasks using medical data. The paper is well-written, presenting its ideas coherently.", "weaknesses": "Overall, the paper is of a standard quality. Considering the extensive research on optimization problems with low-rank matrix factorization structure, the newly proposed techniques appear nontrivial yet not entirely surprising. My primary concern lies with the claim in the abstract (and other sections of the paper, e.g., L55) that the proposed method \"provably converges exponentially fast to a global minimizer of the objective.\" Here are my questions:\n\n* [Clarification question A] It appears that the global convergence result (Theorem 3.5) holds only when a low-rank stationary point Z^* exists, where \"stationary\" is defined as first-order optimality with respect to F under the convex constraint \\Theta. In such a case, the objective function seems to be already strongly convex, and the mentioned existence might imply that the low-rank constraint could be eliminated (please correct me if I am mistaken). This assumption is quite strong and essentially transforms a non-convex problem into a convex one. I believe this could be a significant limitation of the analysis.\n\n* [Clarification question B] I acknowledge that the authors present a general convergence result in Theorem D.1. My second question pertains to Theorem D.1(ii). It seems that, under the \"possibly misspecified case,\" Theorem D.1(ii) cannot even guarantee the sequential convergence of {Z_t} since the residual term (second term on the right-hand side of Equation (91)) might not be zero even at the global optimal Z^*. I fail to see (please correct me if I am mistaken) why the gradient mapping would be zero at optimal Z^*, as this mapping does not consider the normal cone (defined in a certain generalized sense) of the rank constraint.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**`Clarification question A`**\n\n**Response:** We would like to express our gratitude for the thoughtful comments provided by the reviewer. To begin, we wish to emphasize that our theoretical analysis comprises three parts:\n> 1. (Thm. C.2) Establishing exponential convergence results for the general low-rank projected gradient descent (LPGD) algorithm for $r$-restricted strongly convex (RSC) and $r$-restricted smooth (RSM) objectives.\n\n> 2. (Thm. 3.5, Thm. D.1) Conducting an extensive second-order analysis of the SDL objective to verify that the reformulated SDL problems satisfy the hypotheses of Thm. C.2.\n\n> 3. (Thm 4.1) Establishing statistical estimation guarantee for SDL under generative models by using previous computational guarantees (Thm D.1).\n       \nThe reviewer's concern appears to stem from mixing the two statements in Thm 3.5 (on SDL) and Thm C.2 (on LPGD). Indeed, the objective functions in Thm 3.5 correspond to the lifted SDL objectives in eq. (8) and (9). Notably, these objectives are not inherently ($r$-restricted) strongly convex. In fact, verifying that the lifted SDL objectives do satisfy the RSC/RSM properties is the content of the contribution 2 mentioned above. \n\nFurthermore, we emphasize that the proof of Thm. C.2.(i) for the correctly specified case, where the assumed RSC/RSM objective has a rank-$r$ stationary point, is nontrivial, since the strong convexity holds only when we restrict the objective on matrices with rank at most $r$. Thus one cannot eliminate the low-rank constraint, as opposed to the reviewer's comment. \n\nMore specifically, even if the objective function admits a low-rank stationary point $\\mathbf{Z}^{*}$, there could be many other stationary points among matrices with rank $>r$. Moreover, since our LPGD algorithm in general goes in (after low-rank projection) and out (after gradient descent in the ambient space) of the low-rank space, one needs to carefully address that the possibly wild landscape outside of the low-rank space does not significantly impact the convergence. \n\nIn addition, Theorem 3 of Zhu et al. '18 (see below) shows that if an objective function with a matrix input is both RSC and RSM, with a condition number $L/\\mu<3/2$, and if it admits a low-rank critical point (without additional convex constraint), then there are no spurious local minima in the factored parameter space: \n\n> [1] Zhihui Zhu, Qiuwei Li, Gongguo Tang, and Michael B Wakin. *Global optimality in low-rank matrix optimization.* IEEE Transactions on Signal Processing, 66(13):3614–3628, 2018.\n\nThis result also indicates that minimizing an RSC function with a low-rank critical point is a nontrivial problem. To put in context, our Thm C.2(i) shows that under the hypothesis of a condition number $L/\\mu<3$ (weaker than in Zhu et al.) and with a convex constraint, the assumed low-rank stationary point $\\mathbf{Z}^{\\star}$ is the global minimizer of the objective among the low-rank matrices and that the LPGD algorithm converges to $\\mathbf{Z}^{\\star}$. \n\n\n**`Clarification question B`**\n\n**Response:** Thank you for your thoughtful question. \n\n We do see the reviewer's concern regarding the claimed global exponential convergence to the global minimizer of the SDL problem. The exponential convergence technically holds for the correctly specified case, and for the misspecified case, the exponential convergence is up to an additive error that depends on the extent of misspecification. We will revise the statement in the abstract and throughout the paper to clarify this point. We would also like to bring the Reviewer’s attention to two relevant points below.\n\nFirst, even in the presence of a nonzero misspecification error (which bounds the unnormalized estimation error $\\lVert \\mathbf{Z}^{\\star}-\\mathbf{Z}_\\text{$\\infty$} \\rVert_F$, $\\mathbf{Z}^{\\star}\\in \\mathbb{R}^{p\\times n}\\times \\mathbb{R}^{q\\times n}$), our Thm. 4.1 demonstrates that, under natural generative models for SDL, this error becomes vanishingly small with high probability with noise variance $\\sigma^{2}=O(1/n)$ for SDL-$\\mathbf{W}$ and $\\sigma^{2}=o(1/\\sqrt{n})$ for SDL-$\\mathbf{H}$. Roughly speaking, these results indicate that the generative SDL models are nearly correctly specified with high probability. As a result, our algorithm achieves exponential convergence to the correct parameters for the generative SDL model up to a statistical error that vanishes as the sample size $n$ tends to infinity.\n\nSecondly, it seems to be challenging to exactly recover the global minimizer of an RSC function under a low-rank constraint even when there is no additional convex constraint. To the best of our knowledge, all existing works on similar low-rank matrix estimation problems without additional assumptions (e.g., incoherence assumption for matrix sensing problems) recover global optimum up to a misspecification error, which is zero when there is a low-rank stationary point (correctly specified) but is nonzero otherwise. This is the case, for example, in the following references:\n   \n>[2] Lingxiao Wang, Xiao Zhang, and Quanquan Gu. *A unified computational and statistical framework for nonconvex low-rank matrix estimation.* In Artificial Intelligence and Statistics, pages 981–990. PMLR, 2017.\n\n>[3] Dohyung Park, Anastasios Kyrillidis, Constantine Caramanis, and Sujay Sanghavi. *Finding low-rank solutions via nonconvex matrix factorization, efficiently and provably.* SIAM Journal on Imaging Sciences, 11(4):2165–2204, 2018.\n\n>[4] Sahand Negahban and Martin J Wainwright. *Estimation of (near) low-rank matrices with noise and high-dimensional scaling.* The Annals of Statistics, 39(2):1069–1097, 2011."}, {"review_id": "RrQiUomA9c", "reviewer": "Reviewer_bLMw", "summary": "The work focuses on a generic supervised dictionary learning formulation, which is convex in each of the variable-blocks but not overall. The idea is to stack up the matrix-variables  and obtain an equivalent low-rank optimisation problem with a convex objective. And, then this is solved using a projected gradient descent algorithm (algo1).\n\nInterestingly, under mild conditions, exponentially fast convergence to global minimiser for algo1 is proven. This is a far stronger guarantee than the popular plain block-coordinate descent (BCD), which (in general case) guarantees only stationary-point convergence.\n\n As a side-result, assuming a generative model for the data and under low noise conditions (<O(1/n)), statistical consistency of the sample based-optimal solution is also presented. \n\nThe convergences rates are empirically verified on synthetic data. Interestingly in fig3a, it is shown that just by using the proposed pgd solver instead of the popular BCD, one may achieve considerable improvement in application performance.", "questions": "Following points may improve the presentation of this already good submission\n\n1. I guess the order in which we perform the projection \\pi_\\Theta and \\pi_r matters for the rate of convergence. Also, as noted in the paper, works like [12] alternate between projections. Any discussion on this may help the reader appreciate the algorithm better. For example if one does \\pi_r followed by \\pi_\\Theta, then what may happen? Why not alternate between projections? do these make the convergence analysis difficult or is it that the convergence actually slows down if the order of projection is tampered with ?\n\n2. Will a distribution free analysis in section 4 be more interesting ? With feasibility set being bounded and the objective in terms of the data is lips.conts., this should be possible.\n", "limitations": "Limitations are clearly mentioned in the paper.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. While BCD style of algorithms are popularly studied in this context, I think studying the lifted-pgd is interesting.\n2. Theorem 3.5 result is interesting and potentially can be applied elsewhere.\n3. Improvement over BCD in fig3a is considerable, highlight the impact of the study.", "weaknesses": "1. Though theoretical guarantees are proven, from simulations it is not clear how much is the improvement of pgd vs baseline (BCD) in terms of time. For e.g. in figure 2, it would have been very helpful if BCD is also included.\n\n\n\nVery minor comment:\n1. The presentation of Fig 2 & Eqn 9 could be adjusted to not interfere with the text.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate the reviewer’s overall positive comments. \n\n> Though theoretical guarantees are proven, from simulations it is not clear how much is the improvement of pgd vs baseline (BCD) in terms of time. For e.g. in figure 2, it would have been very helpful if BCD is also included.\n\n**Response:** We sincerely thank your suggestion here! We added comparison with BCD in Figure 2 in the revision. We also included the revised Figure 2 in the one-page supplementary PDF for rebuttal so that the reviewer can check. \n\n\n> Very minor comment: The presentation of Fig 2 & Eqn 9 could be adjusted to not interfere with the text.\n\n**Response:** Fixed.\n\n\n>Questions: I guess the order in which we perform the projection $\\Pi_{\\Theta}$ and $\\Pi_r$ matters for the rate of convergence. Also, as noted in the paper, works like [12] alternate between projections. Any discussion on this may help the reader appreciate the algorithm better. For example if one does $\\Pi_r$ followed by $\\Pi_\\Theta$, then what may happen? Why not alternate between projections? do these make the convergence analysis difficult or is it that the convergence actually slows down if the order of projection is tampered with ?\n\n\n**Response:** Thanks for these great questions. First, in principle, one could alternate between the two projections $\\Pi_{r}$ and $\\Pi_{\\Theta}$ at every iteration after a gradient descent step until convergence, similarly to the alternating projection in [12]. However, this would make each iteration of the algorithm prohibitively expensive as this requries to perform rank-$r$ SVD until convergence at *every iteration*. The problem in [12] is much simpler than ours as the objective function is simply the Frobenius norm between the target and the estimated low-rank constrained matrix. \n\nSecond, it is a great insight to switch the order of two projections $\\Pi_{r}$ and $\\Pi_{\\Theta}$. Our proposed LPGD algorithm performs the convex projection $\\Pi_{\\Theta}$ first and then applies the low-rank projection $\\Pi_{r}$. The key inequality we derive in the proof of Thm. C.2 is (eq. (60) in the submission)\n$$\n\\lVert \\mathbf{Z}_\\text{$t$} - \\mathbf{Z}^{\\star}  \\rVert_F  \\le 2\\eta  \\lVert \\mathbf{Z} _{t-1} - \\mathbf{Z}^{\\star}\\rVert_F + \\lVert \\Pi_t (\\tau \\Delta _{\\Theta} \\mathbf{Z}^{\\star} )  \\rVert_F,\n$$\n\n\n\nwhere $\\tau \\Delta_\\text{$\\Theta$}\\mathbf{Z}^{\\star}:=\\mathbf{Z}^{\\star} - \\Pi_\\text{$\\Theta$}(\\mathbf{Z}^{\\star}-\\tau \\nabla f(\\mathbf{Z}^{\\star}))$ denotes the gradient mapping at $\\mathbf{Z}^{\\star}$ w.r.t. the convex constraint $\\Theta$, and $\\Pi_\\text{$t$}$ is a linear projection onto a $3r$-dimensional linear subspace that depends $\\mathbf{Z}^{\\star}$, $\\mathbf{Z} _{t}$, and $\\mathbf{Z} _\\text{$t-1$}$. The last error term above can be bounded above uniformly in $t$ using $\\lVert \\Pi _\\text{$t$}( A)\\rVert _\\text{$F$}\\le \\sqrt{3r} \\lVert A \\rVert _\\text{2}$. So we can apply the above inequality recursively to obtain the desired result. \n\nNow if we consider an alternative algorithm that uses the  low-rank projections $\\Pi_{r}$  first and then the convex projection $\\Pi_{\\Theta}$, then we can derive a corresponding key inequality: \n$$\n\t\t\\hspace{-2cm}(*)\\cdots  \\qquad \\qquad \t\\lVert \\mathbf{Z} _\\text{$t$} - \\mathbf{Z}^{\\star}  \\rVert _\\text{$F$}  \\le  2\\eta \\, \\lVert  \\mathbf{Z} _\\text{$t-1$} - \\mathbf{Z}^{\\star}\\rVert _\\text{$F$}+ \\lVert  \\tau\\Delta^{t} \\mathbf{Z}^{\\star} \\rVert _\\text{$F$},\n$$\n       \nwhere $\\tau \\Delta^{t}\\mathbf{Z}^{\\star}:=\\mathbf{Z}^{\\star} - \\Pi _{t}(\\mathbf{Z}^{\\star}-\\tau \\nabla f(\\mathbf{Z}^{\\star}))$ denotes the gradient mapping at $\\mathbf{Z}^{\\star}$ w.r.t. the *virtual* linear constraint that we constructed during the proof to approximate the low-rank constraint. Indeed, the above inequality can be derived by modifying the argument in (42)-(47) in the submission assuming the reverse order of projections. We omit the details due to the space constraint. \n        \n  \n        \nThen by recursively applying the inequality $(*)$, we can obtain \n$$\n\\lVert \\mathbf{Z} _\\text{$t$} - \\mathbf{Z}^{\\star}  \\rVert _\\text{$F$}  \\le  (2\\eta)^{t}   \\lVert  \\mathbf{Z} _\\text{$0$} - \\mathbf{Z}^{\\star}\\rVert _\\text{$F$}  +  \\sum _\\text{$k=1$} ^{t} (2\\eta)^{t-k} \\lVert  \\tau\\Delta^{k} \\mathbf{Z}^{\\star} \\rVert _\\text{$F$}.\n$$\n       \nHence the rate of convergence we would get is the same as the original algorithm, but the additive error takes a different form. Since the \"low-rank gradient mapping\" $\\Delta^{k}\\mathbf{Z}^{\\star}$ depends on the iterates $\\mathbf{Z} _\\text{$k$}, \\mathbf{Z} _\\text{$k-1$}$, we find it easier to control the gradient mapping with respect to the convex projection that comes out from the analysis of the original algorithm. \n        \nWe absolutely agree with the reviewer that this discussion will be helpful for the readers. We will add this discussion as a remark in the revision. \n\n\n\n> Question: Will a distribution free analysis in section 4 be more interesting ? With feasibility set being bounded and the objective in terms of the data is lips.conts., this should be possible.\n\n**Response:** This is an excellent suggestion. The only distributional assumption we have made in Section 4 is that the noise terms $\\epsilon_{i}$ and $\\epsilon_{i}'$ follow Gaussian distributions with a mean of zero. In fact, our current analysis holds when assuming only a sub-Gaussian distribution, as the only place we used Gaussian distributional assumption was in the concentration inequalities in Lemmas D.2 and D.3 (in order to obtain tail bounds in eq. (162) and (174)-(175)), which are already stated for the sub-Gaussian case. If we restrict the support of noise to be bounded, then similar concentration inequalities hold without any distribution assumption (e.g., by matrix Burnstein inequality). Such bounded noise assumption might be reasonable when the feasibility set is bounded, as the reviewer pointed out. We will add this discussion in the revision."}, {"review_id": "66DJIallqp", "reviewer": "Reviewer_ucPD", "summary": "This paper proposed a novel supervised dictionary model, with two variations, one feature-based and one filter-based. The problem setting is on classification tasks with both high- and low-dimensional feautures, where the high dimensional features are learnd through dictionary learning, and then intergrated in to multinomial logistic regression. \n\nThe problem estimation is formulated to use a low-rank projected gradient descent algorithm. Exponential convergence guarantee is provided for the proposed algorithm. ", "questions": "A few typos:\n\n-   Line 151, \"is scalar and continuous\"\n-   Line 156, the equation should be the one mentioned on line 154 instead of (6).\n\n", "limitations": "The authors have sufficiently addressed the limitations. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper, overall, is written clearly and it is a pleasure to read. \n\n1.  The proposed algorithm provides theoretical guarantees in terms of convergence.\n2.  The proposed algorithm demonstrates improved performance compared to prior algorithms, as well as selected standard classification algorithms. It additionally exhibits good interpretability.", "weaknesses": "It would be beneficial to conduct benchmarking against deep neural networks in order to gain insights into gaps, if any, in performance and to understand potential trade-offs between interpretability and classification accuracy.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your positive feedback on our submission. The reviewer has provided the following comments: \n\n**`Q1. It would be beneficial to conduct benchmarking against deep neural networks in order to gain insights into gaps, if any, in performance and to understand potential trade-offs between interpretability and classification accuracy.`**\n\n**Response**: Thank you very much for the suggestion. In the revision, we have included an additional benchmarking analysis involving deep convolutional neural networks (CNN) and feedforward neural networks (FFNN) to provide deeper insights into the performance of our proposed methods. (We report this result in the supplementary 1-page author response PDF). \n       \nFor the task of classifying microarray data into cancer classes, we compared the performance of our method with both CNN and FFNN. Specifically, the CNN architecture was designed with a convolutional layer with 32 filters and a kernel size of 3, followed by an average pooling layer with a pool size of 2. Subsequently, a second convolutional layer with 64 filters and a kernel size of 3 was integrated, further followed by another average pooling layer with the same pool size. The architecture was finalized with a flatten layer, a fully connected layer of 128 neurons activated by ReLU, a dropout layer with a rate of 0.5, and finally a fully connected layer with a sigmoid activation. On the other hand, the FFNN consists of a fully connected layer featuring 64 neurons with ReLU activation, followed by a dropout layer with a regularization rate of 0.5. A subsequent fully connected layer with 32 neurons activated by the ReLU was incorporated, followed by a fully connected layer with a sigmoid function. This comparative analysis was repeated five times, consistent with the procedure outlined in the main paper. \n      \nIn our experiment,  the CNN achieved an average accuracy of 0.769 (0.07) (here 0.07 is the standard deviation, and we use this notation hereafter) on the pancreatic cancer dataset, and an average accuracy of 0.854 (0.06) on the breast cancer dataset. In contrast, the FFNN yielded an average accuracy of 0.816 (0.04) for pancreatic cancer and 0.890 (0.02) for breast cancer. We also included the revised table with all methods in the one-page supplementary PDF for rebuttal so that the reviewer can check.\n\nAn intriguing observation emerges from our benchmarking analysis. While the FFNN's performance on the breast cancer dataset is comparable to ours (The LPGD algorithms for SDL), the overall performance of CNN is notably inferior to ours. This disparity can primarily be attributed to the small sample size of the training set (145 samples for breast cancer and 25 samples for pancreatic cancer) in comparison to the substantial dimensionality of gene features (exceeding 30,000 features). We note that obtaining a substantial volume of biomedical data for cancer research is very expensive, making it challenging to feasibly train complex models such as deep neural networks. The significance of our approach becomes evident in its ability to retain robust performance even when facing the challenges posed by a restricted sample size and a complex high-dimensional feature landscape. Moreover, our method augments this resilience with the advantage of interpretability. We will include the comparison with CNN and FFNN, as well as the above discussion, in the revision.\n\n**`Q2. A few typos: Line 151, \"is scalar and continuous\"`**\n\n**Response**: Thank you for this comment. We fixed this in the revision. \n\n**`Q3. Line 156, the equation should be the one mentioned on line 154 instead of (6).\"`**\n\n**Response**: Thank you for this comment. We fixed this in the revision. "}]

Author rebuttal:
We submit an optional 1-page PDF to show the revised Figures 2 and 3 with captions. In Figure 2, we present a comprehensive comparison between our LPGD algorithm and BCD. In Figure 3, we include additional experimental details, focusing on breast cancer classification, which successfully identifies well-known oncogene and cancer-associated genes (prognostic markers). Moreover, Figure 3 now contains further benchmarks that compare our approach to convolutional neural networks and feed-forward neural networks.

Meta-review: The reviewers commend the paper's approach in DL, which combines supervision with reconstruction to facilitate the learning of a low-rank matrix. The transformation of this approach into a low-rank optimization problem with a convex objective is particularly noteworthy. The inclusion of novel theoretical results is also acknowledged.

Certain reviewers initially expressed concerns about the paper's focus in relation to the extensive studies within DL, particularly those emphasizing overcomplete dictionaries with sparse regularization. However, these concerns, along with others, were effectively addressed during the authors' engaged response.

Taking into account the recommendations from the reviewers and the productive discussions during the rebuttal, I am inclined to recommend its acceptance.

For the revised version, I strongly urge the authors to fully integrate the constructive suggestions provided by the reviewers. This should encompass clarifying the distinctions between their approach and the dominant dictionary learning paradigm involving sparse regularization on overcomplete dictionaries. Additionally, revising the terminology from "supervised dictionary learning" to "supervised low-rank matrix factorization" would enhance clarity. It would be beneficial to make the source code available and to augment the paper with more comprehensive experimental results and comparisons.
Decision: Accept (poster)