Rethinking PCA Through Duality
Authors: Jan Quan
Johan Suykens
Panagiotis Patrinos
Venue: NeurIPS 2025 poster ()
Keywords: PCA
DCA
Non-Convex Optimization
Unsupervised Learning
Primary area: general_machine_learning

Abstract: Motivated by the recently shown connection between self-attention and (kernel) principal component analysis (PCA), we revisit the fundamentals of PCA. Using the difference-of-convex (DC) framework, we present several novel formulations and provide new theoretical insights. In particular, we show the kernelizability and out-of-sample applicability for a PCA-like family of problems. Moreover, we uncover that simultaneous iteration, which is connected to the classical QR algorithm, is an instance of the difference-of-convex algorithm (DCA), offering an optimization perspective on this longstanding method. Further, we describe new algorithms for PCA and empirically compare them with state-of-the-art methods. Lastly, we introduce a kernelizable dual formulation for a robust variant of PCA that minimizes the $l_1$-deviation of the reconstruction errors.

Rethinking PCA Through Duality
Jan Quan' ESAT-STADIUS Leuven.AI KU Leuven, Belgium jan www.nuuuuuuueeen.n....
Johan Suykens
Panagiotis Patrinos
ESAT-STADIUS & Leuven.AI
ESAT-STADIUS & Leuven.Al
KU Leuven, Belgium
KU Leuven, Belgium
johan www...nrwwnnkulleennnnn.no be
p pat.sppssetnoeeeeeeeeee. be
Abstract
Motivated by the recently shown connection between self-attention and (kernel) principal component analysis (PCA), we revisit the fundamentals of PCA. Using the differencceof-cconvex (DC) framework, present several novel formulations and provide new theoretical insights. In particular, we show the kerneiizability and out- of- of-sample applicability for PCA-like family problems. Moreover, we uncover that simultaneous iteration, which is connected to the classical QR algorithm, is an instance of the differenceeof-convee algorithm (DCA), offering an optimization perspective this longstanding method. Further, we describe new algorithms for PCA and empirically compare them with state-of-the-ard methods. Lastly, introduce kernelizable dual formulation for robust variant of PCA that minimizes the lı-deviation the reconstruction errors.
Introduction
Principal component analysis (PCA) is one of the most fundamental dimensionality reduction techniques and has seen widespread use in various fields, including machine learning [20, 10, 38, 23], signal processing [4, 24], and computer vision [52, 78, 82], due its ability to efficiently extract meaningful features from high-ding -igg-dimensional data. Its nonlinear extension, kernel PCA [51], extends this method through feature mapping into high-dimensional space (possibly enabling the discovery of nonlinear patterns. Crucially, this is achieved without ever computing the transformed data explicitly, using the kernel trick [11].
In recent years, (kernel) PCA has regained interest due to its connection to the self-attention mecha- nism [79]. which lies the heart of the widely successful transformer architecture. Transformers have revolutionized a wide range of applications such as natural language processing [13, 6, 42], computer vision [25, 59, 49, 76], and speech processing [41, 50]. This success has led to many attempts to understand the underlying mechanisms from more principled viewpoint.
In particular, [69] establishes direct link beween self-attention and kernel PCA by showing that the attention outputs are projections of the query vectors onto the principal components axes of the key matrix in feature space. In similar vein, [18] recovers self-attention from asymmetric kernel singular value decomposition [67], which can seen as generalization of kernel PCA to multiple data sources. Moreover, it has also been shown that models based low-rank matrix decompositions can achieve results on par with modern transformers when trained correctly [32].
"Corresponding author.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Motivated by these modern connections to PCA, we aim to revisit the fundamentals of PCA. Our main tool is the concept of difference-of-conver (DC) duality, which provides general framework for analyzing nonconvex problems, as it is for example known that any C2 function admits difference-of-conver formulation [34]. One parameter-free algorithm that suited for solving these types of problems is the difference-of-convee algorithm (DCA) [68], of which the convex-concave procedure [81] is a special case. This algorithm, which fits into the injorizatiin-minimzzaln framework [66], has found many successful applications in machine learning including expectation- maximization [2I], successive linear approximation [12], and various clustering methods [3, 56].
Contributions Our contribution is threefold.
1. We derive three novel DC dual pairs for PCA. For these formulations, we compare some simple gradient methods with existing solvers. Moreover, we show that if one of the functions in the primal is unitarily invariant, its DC dual is kernelizable and out-of-sample applicable.
2. We show that simultaneous iteration [77, Algorithm 28.3] is an instance of DCA applied to the variance maximization objective, which in turn is related to the famous QR algorithm for computing eigenvalues of dense matrices. This result gives a new connection between numerical linear algebra and optimization, akin to the connection between (linear) conjugate gradients [35] and accelerated gradient descent.
3. Based on least absolute deviation formulation for the robust subspace recovery problem, we provide a novel DC formulation that is kemelizable. Moreover, we show that DCA is related to approaches based on iteratively reweighted least squares for this problem.
Related work The idea of using DC duality on the variance maximization formulation of PCA is not new. In particular, it has been studied for one component in [9, 36, 5J while more recent work extended it to multiple components [74]. Nevertheless, we propose three DC pairs that are based on other formulations of PCA and have, to the best of our knowledge, not been considered before.
The QR algorithm [31] is a classical algorithm for computing the cigenvalues of dense matrices. It known as one of the top ten algorithms of the 20th century [19] and still the de facto industry standard in high-performance computing libraries and software such as MATLAB. The connection between the power method without normalization and DCA is described in [71, Prop. 4]. Notably, this connection is only made for the leading principal component. In contrast, we show that DCA applied to the variance maximization objective of PCA yields a fundamental link with simultaneous iteration [77, Algorithm 28.3], which in turn is related to the QR algorithm. Some more advanced algorithms for solving the same formulation can be found in [2, 64, 1, 27].
Beyond standard PCA, there has been significant interest in incorporating additional structure such as sparsity, or enhancing robustness to noise and outliers [15, 53, 72, 40, 28, 29, 71, 48]. In the kernel setting, [74] extends DC duality to robust PCA by considering other variance-like objectives, though this complicates interpretation and the out-of-sample extension was not considered.
Preliminaries
2.1 Notation
We denote by () the Euclidean inner product for vectors and the Frobenius inner product for matrices. Let T E Rmxx. We denote by Te, the ith row of T. The Schatten p-norm of T with E +o) is defined through IIT'I', = Emin Xminn (T, where (T') E mimi((..)) denotes the vector of singular values of T in nonincreasing order. Szo is defined as the spectral norm, i.e., the largest singular value of its argument. The Schatten 1-norm is also known as the nuclear norm or trace norm. The Schatten 2-norm is also known as the Frobenius norm or the Hilbert-Schmin norm. (full) SVD of T is denoted as U (Tag(a(T))V where Diag(o(T)) E Rm is a rectangular diagonal matrix and U. V are real orthogonal. If T is square, X(T) denotes the eigenvalues of T in any order. We use Diag (X(T)
denote a square diagonal matrix with X(T) on the diagonal. The (closed) unit ball of a norm ||- || on X' is denoted by B1 := {z E X |||||| S 1}. The extended reals
denoted by E. = RU (o Let f : I + R. The convex subdifferential of f is denoted by

---

0f. whereas its convex conjugate is denoted by f" The function f is called absolutely symmetric if f(7) = f()) for all y and where Ý denotes the vector with components [y| in decreasing order. function F: IR + R is unitarily invariant if F(VXU) F(X) for all X e mR and both U E R V E m are real orthogonal.
Remark 2.1. For clarity of exposition, we perform complete discussion in Euclidean spaces. The extension to nnfinte-dimensional Hilbert spaces for kernel methods only needs some additional
technical details. More concretely, let R -> H be some feature mapping and {2;] i=1 C Rd
the given data. Any time formulation contains the data matrix X (z]
2 [N) it can be formally replaced by I : H + R where (Tw); = (p(az;) w) for all w E H. and (,, 1) denotes the inner product of the Hilbert space. The kernel matrix is then We also employ the following informal definition to denote formulations which allow for practical implementations.
DT = (((()))) p(ay) )
Definition 2.2 ((Informal) kernelizability). Let X E IR be the data matrix. A problem formulation is said to be kernelizable if it can be written solely in terms of the kernel matrix K = XX" E XNN a problem kernelizable, we moreover call it sut-oe--ample applicable for a new data point á E RR, the "output' of the problem can be computed from only K and X
2.2 Difference-of-convee duality
Our primary tool for deriving new formulations is difference-of-conwe (DC) duality, also known as Toland duality, which yields a pair of primal-dual problems between which strong duality holds.
Tgdx* + R. F Nx - T. be two convex, closed and proper functions, and X Rdxe N x a linear mapping. Then, the following pair of primal-dual problems
Proposition 2.3 (Toland duality [73D. Let G
minimize G(W) - F(XW) WER
(P)
minimize (H) - GX HERN
(D)
have the same infimum, i.e., strong duality holds. If W solution of (P), then any H* OF(XW*) is solution of (D). If H* is a solution of (D), then any W* E 0G (X H*) is solution of (P).
Since this DC duality is a special case of the duality described in [62, 11H). the way the primal problem is perturbed/modified can lead to vastly different dual problems, as will become clear from comparing DC dual pairs (1)(m) and (n)-(o) in Figure 2.
2.3 Difference-of-convey algorithm
effective and widely used algorithm solve DC problems of the form (P) is the difference-of- convex algorithm (DCA) [68]. The iterates of DCA are given by
(h)) (F(XW(())); w(++1) E OG (X H())
arg min G(W) (W, X H(*)). (DCA)
WeR4xa
To make sure the iterates are well-defined, we also require the following constraint qualification. Assumption 2.4. X dom(aG) S dom(OF) and X range(OF) S. rang
We have the following relation between the iterates of applying (DCA) to (P) and (D) respectively. Proposition 2.5. Let (W()) Y be the sequence of iterates obtained from applying (DCA) (P) starting from some W(O) E Bdxx Then, there exists a sequence of iterates (H( [== obtained from applying (DCA) to (D) starting from (V) E (F(X)(() such that W(k+1) E OG* (XT (k)).
The convergence properties for (DCA) are well-established in the literature [68, 63, 55]. We sum- marize the basic results for completeness. Consider the problem (P) and suppose that the iterates [W o are generated by (DCA).
Since (DCA) fits into the miniiiziiioommiiiiiztttooooo framework [66], it is a descent method (without linesearch), i.c., G(W((++1) - ((XW(++1) < G()(+)) F(XW(+))
3

---

If G(W(k+1)
F( then (DCA) terminates at the kth iteration and ((() is a critical point of (P). Here, critical point W of (P) is defined as point satisfying (G(W) X (F((W) o Ó, which is a necessary condition for optimality. If G(W) F(XW) is bounded below, then every limit point of the sequence (W(()) k= is critical point. We note that the optimal value being bounded below is an assumption that is trivially satisfied for all our upcoming PCA formulations.
= G(W(k)) F(XW(++1)
In general, under mild assumptions, the rate for (DCA) is sublinear. However, if certain growth conditions hold, then linear rates are possible, see for example [55, 43].
3 Principal Component Analysis
In this section, we briefly review the classical formulations of PCA before proposing novel DC dual pairs. In particular, we show that PCA-like family with one of the functions the primal being unitarily invariant, is kemelizable and out-of-sample applicable in the dual. Moreover, we look into (DCA) applied to these new formulations and provide a connection to simultaneous iteration.
3.1 Classical PCA
There are many formulations of PCA that are commonly used. Some of the most important ones are shown in Figure 1. Following the exposition of [54], we start from the fundamental low- rank approximation of the data matrix X shown (a). By the classical Eckart-
theorem [26], the optimal solution is obtained by the sum of scaled outer products of left and right singular vectors, corresponding to the top & singular values. This solution unique if the spectrum of the covariance matrix is simple.
By parameterizing A in (a) its (full) SVD EAVA where UA and VA are orthogonal matrices, one obtains the formulations (b) and (c) by taking W = VA and H = UA respectively. It should be noted that these problems do not admit unique solutions. In fact, there are infinitely many non-isolated minimizers since the problem is invariant to orthogonal transformations. This formulation is closely related the -urer-Monteiro factorization for semidefinite programming [14, 47], though one key difference is that we impose orthogonality constraints on one of the factors.
To obtain formulation (c) from (b), it suffices to remark that (b) is a linear least-squares problem in the unconstrained factor B and the unique solution is B X W. The resulting formulation minimizes the reconstruction error of PCA linear autoencoder [58]. The formulation (f) is obtained from (e) in a completely analogous manner. By expanding the Frobenius norm in (c), the formulation (d) is obtained. This is the classical variance maximization objective of PCA. where is common, though not necessary, to assume that each column has zero mean. The formulation (g) is similarly derived from (f) and now contains the kernel matrix X.X (also known as the Gram matrix of the data). The relation between (d) and (g) lies at the heart of the success of many kernel methods.
Formulation (h) is derived by considering the low-rank approximation of the covariance matrix X X, which closely related to the formulations from the first column in Figure 1. To obtain (i), we impose the additional constraint that A is positive semi-definite in (h) since it does change the solution. A SVD A can now be written as (U^>)/3)(U) (U which indeed yields the desired formulation. The last two formulations (j) and (k) follow immediately by expanding the square and using the definition of the Schatten 4-norm. should be noted that four more formulations can be derived in an analogous manner by starting from the low-rank approximation of the kernel matrix X instead.
3.2 DC dual formulations
This subsection is devoted to describing novel DC dual formulations of PCA. Our main result is the following.
Theorem 3.1 (Fundamental PCA DC pairs). Let X Nxd be the data matrix and S rank(X). Then, the three DC dual pairs in Figure 2 hold. Moreover, (1) and (n) share the same minimizers, as do (o) and (q).

---

(a)
(h)
minimize AER MERM*E,rank(AA rrank(A)= =)==
| All32
minimize AERENÉ (ankk(A) -)
1|XTXX - -ll
= 1 DRUNN
A" = Σi= afuio
(b)
(c)
minimize
X| BW 113
minimize
X HC 13
BER BERMx,.WER WER"X
HER .CERt 11=,
W W=1,
minimize
I| X WW 133
B' = (aiu) ... orUL)
=1== (ui ... us)
WER*X
W (al0)
W* = (01 va)
C (010)
$
(c)
(f)
ii)
minimize
WER*X
|| - XWW 11
minimize
HERNK* H-1.
|| - HH x xl3,
minimize
WgR-*+ 2
1
WWW 11 X
W=L.
W (0101
o.v.)
W
(0)
v.)
H
(ui
u.)
$
(d)
(g)
(k)
minimize WER
minimize #g
x
m
minimize WE
[
2 W (aiv) a.v.)
W' W==,
H -=1,
W (0)
02)
H =- (w1
u.)
Figure 1: Classical PCA formulations for data matrix X N and rank(X) principal components. The starred variables denote (not necessarily unique) global minimizers of the associated formulation. A (full) SVD of X is given by U Diag(o)V The arrows <> denote that the two formulations are 'equivalent', in the sense that one can easily obtain the minimizers of one problem by solving the other.
(n)
(p)
minimize We [W]laws] $1
XWI
minimize WE
||XW||s
minimize WERlx
||
|W <1
DC
I dual
DC 1 dual
DC dual
(m)
(0) minimize "ER [Rls,,
(q)
j1m TH|s
||X T
ls.
minimize HeR S [ 18151
X H|]
minimize HERRXA
Figure 2: Three fundamental DC dual pairs for PCA.
To characterize the global minimizers of each problem from Figure 2, it suffices to combine the following known result with Proposition 2.3.
Proposition 3.2. Let X N be the data matrix and S rank(X). Then, the optimal value of (l) is -1 Die (X)2. Consider any eigenvalue decomposition X X WAWT, where A = Diag(A1 ...,,,)) with X 2 A2 ... 2 Xd and W is a real orthogonal matrix. Then, the matrix W E I consisting of the first columns of W achieves the optimum of (l).
The DC dual pair (1)(m) well known for one component [9, 36, 5], and has recently been extended to multiple components in [74, Eq. (1) and (6)]. It should be noted that (1) is equivalent to (d) by
5

---

observing that the variance maximization formulation maximizes (non-constant) convex function such that the constraint set can be relaxed its convex hull, since the solutions necessarily lie on the boundary [61, Cor. 32.3.2]. Moreover, the constraint set is identified as the Stiefel manifold and its convex hull is the closed unit ball of the spectral norm [39, $3.4] which is exactly what appears in (l).
The two remaining DC dual pairs (n)-(o) and (p)-(q) are, to the best of our knowledge, novel. Their derivations rely on the simple fact that squaring/taking the square root and (positive) scaling does not change the solution sets of the minimization problems. Nevertheless, since the conjugates of the resulting transformed problems are quite different, so are their DC duals. In particular, the DC dual pair (n)-(o) is arguably more fundamental due to its symmetry. Each formulation captures a specific trade-off between Schatten p-norms.
Another interesting and novel dual problem follows from starting from (k) in Figure 1. Proposition 3.3. Let X E IR be the data matrix and < rank(X). Then, the following DC dual pair holds.
minimize W
|W
[XWISs 2
DC dual
minimize HER x* 2
|
3 |
H|l
4/3 Saya"
4
Remark 3.4. By replacing X with X and vice versa (cf. (d) and (g)), four more analogous DC dual pairs can be formulated as the ones from Theorem 3.1 and Proposition 3.3.
In the Hilbert space setting, formulations involving the kernel matrix K == X are needed, as they enable computations via the kernel trick. This is the case when G the primal unitarily invariant as we show in the following proposition.
Proposition 3.5 (Kernelizability for unitarily invariant G). Let X E Nxd be the data matrix and S rank Consider the primal problem (P) and suppose G = 9 o is unitarily invariant. Then, the corresponding DC dual (D) is kernelizable. particular, (D) can be written as
(*) -g" (VALHT KH)
minimize HERN x
where v is taken elementwise and X()) denotes the vector of eigenvalues of its argument (in any order).
Remark 3.6. The decomposition G g a always exists for unitarily invariant functions, where moreover is absolutely symmetric, see Fact A.4. This fact also ensures that g" does not depend on the order of its arguments.
Remark 3.7. Note that while the objective function contains an eigenvalue decomposition, is of an X matrix, which is small in practice, as typically only a few principal components are required. In particular, when s = 1, we have that A(H KH) = HH" KH e R.
The previous result allows us to characterize wide variety formulations that are kernelizable the DC dual, making the computations tractable. In particular, all the previously described formulations are kernelizable, since they are based on Schatten p-norms, which are unitarily invariant.
Another important requirement in practical upplications is handling new, unseen data. The following result shows how one can perform feature extraction of new sample by projecting onto the associated primal fower-dimensionas subspace using only operations that are compatible with the kernel trick, thereby yielding out-of-sample applicability.
Theorem 3.8 (fut-of-sample applicability for unitarily invariant G). Let X E N be the data NNN the kernel matrix, & S (ank(X) the number of principal components
matrix, K XXT T E
and ã E RRd a (new) data point. Consider the primal problem (P) and suppose G g ooo a is unitarily invariant, Suppose moreover that H E T2 Nxx is minimizer of the corresponding DC dual (D) and H*T T KH* is nonsingular. Denote the eigenvalue decomposition of H* K E Raxs as V DXag Then, there exists W* which a solution of (P) such that
W & VDiag(u X X H² X
where E Og" (V)), @ denotes the elementwise product, and both VA and 1-12 are to be understood elementwise.
Remark 3.9. In practical implementation, after finding an optimal solution H² based on the training data, the matrix VDiag(u C T H* e R only needs to be calculated once.

---

3.3 DC algorithms
We now derive the algorithm (DCA) for multiple formulations from Subsection 3.2. To this end, subgradient of a unitarily invariant function G = g o a is needed, which we can calculate from the following known result.
Proposition 3.10 ([46, Cor. 2.5]). Let X Rmxx and let g Rmin((m.)) R be proper and absolutely symmetric function.. Suppose moreover that o o S domg. Then,
(x)) = (iag(u)) I Og(o(()), U Diiag(o(X)) -SV((X))
By now applying (DCA) to the formulations (d) and (g) with the constraint set relaxed to its convex hull, one obtains Algorithms and Here SVD denotes a compact SVD, i.e. for [U, Σ VT] SVD(X) with X E Rmxn, we have U RR and VT E Bexx where denotes the rank of X, and U I, VTV. The derivation of these algorithms, as well as the reason why a compact SVD can be used here, is detailed in Appendix C.
Algorithm (DCA) for (d)
Algorithm (DCA) for (g)
Require: Matrix X Nxd
Require: Matrix K = XXT E RNN
I: Initialize (( E Rdxa s.t. || =
1: Initialize ( E R Nx s.t. |s =
2: Initialize & 0
Initialize ls +
3: repeat
3: repeat
4: [U(K), E0,) ((K) 1 SVD(X (VD(XTXW(1))
(U(K), g(k), V()) (+
5: W(++T) + v(k)
5: l(k+1) +- ((A)
6: k+
6: k- k+1
7: until convergence
until convergence
Observe that when 8 1, the SVD is applied to a vector, resulting in its normalization, and the iterations correspond to the power method. This fact was also observed in [71, Prop. 4] after removing the sparsity constraint. In the next result, we show that these algorithms are connected to simultaneous iteration [77, Algorithm 28.3] which is strongly related to the QR algorithm [77, Theorem 28.3].
Theorem 3.11.
Algorithms and correspond to simultaneous iteration applied to X X and X respectively, up orthogonal transformation.
Although the previous result states that (DCA) identifies the same subspace as the simultaneous iteration method, whose iterates span the leading eigenvector basis, (DCA) does not necessarily converge to the eigenvectors themselves, as the subspace admits infinitely many (orthonormal) bases. For the purpose of obtaining lower-dimensional representation with meaningful features, the projection onto this subspace is often sufficient. If the principal vectors are specifically required, they can be recovered by performing an eigenvalue decomposition of smaller X matrix.
An immediate corollary of this result and the known convergence result for simultaneous iteration [77, Theorem 28.4] is that Algorithms and 2 inherit the linear rate well the convergence to global optimality with probability one, starting from random initialization. Similarly, the proximal gradient method that
in Section 5 has well-known interpretation of (DCA) applied to (G + (1/2)ll I12) (F + [1/ I1?), where G denotes the indicator of a convex set. Therefore, identity-shifted covariance/kernel matrix (i.e., naturally incorporates regularization), such that again convergence to global optima is guaranteed.
these can be connected to simultaneous iteration
Applying (DCA) to the problems (m) through (q) as well as their transposed variants (cf. Remark 3.4)
deferred to Appendix C. These algorithms are cquivalent up to normalization and simple operations, as is to be expected due to the elementary problem transformation. A more interesting algorithm follows from considering (DCA) for the problems from Proposition 3.3. These are presented in Algorithms 3 and 4. The iteration is similar but now also uses information from the singular values through some nonlinear preconditioning. Here EIG denotes an eigenvalue decomposition of its argument but only keeping the columns corresponding to nonzero eigenvalues. Note that we are able to obtain kernelizable formulation in the dual. The derivation is deferred to Appendix C.
7

---

Algorithm 3 (DCA) for Proposition 3.3 (primal) Algorithm (DCA) for Proposition 3.3 (dual)
Require: Matrix X RR
Require: Matrix K XX" R
1: Initialize (( E Rdxx
1: Initialize H E Nx
2: Initialize &
Initialize k +-
3: repeat
3: repeat
4: [U(A), 2X1), v(k t SVD(X xw()) SVD(XTxw())
[V(#), A (K) (HG(H KH(M)
5: W(++)) U 1/8 ((()
H(++1) < (A
6: ++++1
6: k
7:
until convergence
7: until convergence
Robust PCA
In this section, we start from formulation (c), the reconstruction error interpretation of PCA linear autoencoder, and note that it can be equivalently written as
minimize WER
Σ I WW x1,1l
WnT,
where || denotes the Euclidean norm. is well known that measuring the error through the squared loss leads to non-robust models. A popular approach to obtain robust version is to remove the square, thereby considering l penalty on the reconstruction errors instead. This known in the robust subspace recovery literature as the least absolute deviation formulation [45]. We have the following DC dual pair, which was considered in [9] for only one component.
Proposition 4.1. Let X N xd be the data matrix and S &rank(X). Then, the following DC dual pair holds:
minimize WER
Σ IX.: - WW x
minimize HERN Σ IIx..Ill11 || - -11X H|lss
[[W]s Poc 51
The dual formulation kernelizable, as it can be written as
NV Σ CV(A (i, + [|4...19)) (1 - tr (VHTXXTH)
minimize HERNK X
Moreover, since the associated G' is the indicator function of the spectral unit nomn ball, which is unitarily invariant, it satisfies the assumptions of Theorem 3.8 and therefore the problem out- of-sample applicable. By applying (DCA) to both formulations from Proposition 4.1, obtain Algorithms 5 and 6. In Appendix D, we show how the primal algorithm related to an iteratively reweighted least squares scheme.
Algorithm 5 (DCA) for Proposition 4.1 (primal)
Algorithm 6 (DCA) for Proposition 4.1 (dual)
Require: Matrix X E Nxd
Require: Matrix K =XXT NNN
1: Initialize (0)) E IRdx s.t. || Iss..
Initialize ()) E NR
2: Initialize k + 0
Initialize 4-0
3: repeat
3: repeat
4: for i
1 to N do
4: [V(K), A(0 + EIG (H()))
(XW(2)) bi
5: Y
(R)
(*) (H(M)((N)) (A()))11/(((
VIIX,. - (((
6: for i = 1 to N do
6: end for
(U()) {()) ( SVD
(((k)
7:
VK,, (I(Y
W(++1) (())
8:
8: end for
9: kt k+1
9: k+ *+k+1
until convergence
10: until convergence

---

Table 1: Timing results for various methods applied PCA formulations. The problem setting
(N, d, S, e) denotes data matrix X E INxd with entries sampled from standard normal distribution. s denotes the computed number of principal components and & the stopping criterion tolerance. All timings are in milliseconds, and timings longer than 5 seconds are not displayed.
Method
(4000, 2000, 20, 10 30
(2000, 4000, 20, 10 10
(4500, 4500, 20, 10
ZeroFPR (1)
2828. + 376.5
1949.3 + 97.9
8719.1
ZeroFPR (n)
210.4 + 54.9
196.3 49.1
441.0 - 86.3
ZeroFPR (o)
392.6 + 100.7
463.0 + 121.5
1130.7 7 = 143.3
ZeroFPR (q)
2341.7 + 310.4
2021.0 + 205.4
PG (1)
PG (n)
207.5 + 14.0
139.1
42.4
309.1 26.1
PG (o)
203.6 + 57.4
216.4 + 29.1
689.7 = 53.4
PG (
Algorithm
640.1 + 154.3
2630.9
33.0
3696.4 115.2 3294.2 = 359.8
Algorithm 2
2218.5 + 35.2
635.3 + 85.1
Algorithm 3
3848.8 + 393.3
Algorithm 4
1868.2 + 182.6
SVDS
808.0 0 + 37.6
808.2
29.6
2985.8 128.3
KrylovKit
1130.4 + 21.6
940.1 + 13.6
3844.6 + 253.2
10s
Time (ms)
Time (ms)
10a
ZF (n)
ZF (n)
10a
PG (n)
10a
PG (n)
SYDS
SVDS
50
100
150
200
0
500 1,000 1,500 2,000 2,500
d
Figure 3: (left) Timings for the problem setting (2000, 2000, 8, 10-3) as defined in Table 1, with vary- ing and one standard deviation error bars. (right) Timings for the problem setting (2000, d, 30, 10 10-3) as defined in Table with varying d and one standard deviation error bars. In both figures, ZF is shorthand for ZeroFPR.
Experiments
We compare several simple gradient- based methods various formulations from Theorem 3.1 and Proposition 3.3. All experiments are implemented in Julia 1.11.1 on a machine with an AMD Ryzen 7 Pro 5850U processor and 32 GB RAM. The timings taken using Benchmark [17]. Our results are presented in Table 1. The code for reproducing all experiments is publicly available?.
For the constrained problems, apply both the proximal gradient (PG) method and ZeroFPR [70] using the implementation from [65]. Regarding the convergence properties these algorithms, the standard assumptions to minimize objectives of the form f(x) + ((x) using these methods require f to be L-smooth. In general, the (global) rates of PG and ZeroFPR to stationary points are both sublinear [70]. We note that, strictly speaking, for some formulations not pose a problem since the objectives are concave. Therefore. the Euclidean descent lemma [8, Lemma 5.7] holds with / 0, and both PG and ZeroFPR enjoy the same theoretical guarantees as described above. In fact, we may choose arbitrarily large stepsizes, despite nonsmoothness. As baseline, compare formulations against classical eigensolvers. Specifically, used the classical ARPACK implementation of SVDS, which is based on the implicitly restarted Arnoldi
is nonsmooth. However, this does
hnttps://githut http:///tthhuu..m////////// wwtp////////aaduuuuaaaao duality

---

method [44]. Moreover, we also compare with KrylovKit.jl, a state-of th-the-artt package containing efficient and stable implementations of various Krylov subspace methods [33].
In the setting of this experiment, we conclude that the performance of PG and ZeroFPR varies depending on the formulation that is considered, thus motivating the theoretical framework presented in Section 3. In particular, state-of-the-an thrrttrrttttt performance is achieved when using the newly proposed formulations (n)-(o). Moreover, we observe that depending on the shape of the data matrix, either the primal or the dual formulation is preferred. This is to be expected since the number of decision variables varies accordingly. We also remark that for classical variance-barc -ased formulations, these first-order methods do not perform well.
We also investigate the effect of the number of principal components and the sizes of the matrices involved for the best-performing methods in Figure 3. It can be seen that the first-order methods consistently outperform the classical methods for varying number of principal components, provided that the dimensions of the involved matrices are sufficiently large.
Additional experiments showcasing the performance of these first-order methods various formula- tions as well as toy experiments for our robust kernel PCA can be found in Appendix F.
6 Limitations
We showed that the dual is kernelizable if one of the functions in the primal is unitarily invariant. While this may appear restrictive, many interesting functions such the magnitude of the determinant, Ky Fan norms and Schatten p-norms all fall under this category, thereby leading to a whole new family of kemelizable formulations, derived in principled way. Another issue that many kernel-based methods have is scalability, since the size of the kernel matrix scales as 2 Nevertheless, many mitigations have been proposed to deal with this problem, such Nyström approximations [80, 30] or random Fourier features [60], which are also applicable to our framework.
Another limitation is that (DCA) does not yield the eigenvectors but only (orthonormal) basis for the eigenspace. However, this is not really an issue in practice since only a small number of principal components are typically required and the eigenbasis can be recovered in a postprocessing step.
7 Conclusion and future outlook
In this paper we revisited PCA under the light of difference-of-connvx duality. First, we proposed several novel DC pairs that yield more insight into the inner workings of PCA. Moreover, we showed that when one of the terms in the primal is unitarily invariant, the corresponding dual is not only kernelizable but also supports out-of-sample extensions which is essential for modern machine learning applications. Further, showed that applying DCA the variance maximization is related to simultaneous iteration, revealing a novel connection between optimization and numerical linear algebra. In addition, derived new kernelizable DC dual for robust kernel principal component analysis and an associated algorithm, which in turn is connected to iteratively reweighted least squares. Lastly, our experimental results showed that for the correct formulation, simple first-order optimization methods can outperform state-of-the-aa thr- solvers, depending on the required accuracy.
Several interesting research directions for future work remain open. One promising avenue is the exploration of other unitarily invariant objectives in the primal by incorporating domain knowledge or structural priors of the problem at hand. An alternative research direction is to consider deep variants by stacking (kernel) PCA layers [75] where each layer can be described with different formulations. In addition, further line of research could involve developing specialized algorithms for our new DC dual pairs, as currently most of the algorithms are derived for (I) in Figure 2.
Acknowledgments and Disclosure of Funding
This work was supported by the Research Foundation Flanders (FWO) PhD grant 11A8T26N and research projects GO81 I 122NN, G033822N, GOA0920N; Research Council KUL grant C14/24/103, iBOF//3/064; Flemish Government AI Research Program.
The authors thank Alexander Bodard for assistance with the experiments.
10

---

References
[1] P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2008
[2] Foivos Alimisis, Yousef Saad, and Bart Vandereycken. Gradient- subspace iteration methods for the symmetric eigenvalue problem. SIAM Journal on Matrix Analysis and Applications, 45(4):2360-238) 2386, 2024.
[3] Le Thi Hoai An, M Tayeb Belghiti, and Pham Dinh Tao. A new efficient algorithm based DC programming and DCA for clustering. Journal of Global Optimization, 37(4):593-608 2007.
[4] Paolo Antonelli, HE Revercomb, LA Sromovsky, WL Smith, RO Knuteson, DC Tobin, RK Gar- cia, HB Howell, H-L Huang, and FA Best. A principal component noise filter for high spectral resolution infrared measurements. Journal of Geophysical Research: Amoosphrrs, 109(D23), 2004.
[5] Giles Auchmuty. Dual variational priinppppess for eigenvalue problems. In Symposia in Pure Mathematics, volume 45, pages 55-71, 1986.
[6] Alexei Baevski and Michael Auli. Adaptive input representations for neural language modeling. arXiv preprint arXiv: 10853, 2018.
[7] Heinz H Bauschke and Patrick L Combettes. Convex Analysis and Monotone Operator Theory in Hilbert Spaces. Springer, 2017.
[8] Amir Beck. First-order methods in optimization. SIAM, 2017.
[9] Amir Beck and Mare Teboulle. Dual randomized coordinate descent method for solving class of nonconvex problems. SIAM Journal on Optimization, 3177)_1877-1866 2021.
[10] Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning, volume 4. Springer, 2006.
[11] Bernhard E Boser, Isabelle M Guyon, and Vladimir N Vapnik. A training algorithm for optimal margin classifiers. In Fifth Annual ACM Workshop on Computational Learning Theory, pages 144-152, 1992.
[12] Paul S Bradley and Olvi Mangasarian. Feature selection via concave minimization and support vector machines. In International Conference on Machine Learning, volume 98, pages 82-90, 1998.
[13] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. Language models are few-shot learners. Advances in Neural Information Processing Systems, 13:187711901, 2020.
[14] Samuel Burer and Renato DC Montciro. A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Mathematical Programming, 92(2):329-357 2003.
[15] Emmanuel Candès, Xiaodong Li, Yi Ma, and John Wright. Robust principal component analysis? Journal of the ACM, 58(3):1-37, 2011.
[16] Riley Carlson, John Bauer, and Christopher D Manning. A new pair of gloves. arXiv preprint arXiv:2507.18/00 2025.
[17] Jiahao Chen and Jarrett Revels. Robust benchmarking in noisy environments. arXiv preprint 1608 04295, 2016.
[18] Yingyi Chen, Qinghua Tao, Francesco Tonin, and Johan Suykens. Primal-attention Self- attention through asymmetric kernel SVD in primal representation. Advances in Neural Infor- mation Processing Systems, 36:65088 65101, 2023.
[19] Barry A Cipra. The best of the 20th century: Editors name top 10 algorithms. SIAM news, 33(4):1-2, 2000.

---

[20] John P Cunningham and Zoubin Ghahramani. Linear dimensionality reduction: Survey, insights, and generalizations. The Journal of Machine Learning Research, 16(1):2859-2900, 2015.
[21] Arthur Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Stafistical Society: Series B, (9)::1-2 1977.
[22] Li Deng. The MNIST database of handwritten digit images for machine learning research [best of the web]. IEEE Signal Processing Magazine, 29(6):1411142 2012.
[23] Konstantinos Diamantaras and Sun Yuan Kung. Principal component neural networks: theory and applications. John Wiley & Sons, Inc., 1996.
[24] Chris Ding and Xiaofeng He. K-means clustering via principal component analysis. In International Conference on Machine Learning. 2004.
[25] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al. An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint 1XXi:20010.11928 2020.
[26] Carl Eckart and Gale Young. The approximation of one matrix by another of lower rank. Psychometrika, 1(3):211-218 1936.
[27] Alan Edelman, Tomás Arias, and Steven T Smith. The geometry of algorithms with or- thogonality constraints. SIAM Journal on Matrix Analysis and Applications, 20(2):303-353, 1998.
[28] N Benjamin Erichson, Peng Zheng, Krithika Manohar, Steven Brunton, J Nathan Kutz, and Aleksandr Y Aravkin. Sparse principal component analysis via variable projection. SIAM Journal on Applied Mathematics, 80(2):97--0022, 2020.
[29] Jicong Fan and Tommy WS Chow. Exactly robust kemel principal component analysis. IEEE transactions on neural networks and learning systems, 311)7749-761, 2019.
[30] Shai Fine and Katya Scheinberg. Efficient SVM training using low-rank kernel representations. The Journal of Muchine Learning Research, 2:243-264, 2001.
[31] John GF Francis. The QR transformation a unitary analogue to the LR transformation -part 1.
The Computer Journal, 4(3):265-271, 1961.
[32] Zhengyang Geng, Meng-Hao Guo, Hongxu Chen, Xia Li, Ke Wei, and Zhouchen Lin. Is attention better than matrix decomposition? arXiv preprint arXi::2109.04553, 2021.
[33] Jutho Hacgeman. KrylovKit, March 2024.
[34] Philip Hartman. On functions representable as a difference of convex functions. 1959.
[35] Magnus R Hestenes, Eduard Stiefel, et al. Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards, (166):40-4366 1952.
[36] J-B Hiriart-Urruty Generalized differentiabllity/duulttt and optimization for problems dealing with differences of convex functions. In Symposium on Convexity and Duality in Optimization, pages 37-70. Springer, 1985.
[37] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012.
[38] Ian T Jolliffe. Principal component analysis for special types of data. Springer, 2002.
[39] Michel Jourée,, Yurii Nesterov, Peter Richtárik, and Rodolphe Sepulchre. Generalized power method for sparse principal component analysis. The Journal of Machine Learning Research,
1 11(2), 2010.
[40] Cheolmin Kim and Diego Klabjan. simple and fast algorithm for L.-norm kernel PCA. IEEE Transactions on Pattern Analysis and Machine Intelligence, 42(8): 1842):184211855 2019.
12

---

[41] Sehoon Kim, Amir Gholami, Albert Shaw, Nicholas Lee. Karttikeya Mangalam, Jitendra Malik, Michael W Mahoney, and Kurt Keutzer. Squeezeformer: An efficient transformer for automatic speech recognition. Advances in Neural Information Processing Systems, 35:9361-9373, 2022.
[42] Takeshi Kojima, Shixiang Shane Gu, Machel Reid, Yutaka Matsuo, and Yusuke Iwasawa. Large language models are zero-shot reasoners. Advances in Neural Information Processing Systems, 35:22199
2222222212 13, 2022.
[43] Hoai An Le Thi, Van Ngai Huynh, and Tao Pham Dinh. Convergence analysis of difference- of-convex algorithm with subanalytic data. Journal of Optimization Theory and Applications, 179(1):103-126, 2018.
[44] Richard B Lehoueq, Danny Sorensen, and Chao Yang. ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods. SIAM, 1998.
[45] Gilad Lerman and Tyler Maunu. An overview of robust subspace recovery. Proceedings of the IEEE, 106(8))
2018.
[46] Adrian S Lewis. The convex analysis of unitarily invariant matrix functions. Journal of Convex Analysis, 2(1 2(1):173-183, 1995.
[47] Xingguo Li, Junwei Lu, Raman Arora, Jarvis Haupt, Han Liu, Zhaoran Wang, and Tuo Zhao Symmetry, saddle points, and global optimization landscape of nonconvex matrix factorization. IEEE Transactions on Information Theory, (100:3389-33144 2019.
[48] Jie Lin, Ting-Zhu Huang, Xi-Le Zhao, Teng- Yu Ji, and Qibin Zhao. Tensor robust kernel PCA for multidimensional data. IEEE Transactions Neural Networks and Learning Systems, 36(2):66222677,, 2024.
[49] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo. Swin transformer: Hierarchical vision transformer using shifted windows. In International Conference on Computer Vision, pages 10012-10022. IEEE, 2021.
[50] Pingchuan Ma, Stavros Petridis, and Maja Pantic. End-to-end audio-visual speech recognition with conformers. In International Conference on Acoustics, Speech and Signal Processing, pages 7613-7617. IEEE, 2021.
[51] Sebastian Mika, Bernhard Schölkopf, Alex Smola, Klaus-Rohert Müller, Matthias Scholz, and Gunnar Rätsch. Kernel PCA and de-noising feature spaces. Advances Neural Information Processing Systems, 11, 1998.
[52] Nikhil Naikal, Allen Y Yang. and Shankar Sastry. Informative feature selection for object recognition via sparse PCA. In International Conference on Computer Vision, pages 118-825. IEEE, 2011.
[53] Minh Nguyen and Fernando Torre. Robust kernel principal component analysis. Advances in Neural Information Processing Systems, 2008.
[54] Feiping Nie, Jianjun Yuan, and Heng Huang. Optimal mean robust principal component analysis. In International Conference on Machine Learning, pages 1062-1070. PMLR, 2014.
[55] Konstantinos Oikonomidis, Emanuel Laude, and Panagiotis Patrinos. Corward-backwer splitting under the light of generalized convexity. arXiv preprint arXiv:2503. 180988 2025.
[56] Wei Pan, Xiaotong Shen, and Binghui Liu. Cluster analysis: unsupervised learning via su- pervised learning with a non-convex penalty. The Journal of Machine Learning Research, 14(1):1865-1889 2013.
[57] Jeffrey Pennington, Richard Socher, and Christopher D Manning. Glove: Global vectors for word representation. In Empirical Methods in Natural Language Processing, pages 1532-1543, 2014.
[58] Flad Plaut. From principal subspaces to principal components with linear autoencoders. arXiv preprint 1rXiv:1804.10253, 2018
13

---

[59] Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, et al. Learning transferable visual models from natural language supervision. In International Conference on Machine Learning. pages 1748-8763 PmL.R, 2021.
[60] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in Neural Information Processing Systems, 20, 2007.
[61] R Tyrrell Rockafellar. Convex analysis, volume 28. Princeton university press, 1997.
[62] R Tyrrell Rockafellar and Roger J-B Wets. Variational analysis, volume 317. Springer Science Business Media, 2009.
[63] Teodor Rotaru, Panagiotis Patrinos, and François Glineur. Tight analysis of difference-of fiference-off-coverr oo-vce algorithm (DCA) improves convergence rates for proximal gradient descent. arXiv preprint arXiv:2503.080 04486, 2025.
[64] Ahmed Sameh and Zhanye Tong. The trace minimization method for the symmetric generalized eigenvalue problem. Journal of Computational and Applied Mathematics, 123(1-2):155-177 2000.
[65] Lorenzo Stella. ProximalAlgorithmms]] Proximal algorithms for nonsmooth optimization in Julia.
[66] Ying Sun, Prabhu Babu, and Daniel P Palomar. processing, communications, and machine learning. IEEE Transactions Signal Processing, 65(3):794-816, 2016.
algorithms in signal
[67] Johan Suykens. SVD revisited: A new variational principle, compatible feature maps and nonlinear extensions. Applied and Computational Harmonic Analysis, 40(3):600-66 609, 2016.
[68] Pham Dinh Tao and LT Hoai An. Convex analysis approach to DC programming: theory, algorithms and applications. Acta Mathematica Vietnamica, 22(1):289-355, 1997.
[69] Rachel SY Teo and Tan Nguyen. Unveiling the hidden structure of self-attention via kernel principal component analysis. Advances in Neural Information Processing Systems, 37:101393- 101427, 2024.
[70] Andreas Themelis, Lorenzo Stella, and Panagiotis Patrinos. Forward-backward envelope for the sum of two nonconvex functions: Further properties and nonmonotone linesearch algorithms. SIAM Journal on Optimization, 28(3):2274-233,, 2018.
[71] Mamadou Thiao, Pham D Tao, and Le An. A DC programming approach for sparse eigenvalue problem. In International Conference Machine Learning, pages 1063-1070, 2010.
[72] Michael Tipping. Sparse kernel principal component analysis. Advances in Neural Information Processing Systems, 13, 2000.
[73] John Toland. A duality principle for non-convex optimisation and the calculus of variations. Archive for Rational Mechanics and Analysis, 71 71:41-61, 1979.
[74] Francesco Tonin, Alex Lambert, Panagiotis Patrinos, and Johan Suykens. Extending kernel PCA through dualization: sparsity, robustness and fast algorithms. In International Conference on Machine Learning, pages 34379-344 333993. PMLR, 2023.
[75] Francesco Tonin, Qinghua Tao, Panagiotis Patrinos, Johan Suykens. Deep kernel principal component analysis for multi-level feature learning. Neural Networks, 170:578-595, 2024.
[76] Hugo Touvron, Matthicu Cord, Matthijs Douze, Francisco Massa, Alexandre Sablayrolles, and Hervé Jégou. Training data-efficient image transformers & distillation through attention. In International Conference on Machine Learning, pages 10347-10357. PMLR, 2021.
[77] Lloyd N. Trefethen and David Bau. Numerical Linear Algebra. SIAM, 1997.
14

---

[78]
Matthew A Turk, Alex Pentland, et al. Face recognition using eigenfaces. In Conference on Computer Vision and Pattern Recognition, volume 91, pages 586-591, 1991.
[79] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in Neural Information Processing Systems, 30, 2017.
[80] Christopher Williams and Matthias Seeger. Using the Nyström method to speed up kernel machines. Advances in Neural Information Processing Systems, 13, 2000.
[81] Alan L Yuille and Anand Rangarajan. The concave-convex procedure (CCCP). Advances in Neural Information Processing Systems, 14, 2001.
[82] Jun Zhang, Yong Yan, and Martin Lades. Face recognition: eigenface, elastic matching, and neural nets. Proceedings of the IEEE, 85(9):1423-144 1435, 1997.
15

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract contains the claims made in the paper and the contributions highlighted in separate section the introduction.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
The paper has a separate section discussing the limitations.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model sell-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't It acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
16

---

Justification:
The paper states all the made assumptions and the proofs are complete and correct.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: The paper discloses all the information to reproduce the main experimental results.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
17

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
The code for the experiments is provided as supplementary material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https:
wtips.cc/ppblicc//eecccc gwwcc//gueee..ee CoodeSummii agionPoiicy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The paper discloses the complete experimental setting.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: The paper provides the standard deviation of benchmark timings.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
18

---

The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
OK to report 1-sigma error bars, but one should state The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative eтror rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The paper provides sufficient information on the computer resources.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g.. preliminary or failed experiments that NeurIPS Code of Ethics https://neurip:: http://///euupps..cppcppe/////ea httpp//e/p/eeeeeeieceeeeeeeeeell
didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the
Answer: [Yes]
Justification:
The research conducted in the paper conforms to the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This paper presents work whose goal is advance the field Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
19

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The assets used in the paper propely credited.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com com/dataseta has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
20

---

If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmoel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question:
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
21

---

Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that if the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs any important, original, or non-standard components.
Guidelines:
The answer NA
that the method development in this research does involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https http://ies...//eeeeeeeeee httpe//ris..e/eeoeeeeeeeeeeeeeee c for what should or should be described.
22

---

A Facts from convex analysis
We include summary of relevant facts from convex analysis in this Appendix for completeness. Fact A.1 (Equality case of Fenchel- inequality, [7, Theorem 16.29).. Let H be Hilbert space with inner product (s HH. Let f H. + R be convex, closed and proper function. Let E H. then the following are equivalent:
0f(z) + f(x) + f')) (2,u)N = IE 1f0f(())
Fact A.2 (Conjugate separable sum, 18, Theorem 4.12)). Let f R R X ... RRN
E. be given by (2...2.......... IN) Σ f((z;) where each fs R - R is proper. Then
yN)) = Ei= (y:). (1,,,,2,..... 02,...,
Fact A.3 (Conjugates). The conjugate pairs in Table hold. Recall that the indicator function dc of set C. is given by
5o
x EEC, otherwise.
do(x)
1+00
Table 2: Conjugate functions [8, Appendix B]. || ||, denotes the dual norm of || the 1n-norm.
and ||
-lpp denotes
dom(f)
f
Assumption
2 ||
R"
2 Illlll?
0
(z)
B1ll
|lvll
(2)
1 ||
R"
1 ||30
p 1>1, + 11
(3)
p
q
2 ||0||
Ba-1H-I1
a 00
(4)
Fact A.4 (Equivalence unitarily invariant and absolutely symmetric, [46, Prop. 2.2]). function F: gmxn R is unitarily invariant if and only if it can be written as f a where f IR RR is absolutely symmetric and a mR + R29 denotes the singular value map with =min(m, n). iee. a(X (X) )
the vector with components ai (X) Z xx(X) ... (X) 0, the singular values of X.
Fact A.5 (Spectral conjugate, [46, Theorem 2.4]). Let g RR9 - R be absolutely symmetric. Then
(g og) = 9" oa
where a is defined as in Fact A.4.
Corollary A.6 (Conjugates spectral functions). The conjugate pairs in Table 3 hold.
Table Conjugates of spectral functions.
f
dom(f)
f
Assumption
8 S
(X)
||
(5)
||
XX|ss
mmxx
8
(Y)
p> 1, +
1-1
(6)
1p
||
m
||
(7)
5,0
mmx
- 111111,,
p >>
(8)
Fact A.7 (Subdifferential composition with norm, [7, Example 16.73]). Let h R + R be convex and even, and let ff = o || Then, @h( = I-p: for some E R4 and
| r a E an([1=
if + 0,
0f(x)
B
B
if = 0.
23

---

B Missing proofs
B.1 Proof of Proposition 2.3
We recall the proof from [73, Sec. 3.1] for completeness.
Proof. Let Hi Rdx= H₂ R- and denote (1-))11 (- for the corresponding inner products. Since F is convex, closed and proper, we have that F = F [7, Cor. 13.38). Combining this fact with the definition of the convex conjugate yields
inf G(W) F(XW)= inf G(W) F* (XW( WER,
WEH,
inf G(W) - sup {H, ( -F)(H)
WEH1
HEHz
inf
inf G(W) + FH) - XH,
WEH: HeH2
iinf (* sup {(X) H.W',, - G(W))
HEH2
Well
inf HEN,
(H) - G -x
which shows the strong duality of the two problems. Moreover, let W be solution of (P), i...
G(W*)-FFT F(XW*)
min G(W) F(XW) WEH:
F' G (X HER2
inf
Thus we have
G(W") - F(XW") S G H for all H E H2. Let H' e @FF(X"*), then by Fact A.1
F(XW*) + F*(("))) (XW*, H (W*, X H* Hi
which we substitute in (9) with H = 11* to obtain
G(WT) (W*, x 11" H, < G-6 (x' 11") - G(W*) + (X (1TT(XI1))))))) S (W*, mm...
The Senchel-Young inequality [7, Prop. 13.15] yields
G(W*) + (X H') 2 (W*,X) x H² "T,,)
such that we must have equality throughout and therefore
G(W*) -
which shows that H' solves (D) since strong duality holds. The proof that any W* E (x H*)
is a solution of (P) whenever H' is solution of (D) is completely analogous.
B.2 Proof of Proposition 2.5
Proof. The iterates for (DCA) applied (P) and (D) are respectively
H(k) E (F(XW(()); (7(h) E Ox (()));
((++1) e 0G* (X H()),
A((+1) E (XTV()))) (F
where the latter equality follows from the fact that F is convex, closed and proper and therefore F = F [7, Cor. 13.38]. We now proceed by induction on k. The base case follows immediately by our choice of (0) (F(Xw())) Suppose that the statement holds for A: 2 0, we show that it also holds for l: + 1. Indeed, since (((+1) E (x (x)) by the induction hypothesis, we may take (K) (k+1) and
(k+1) H (k+1) E @F((XW(tt))
which implies
((++2) (x ((++1)) 8G" (XTT(XIT
as desired.
24

---

B.3 Proof of Theorem 3.1
Proof. The problem (1) can be written the form (P) via = bB 15
and F 211- The DC
dual follows from Proposition 2.3, (5) and (8).
The problem (n) is equivalent to (l) since multiplying by two and taking the square root does not change the set of minimizers. The problem (n) can be written in the form (P) via G = 1
and
F = ||s.- The DC dual follows from Proposition 2.3, (5) and (6).
The problem (0) is equivalent to (q) since squaring the norm and dividing by two does not change the set of minimizers. The problem (p) can be written in the form (P) via G = 211 -11 and F = Iss The DC dual follows from Proposition 2.3, (7) and (6).
B.4 Proof of Proposition 3.2
Proof. Since the solution of (I) lies at the boundary by the discussion after Proposition 3.2, we may consider the equivalent problem
minimize F(W) 1| WER
tr(W x XW).
(10)
W W==
Writing eigenvalue decomposition X X W (A))(A)W with W real orthogonal, and assum- ing that the (necessarity nonnegative real) eigenvalues are in descending order, find that Diag (A) contains the entries 1 (X)² 2 2 (x) 2 0 on the diagonal by [37, Theorem 2.6.3(b)). Letting
W be the matrix consisting of the first columns of W, and using the orthonormality of the columns of W, we obtain that
F(W*) = (oig>====W"""iii((ooo (X)) (x).,)((()).....ox)))))) W - Σ (x)²
i=1
so the minimal value of (10) is less than or equal to F(W*), On the other hand, from [37, Theorem 4.3.45], we find that F(W)= 2 - - (X)2 such that W* is solution of (10).
B.5 Proof of Proposition 3.3
Proof. The problem (k) can be written in the form (P) via = + 113, and = -111-113, The DC dual follows from Proposition 2.3 and (8).
B.6 Proof of Proposition 3.5
Proof. Since G is absolutely invariant, can be written as g o a by Fact A.4. Then, by using Proposition 2.3 and Fact A.5, we find that (D) becomes
minimize
F* (H) - (g" o oo)(X H). HERNx*
But the singular value vector oXX H) is equal to the elementwise square root of the nonincreasing eigenvalue vector A(H KH), ie.. oXX H) VX(H KH). Therefore, the DC dual (D) takes the form
minimize HERNx*
F (HH 9" (VACH KH)
(A(HTKH)
where A(X) denotes the vector with the eigenvalues of X (necessarily square) in any order and v- is taken elementwise.
B.7 Proof of Theorem 3.8
Proof. Let full SVD of T be given by DDiag(a H*)VVT where U E Rdxd V Raxs are both orthogonal and a(X 11** pmiu(d,s) RR since we assumed rank(X) s
min(N, d) & d. Let ú be after dropping the last d-a s columns, we can then equivalently write
X H UDiaglo(X* HT)
(11)
25

---

where Diaglo(X H*) E R is now square diagonal matrix with oxX H*) on the diagonal.
Since H a solution of (D) by assumption, we know that there cxists a W E 0G (X H*) which is solution of (P) by Proposition 2.3 and Assumption 2.4. By Proposition 3.10, this solution has the form (Diag(u) where Og" (o((" H))) where we again omitted the columns of U that are multiplied with zero.
Therefore, is required to calculate
W (1)a &.
(12)
This is however not yet kerneiizable. To this end, we first note that (11) is equivalent to
Ü == XT (Diaggg(a H²)))-
where we used the assumption that = KHH nonsingular so that the singular values of X I1* are nonzero. Secondly, we have an eigendecomposition
H KH = VDiag(A)VV V (oDiag(o (
such that Diag I)))) Diag(A) where a shorthand for A(H1 (KH* KH*).
Using these two insights and (12), we find the desired result
W+ x VDiag(u @ A 1/3VVT H+ X
where j E 0g* (VX) and both vÀ and A-1/2 1/2
elementwise.
B.8 Proof of Theorem 3.11
Proof. Without loss of generality, we consider the statement for X X. (DCA) for (d) and simultane- ous iteration are shown in Algorithms 7 and respectively, Here, QR denotes compact QR decom- position. The result then follows from noting that QR and SVD compute orthonormal bases for the column spaces of their arguments in U and respectively. Moreover, if (R)) and Q(A) span the same space and v(k) Rxxr has full rank, then we also have col( xUUk)) - col( col(XXX((A) ( where col( A) denotes the column space of A. Therefore, (T(k) and Q(k) form orthonormal bases for the same subspaces every iteration and are related through orthogonal transformation by [37, Theorem 2.1.18].
Algorithm 7 (DCA) for (d)
Algorithm Simultaneous iteration for X X
Require: Matrix X E R
Require: Matrix E Nxd
1: Initialize (')) E Rdxx
Initialize W7() IRdx
2: Initialize & 0
2: Initialize x- 0
3: repeat
repeat
4: (U(6), (k), + SVD(X x
4: [Q(K), R((K) + QR(X XW(K) QR(X)XW(K)
5: ((++1) + ((k)
5: W(K+1) +- +Q(h))
6: k+*+t k+1
k- k +
7: until convergence
7: until convergence
B.9 Proof of Proposition 4.1
Proof. First, the objective from (c) without the square can be written as
minimize WER
(X.. WWT TXX WWT X
WTW W=I,
which is equivalent to
N Σ
minimize 4x
13 - 1
W W=l,
1=1
26

---

where we expanded and used the constraint. This minimization problem can be written as
minimize W
f(((X))..))
W=
where
x 12 2 1
| S |X,,:l
+00
otherwise,
which is convex function (see [61, p. 106]). The separable sum of convex functions is again convex and the objective therefore maximize convex functiin. By the same reasoning in Section 3.2, the constraint set can again be relaxed to its convex hull, being the spectral norm unit ball. The result now follows from taking G = 8.1-1 and F = Ei= fi((-),)), where the conjugates follow by (5), and combining (4) with Fact A.2.
C DC algorithms
This Appendix is devoted to deriving the many algorithms obtained from (DCA). To this end, first recall some known subdifferential results.
Proposition C.I. Let = || Ilsp then
ÚFT E OF(() = {U V 1pe (i = SVD(X))
where || ||| denotes the I--orrm and Ü a compact SVD of X.
Proof. The expression of the subdifferential follows immediately from Proposition 3.10 and the definition of the Schatten 1-norm. To show that ÜVT is subgradient X, from [61, Theorem 23.8] we have that
(X)o( (X)( x (X)) x Ox)(())xmm (X)mnn())))))) (,)))
where we assume that X Nxd Moreover, recall that
[1}
if 0,
0|0
1,1
if 0,
[-1}
if < 0.
Therefore μ sign(a(X)) E 0l]a(X)]l where sign(x: is the elementwise sign function with the convention (0g For this choice of μ. we obtain
ūVT = U Diag(u)VT E OF())
since the singular values are always nonnegative and the columns/rows of U, V corresponding to zero singular values multiplied with zero the sum.
Remark C.2. It should be noted that when one of the singular values zero, possible to choose any value in the interval [-1, 1] for the corresponding entry in μ. While this is generally not encountered in practice, it could lead to some interesting algorithms to deal with degenerate cases.
Proposition C.3. Let F || I|3_ then OF(x) (X).
Proof. The statement follows from the fact that F is differentiable and [8, Theorem 3.33].
Proposition C.4. Let F = || lls,,, then |X
X E ÜF(X) for X + 0.
Proof. The result follows immediately from Proposition 3.10 and the fact that 0|| { { } for * since is differentiable in that case.
Remark C.5. The case X the preceding Proposition can be handled by using the fact that 01[0||2 B 1-1
27

---

Proposition C.6. Let F= 2111 -1s, then tr(E)ÜVT = | UVT OF((X), where USVT T is compact SVD of X.
Proof. Note that F g o H where g(t) 2 max(0, t)2 is nondecreasing, differentiable convex function and H = || lls, convex. By the chain rule subdifferential calculus [8, Theorem 3.47], have that (F (((())))H(X)) The desired subgradient then follows after using Proposition C.I.
Proposition C.7. Let = |ll - S.ys then OF(X) = {ŪΣV/VTT, where ÜEVT is compact SVD X.
Proof. Note that F 21l | where || |l4/3 denotes the 1//3-morm
14//3 R R: (|21 4//3 + |x2|1/3
|z 1/333//4
Moreover, we have that the gradient 11 Ilb,/s given by the elementwise eube root its argument.
The desired result then follows from Proposition 3.10 and the same arguments from the proof of Proposition C.1 use compact SVD instead of full SVD.
C.1 Derivation of Algorithms and 2
Proof. The formulation (I) (which is equivalent (d) after relaxing the constraint), can be written
as (P) with G = 8 and F = 2 -13.1 Algorithm then follows from (DCA) by noting that 0F(X) = } (Proposition C.3) and UVT E (Y) where ÜΣVT SVD(Y) by (5) and Proposition C.1. The derivation Algorithm follows analogously.
C.2 Derivation of Algorithm 3
Proof. The primal formulation of Proposition 3.3 can be written as (P) with G = 1l l1... and F = 2111 113,- Algorithm 3 then follows from (DCA) by using (F = {X} (Proposition C.3) and 0G*(Y) = (US1/STT where UΣVT SVD(Y) by (8) and Proposition C.7.
C.3 Derivation of Algorithm 4
Proof. The dual formulation of Proposition 3.3 can be written as (P) with F = 1 Slyla G= 211-1,, where T takes the role of X. Using (DCA), OG* (X) = {X) (Proposition C.3) and Proposition C.7 yields the updates
4/3
2,(,), (1), v(A) SVD(X H(H),
((++1) XU11)210)///1111
Now note that
((() = X ((0)))
where (k) is invertible since compact SVD is taken. Moreover, a compact eigendecomposition of (()) (*()), yields (*4 (k)) (k)) (H()) where ( {))) The update is therefore equivalent to
[v(k) A + (R) KH(#),
H(+t1) A4) VAA)
which is exactly Algorithm 4.
C.4 Derivation of other DC algorithms
In this subsection, we consider kernelizable dual updates of (DCA) applied to the problems (!) through (q). First, we have the following general result when G' is unitarily invariant.
28

---

Proposition C.8. Let X E RNxd be the data matrix, s S rank(X) and define the kernel matrix K XXT. Consider the primal problem (P) and suppose G = g a is unitarily invariant. Assume
H()) ( invertible for each k. then the iterations of (DCA) applied the dual are [V(M) Diag(A(A)) + (HG( (H)
H(++1) OF (KH() Diag(u Ø v V((i))(())
where ((()) E Og" (v A()). We use the notation @ and v to denote elementwise division and square root respectively.
Proof. The dual problem (D) is
minimize H
(H) - (g'
- oo)(
from Proposition 3.5. By using Proposition 3.10 (note that it contains a full SVD and not a compact SVD), the iterations of (DCA) applied to the dual problem are
(U Diag(o!!)), v(A) T ] SVD(X H()), ((++1) +- OF (1)
where (1) E (a(k)) Since H T (H)) nonsingular, each entry of (kk) is strietly greater than 0 and may use compact SVD and compact eigenvalue decomposition. The result then follows from the same arguments as in Subsection C.3 and the fact that F = F [7, Cor. 13.38]. O
Remark C.9. Instead of the assumption that (k) KH (k) is invertible for each k, the result also holds whenever y/s 0 implies that 0 0g"(2) for each component and 0/0 is taken to be O. This the case for all the functions encountered previously (Propositions C.I through C.7).
Applying this result to the formulations (I) through (q) (or transposed variations) yields the updates in Table 4. It is clear that many of these algorithms are tightly related. For example, (14) is equivalent to Algorithm 2 as well as (16) since the scaling only affects (A) which is discarded. Similarly, (18) the same as (14) up scaling. Moreover, (15) is (13) with additional normalization while (15) is also equivalent to (17) since the multiplication with the trace in (17) is void due to the normalization afterwards. some sense, (DCA) blind to simple problem transformations. Nevertheless, while they may be mathematically tightly related, the practical performance may differ significantly.
C,5 Derivation of Algorithms and 6
Proof. Recall from the proof of Proposition 4.1 that the primal problem can be written (P) with G GB and F == Σm fi(()).). The expression for öG* follows immediately from (5) and Proposition C. while for OF we use [61, Theorem 23.8] 3] to find that the subdifferential a separable sum is the Cartesian product of the subdifferentials. Moreover, note that each fs can be written as h4 o || where
I2 a if aa S ||x otherwise
h
R + I
to0
such that we can use Fact A.7. More concretely, we have that
(a)
a
for [a] < ||X,,.ll and the empty set otherwise. Combining these subdifferentials with (DCA) yields Algorithm 5. Note that due to the initialization as well as the update rule for W the spectral norm of (KF)) is always less than such that the argument of dh, has absolute value |a] < ||xX,.I| and the square root is always well-defined. To derive Algorithm 6, the same subdifferentials can be used in combination with the same arguments as the derivation of Algorithm 4 to make it kernelizable. O
Remark C.10. Whenever (XW())), 0 in Algorithm 5, Y can be set arbitrary vector in some ball of some radius, in accordance with the second case in Fact A.7. Another edge case occurs 2| = IXXX (XX)(()) ).,. practice, this is alleviated by adding some regularization constant е in the denominator. Formally, this is accomplished by modifying the first case in the
whenever || ||
definition of h; to be X a² + 2 for la| VIJ 2 + 2 as in [9].
29

---

Table 4: DC algorithms for formulations (1) through (q) where the data matrix is X E INd Only the kernelizable versions involving K = X are provided.
Problem
(DCA) dual updates (kernelizable versions)
Derived through
[V(M] A (IG(H(() KH())) H(*+1) K (
Proposition C.I
(l)
(13)
Proposition C.3
(U(H)) EX4,, V(H) 1 SVD(KH()) H v(K)
Proposition C.3
(m)
(14)
Proposition C.I
[V(], ((0) + (H(() KH(K))
Proposition C.I
(n)
l((+1) ( (() 1/2 V(K)
(15)
Proposition C.4
H(++1) H 82
[1 (k), v(k) (VD((T(H (H 1/20))1112K11 (12
Proposition C.4
(0)
(16)
H v())
Proposition C.1
[[(++1) (V(M),(0)) (IG( KH())
Proposition C.6
(p)
1/2 V(A)
(17)
Proposition C.4
H(++1) MR+b S2
go+u
(0(0) E(4), 1+ SV H(++1) tr(11()) 11/01)-11/ ((
Proposition C.4
(q)
(18)
Proposition C.6
D IRLS and DCA
In this Appendix, derive iteratively reweighted least squares (IRLS) algorithm for the robust PCA formulation from Section 4 and compare with Algorithm 5. As common in the literature (see [45] and references therein), the IRLS update for the robust subspace recovery problem is of the form
W(++1) E arg min WeRdxa [ 21
Σ
(
IIx WWT x
where
((((())
||x,, (w() Xxall-
This minimization problem is equivalent to
((++1) E arg min ||BX WW BX arg min WeRdx. [W 51
1 BXW
de [W] $1
where a square diagonal matrix with (3,(W()))) (1//2) on the diagonal. The latter is recognized to be (1) such that may Proposition 3.2 obtain a solution of the inner problem. This leads to Algorithm 10, where we also used the fact that = UEVT a SVD of A, then A VΣ T
SVD of AT. Note that for each 1, (xk) (k) I. such that g((h) could equivalently be written || 1| (()xppp)))ll which exactly the denominator on line 5 of Algorithm
The main difference between the two algorithms is that an additional square root is taken of a(() in Algorithm 10, while in Algorithm 9 multiplication with (XW())) )..: is performed instead. Another difference is that an orthogonal factor remains Algorithm which is reminiscent of the connection between Algorithm 7 and simultancous iteration, see Theorem 3.
30

---

Algorithm 9 (DCA) for Proposition 4.1 (primal)
Algorithm 10 IRLS for robust subspace recovery
Require: Matrix X Wxd
Require: Matrix X E R
1: Initialize W E Ixxx s.t. || =
Initialize W(0) E RRxx s.t. || =
2: Initialize k
2: Initialize h <- 0
3: repeat
3: repeat
4: for i
N do
4: for i to N do
(XW())
5: B.(a) |X W(k Xx.ll
5: Y.: + VIIx.,112 II(
6: end for
B(b) Diag( -1/2)
6: end for
7:
( I(K), V( J -=S--==V'" SVD(X
8: [U(K) (1)), ((k)
II ] SVD(X B(h)
8: W(+H)) - v())
9: W(++1) U.(k)
10: ++ :,1:
9:
k- --+1
10:
until convergence
11: until convergence
E DC formulations for sparse PCA
The DC duality framework is very broad. We now discuss two promising avenues for sparse PCA.
First, similar to formulation (q), can consider the DC formulation
1 W - ||X
minimize WER 2
where || lls is sparsity inducing such taking the sum of the absolute value of its elements (l-norm on its vectorization). Note that this not a classical sparse PCA formulation but achieves a similar purpose. The convex conjugate of this function follows from (1).
Second, if we consider only one principal component, then a DC formulation for sparse PCA is described in [9]. This formulation takes the form
minimize (c(w) 1 1 Xull3 wER"
2
where C' is the convex hull of the intersection of the unit norm ball and {w Rd ||wloo S k). The conjugate of this indicator function be computed in closed form. Indeed, (2) max w To by definition. Then, since w has at most & nonzeros, we see that this maximization problem is solved when the support of corresponds to the & largest entries of in absolute value. Taking into account that w is also a unit vector, it follows that (z) is exactly the Euclidean norm of the vector containing the l largest entries of in absolute value.
Both of these formulations are not unitarily invariant and therefore not easily kernelizable through Proposition 3.5 (a sufficient condition for kernelizability) Though, to the best of our knowledge, sparse PCA is generally not used in the kernel setting since the main motivation for sparse PCA is interpretability of the data and principal components. More concretely, by restricting the cardinality of the principal components, this implicitly assumes that they are linear combinations of just a few innput variables, which does not mesh well with kernel methods that inherently nonlinear.
F Additional experiments
In this Appendix, we provide further experimental results.
Higher accuracy Timing results for higher target accuracies than those Section are presented in Table 5. The problem instances are generated in similar manner. We observe that the Krylov-type methods are more suited than the first-order methods this setting. Nevertheless, Zero applied to (n) still outperforms the Krylov methods for larger matrices. We observe that the (DCA) algorithms are quite slow, which is mostly due to the relatively expensive check the stopping criterion each iteration.
31

---

Table 5: Timing results for various methods applied PCA formulations with higher accuracies. The problem setting (N, d, s, E) is the same as in Table 1. All timings are in milliseconds, and timings longer than seconds are not displayed.
Method
(3000, 2000, 20, 10-5)
(2000, 3000, 20, 10 ()
(3500, 3500, 20, 10 10-5)
ZeroFPR ())
2009.0 + 122.5
2202.6 296.0
4269.3 + 57.6
ZeroFPR (n)
735.0 + 295.7
876.1 = 194.8
1269.1 357.0
ZeroFPR (o)
1092.5 + 285.1
1532.8 213.1
2577.7 = 274.3
ZeroFPR (q)
2018.4 + 141.1
2148.0 = 319.6
4245.3 = 715.5
PG (n)
3569.9 + 410.1
PG (o)
4011.7 + 969.8
Algorithm
4709.5 + 3963.0
Algorithm 2
2559.5 + 832.2
Algorithm 3
Algorithm 4
SVDS
625.3 +
10.6
556.9 +
18.7
1564.8 + 31.1 1858.7 + 24.3
KrylovKit
773.8 +
29.7
696.0 17.5
Different spectra The performance of Krylov methods for computing the top singular values is known be highly sensitive the spectrum, and in particular the spectral gaps az Fi41- Table reports timing results for various methods applied to matrices with different fixed spectra. In the first column, where the singular values decay exponentially, all methods perform equally well. In the second column, the singular values decrease linearly from 500 to 100, resulting in small spectral gaps. Here, the Krylov methods are roughly 10 times slower compared the exponential decay case. Notably, ZeroFPR and PG on the formulations (n) and (o) significantly outperform the Krylov solvers. This trend continues in the third column, where the singular values decrease linearly from 500 to 10-8.
Table 6: Timing results for various methods applied PCA formulations for matrices with specified spectrum. All matrices are square with N =3500. A stopping criterion tolerance of 10-3 was used to find s = 20 principal components. All timings in milliseconds, and timings longer than
seconds are not displayed.
Method
a, 100 0.9
0, 500
1)
a; 500 (i 1) 500 10
ZeroFPR ())
490.7 + 123.5
ZeroFPR (n)
171.9 + 64.8
456.2
90.1
455.7 + 80.0
ZeroFPR (o)
260.1 + 68.1
1043.3 + 164.0
1137.7 + 165.6
ZeroFPR (q)
298.5 + 97.2
PG (n)
115.8 + 22.1
730.2 + 55.0
729.5 + 71.1
PG (o)
262.1 40.3
1068.8 + 34.3
1173.2 + 88.5
Algorithm 1
300.2 + 13.9
3633.2 + 184.8
3581.2 + 554.2
Algorithm 2
295.7 + 14.5
3594.1 + 36.0
3866.0 + 59.4
Algorithm 3
429.7 + 19.2
Algorithm 4
387.1 +
16.2
SVDS
288.9 +
2.1
3113.9 + 133.4
2588.5 + 178.1 3653.5 + 309.8
KrylovKit
242.8 +
3.1
4250.3 + 392.9
Real-world datasets Table 7 shows the timing results of the best performing methods the MNIST dataset [22] (N = 60000, d = 784) with a tolerance of E = 10-3 We observe the same behavior in the synthetic experiments: the generic first-order methods are faster than the classical eigensolvers this setting.
Table 8 shows the timing results of the best performing methods on the 100k top words from the 2024 Wikipedia + Gigaword 5, 50d GloVe word embedding dataset [57, 16] (N 100000, 50) with tolerance of 10-3, While the first-order methods are still faster than the classical eigensolvers,
32

---

this difference is less pronounced. We attribute this result to the fact that these word embeddings already form a compressed version of the word vector space, and therefore have favorable spectral characteristics, similar to the results encountered in the first column of Table 6.
Table 7: Timing results for different methods of applying PCA to the MNIST dataset with tolerance e = 10-3 The column headers denote the number of principal components. All timings are in milliseconds, and timings longer than 5 seconds are not displayed.
Method
30
50
100
150
ZeroFPR (n)
286.7 + 92.2
444.7 + 86.1
1734.4 + 247.0
2363.3 + 135.2
PG (n)
215.5 + 35.7
351.9 + 62.4
1254.8 103.1
1772.2
150.5
SVDS
1685.1 + 226.2
2576.5 + 4.5
Table 8: Timing results for different methods of applying PCA to the GloVe word embedding dataset with tolerance of E 10-3, The column headers denote the number of principal components. All timings are in milliseconds.
Method
10
15
20
30
ZeroFPR (n)
131.1
114.5
208.2 + 150.4
271.8
147.2
367.1
198.2
PG (n)
76.0 + 71.4
120.0 + 100.4
165.3 + 111.9
250.4 + 113.2
SVDS
200.6 + 14.6
281.0 + 31.4
319.7 11.5
329.0 34.5
Additional formulations In Table 9, we present some timings for the algorithms described in Table 4, where we used the relative error with respect to the optimal value for the corresponding formulation as a stopping criterion. Wherever possible, we reuse information from the SVDs/EIGs to cheaply evaluate the objective value. We observe that in general, the algorithms using eigenvalue decompositions are faster than those involving singular value decompositions. Moreover, the algo- rithms par with the Krylov solvers whenever the number of features is greater than equal to the number datapoints. This is to expected since the algorithms use the kernel matrix. In the opposite setting, one should use algorithms using the (scaled) covariance matrix instead.
Table Timing results for methods from Table applied PCA formulations. The problem setting (N, s,E) is the same as in Table 1. All timings are in milliseconds.
Method
(3000, 2000, 20, 10
(2000, 3000, 20, 10 10-3)
(3000,3008) 20, 10 10-3)
Iteration (13)
1135.6 + 72.7
520.4 + 58.3
1198.9 +:2 153.3
Iteration (14)
1327.5 + 331.8
567.6 78.6
1179. + 272.4
Iteration (15)
792.4 + 88.1
338.1 + 42.3
800.7 + 79.9
Iteration (16)
1059.4 65.5
556.6 58.7
1182,2 94.0
Iteration (17)
1742.6 : 206.4
774.6 = 84.4
1815.1 + 71.7
Iteration (18)
5926.6 = 149.6
1085.7476.J
4895.1 94.9
SVDS
638.1 + 41.7
522.3 44.4
1164.0 57.5
KrylovKit
755.5 + 26.3
685.4 = 44.3
1215.7 + 43.5
A note on tolerance Previous comparisons used a fixed tolerance across methods, but this is not ideal since (1) the problem formulations differ, and (2) algorithms interpret differently in their stopping criteria.
see this impact, we constructed 10 different 50 x 50 matrices with entries sampled from the standard normal distribution. Then, for each tolerance from 10-1 down to 10-9 on a log-spaced grid, we computed s = 10 components, and evaluated the relative error with respect to the optimal solution on the formulations (n) or (o). We averaged these errors over the 10 different matrices and present the results Figure 4. We observe that KrylovKit.ji reaches machine precision for the function value with e 10 while ZeroFPR and PG perform the worst on (n) and (o).
Robust (kernel) PCA Consider the MNIST dataset [22] with a train-test split of 80-20. To verify the robustness properties of the robust PCA formulation in Section 4, we contaminate 15% of the
33

---

10
10
svns
Krylov
PG (1)
PG (8)
()
ZF (n) ZF (q)
PG (o)
PC (g)
(o)
10
10
10-16
10 10-16
10-1
10
10 5 10-7
10
10-1
10
10-5
10
10
10"
Figure 4: Relative error of the objective value (n) (o) for multiple methods with varying tolerance for the stopping error criterion. f denotes the optimal value while ZF is shorthand for ZeroFPR.
training data with heavy Gaussian noise (ar = 15) and leave the test set untouched. We consider four different settings:
1. Linear PCA on the noncontaminated data (baseline)
2. Linear PCA on the contaminated data
3. Robust PCA the contaminated data using Algorithm 5
4. Robust PCA on the contaminated data using Algorithm with a linear kernel (note that the kernel matrix does not fit in memory but we use its factored representation since the kernel is linear)
For each of these settings, we evaluate the reconstruction error on the noncontaminated test set. These errors are summarized in Table 10. The robust methods outperform standard linear PCA contaminated setting. Further, the similar performance of the two robust methods is to be expected, due to strong duality. Lastly, the top components are less affected by the outliers than the remaining components. This is logical since these latter components explain less 'variance" and it becomes more difficult to distinguish noisc/outliers from data.
Table 10: Reconstruction errors on the test set of various PCA schemes. The column headers denote the number of principal components. The robust PCA formulations outperform standard linear PCA on contaminated dataset.
Method
50
100
150
Linear PCA (noncontaminated)
0.011
0.0057
0.0034
2. Linear PCA (contaminated)
0.062
0.0569
0.0520
Robust PCA (contaminated) with Algorithm 5
0.014
0.0107
0.0090
4. Robust PCA (contaminated) linear kernel with Algorithm 6
0.014
0.0107
0.0090
To further illustrate the strength and usecase of this formulation. in light of Theorem 3.8. To this end, for each of the settings, we train a small multilayer perceptron classifier (1 hidden layer with 20 neurons) on these extracted features. Note that choose very simple classifier so the quality of the features becomes more apparent. Additionally, we also consider the following two settings:
extract robust features
5. Robust PCA on the noncontaminated data using Algorithm
6. Robust kernel PCA the contaminated data using Algorithm 6 with RBF kernel (length scale paramter 0.01) approximated using a Nyström approximation with 500 pivots
The test accuracies of each classifier are displayed in Table 11. always perform better than the non-robust formulations for contaminated data. The robust formulation the noncontaminated data (setting 5) also does not degrade in performance with respect the baseline (setting 1). Lastly, the RBF kernel trained on the contaminated data performs on par with models trained uncontaminated data.
We observe that the robust formulations
34

---

Table 11: Test accuracies of small multilayer perceptron after feature extraction.
Method
50
100
150
1. Linear PCA (noncontaminated)
95.87
95.73
94.84
2. Linear PCA (contaminated)
76.57
79.93
81.09
3. Robust PCA (contaminated) with Algorithm 5
86.89
87.32
85.53
4. Robust PCA (contaminated) linear kernel with Algorithm 6
87.63
84.44
86.92
5. Robust PCA (noncontaminated) with Algorithm
95.73
96.06
95.43
6. Robust PCA (contaminated) RBF kernel with Algorithm 6
95.21
95.43
94.34
As last experiment, we compare the robust PCA formulation with the classical robust PCA [15]. Note that the two approaches of different nature. The classical method decomposes matrix into a low-rank component and a sparse component, which is well-suited for motion segmentation. In contrast, the formulation in Section 4 is based the autoencoder perspective and aims to enforce sparsity on the reconstruction errors, which better aligned with outlier detection. Another major difference is that the classical method is not out-of-sample applicable while ours is through Theorem 3.8. These differences make fair comparisons difficult.
To further illustrate this fact, consider the synthetic experiment from [15, Section 4.1]. We first generate a low-rank component according to the same settings and add sparse matrix where E has limited support, with the nonzero entries of E being independent Berouuli +1 entries. We consider two settings:
1. The support of E is chosen uniformly distributed over all its entries
2. The support of E consists exclusively of full rows (i.e., an image where complete rows are perturbed)
We then compare the classical principal component pursuit (PCP) with alternating directions [15, Algorithm 1] with Algorithm For our experiment, we choose rank(L*) = and L* E 500x000 We assume the size of the support of E' is [E' 12500 and compute both the reconstruction error of the low-rank component and the primal cost of Proposition 4.1 (i.e., the (norm the row-wise reconstruction errors) and summarize different random problems). We observe that the classical PCP performs better in the setting where was proposed while algorithm is better in the outlier setting. Moreover, cost lower for our algorithm both cases, which is to be expected since the PCP algorithm is not designed to minimize this cost.
results in Table 12 (all metrics were averaged over 20
that the primal
Table 12: Comparison principal component pursuit (PCP) and Algorithm for decomposing a matrix into a low-rank component L and a sparse component E.
||-- [
Setting
Method
Primal cost
Uniform perturbation
PCP
3 x 10-8
2488 2364
Algorithm 5
0.013
Full row perturbations
PCP
0.045
559
Algorithm 5
0.009
545
35


Reviews:
[{"review_id": "sOo29vMbmy", "reviewer": "Reviewer_BYYV", "summary": "This paper revisits Principal Component Analysis (PCA) through the lens of difference-of-convex (DC) duality, particularly Toland duality, and presents several contributions:\n1. New DC Dual Pairs for PCA: The authors derive three new dual formulations of PCA and show that classical algorithms (like simultaneous iteration and the QR algorithm) can be interpreted as instances of the Difference-of-Convex Algorithm (DCA).\n2. Robust and Kernelizable PCA: A new DC formulation is proposed for robust PCA, including a dual that is kernelizable and supports out-of-sample extension, which is significant for real-world applications.\n3. Efficient Algorithms: Several DCA-based algorithms are proposed and compared empirically against classical and modern methods, showing strong performance, especially with newly introduced formulations.", "questions": "1. Have you tried your robust kernel PCA formulation on any real-world datasets with known outliers? How does it compare to classic robust PCA variants?\n2. You mention that standard kernel PCA methods suffer from scalability issues. Have you evaluated the performance of your kernelizable DC duals with Nyström approximation or other scalable kernel methods?\n3. You report that formulations (n)-(o) perform best in your experiments. Do you have any theoretical insights into why these perform better than others (e.g., (l)-(m))? Is it due to smoother optimization landscapes or better conditioning?\n4. You discuss the link between PCA and linear autoencoders. Could the DC duality framework also inspire new architectures or training methods for nonlinear or deep autoencoders?", "limitations": "Yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n\nThe reinterpretation of PCA through DC duality, and the discovery of new dual pairs (beyond existing literature), are mathematically insightful. Showing that classical methods like simultaneous iteration are instances of DCA gives new understanding into optimization perspectives in numerical linear algebra. The kernelizable robust PCA formulation is a practical strength that aligns well with modern machine learning use cases.\n\nWeaknesses:\n\nWhile the theoretical results are strong, it's unclear how well the proposed methods extend to large-scale, high-noise, or real-world datasets. The experiments are mostly synthetic or based on standard PCA settings. A demonstration on real-world tasks (e.g., image compression, NLP embeddings) would be helpful.", "was_revised": true, "final_justification": "I appreciate the effort the author put into addressing my concerns, and I’m satisfied with the explanation made by the authors.", "rebuttal": "We thank the reviewer for their positive evaluation of our work and interesting questions.\n\n---\n\nIn the following we address the **Weakness** stated by the reviewer.\n\n> While the theoretical results are strong, it's unclear how well the proposed methods extend to large-scale, high-noise, or real-world datasets. The experiments are mostly synthetic or based on standard PCA settings. A demonstration on real-world tasks (e.g., image compression, NLP embeddings) would be helpful.\n\nWe will add the following experimental results to the manuscript. Consider the following 2 datasets.\n- MNIST dataset (60000x784)\n- The 100k top words from the 2024 Wikipedia + gigaword 5, 50d GloVe word embeddings dataset (100000x50)\n\nWe compare the computational times of classical PCA eigensolvers against our methods with a required tolerance of $\\varepsilon = 10^{-3}$. The results for the best algorithms are summarized in the following table where the column header denotes the number of principal components (no standard deviation is provided for timings larger than 5 seconds).\n\nMNIST:\n|             | 10                       | 30                         | 50                       | 100                        | 150                        |\n| ----------- | ------------------------ | -------------------------- | ------------------------ | -------------------------- | -------------------------- |\n| ZeroFPR (n) | $125.880 \\pm 28.716$ ms | $286.692 \\pm 92.200$ ms   | $444.745 \\pm 86.127$ ms | $1734.389 \\pm 247.000$ ms | $2363.292 \\pm 135.235$ ms |\n| PG (n)      | $111.101 \\pm 9.363$ ms  | $215.460 \\pm 35.660$ ms   | $351.948 \\pm 62.360$ ms | $1254.828 \\pm 103.131$ ms  | $1772.221 \\pm 150.473$ ms |\n| SVDS        | $820.029 \\pm 85.984$ ms | $1685.063 \\pm 226.155$ ms | $2576.462 \\pm 4.456$ ms | $5251.069$ ms                | $7850.230$ ms                |\n\nGloVe:\n|             | 5                        | 10                         | 15                         | 20                         | 30                         |\n| ----------- | ------------------------ | -------------------------- | -------------------------- | -------------------------- | -------------------------- |\n| ZeroFPR (n) | $61.963 \\pm 79.404$ ms  | $131.065 \\pm 114.521$ ms | $208.227 \\pm 150.357$ ms | $271.817 \\pm 147.240$ ms | $367.117 \\pm 198.195$ ms |\n| PG (n)      | $43.312 \\pm 59.676$ ms  | $75.954 \\pm 71.368$ ms   | $119.950 \\pm 100.352$ ms | $165.287 \\pm 111.865$ ms | $250.385 \\pm 113.187$ ms |\n| SVDS        | $136.180 \\pm 19.130$ ms | $200.576 \\pm 14.643$ ms  | $281.008 \\pm 31.438$ ms  | $319.668 \\pm 11.535$ ms  | $328.987 \\pm 34.468$ ms  |\n\nIn both of the experiments, we observe the same main observations as from the synthetic experiments: the generic first-order methods are faster than the classical eigensolvers in this setting. However, we also remark that the difference for the GloVe dataset is quite small. We believe this is the case because the word embeddings have favorable spectral properties such that the results are quite similar to the first column of Table 7.\n\n---\n\nIn the following we answer the posed **Questions**.\n\n> Have you tried your robust kernel PCA formulation on any real-world datasets with known outliers? How does it compare to classic robust PCA variants?\n\nTo demonstrate the robust feature extraction of the robust PCA formulation, consider the MNIST dataset (train-test split 80-20) and contaminate 15% of the training data with heavy Gaussian noise. We then consider the following settings:\n1. linear PCA on the noncontaminated data (baseline)\n2. linear PCA on the contaminated data\n3. robust PCA on the contaminated data using Algorithm 5\n4. robust PCA on the contaminated data using Algorithm 6 with a linear kernel (note that the kernel matrix does not fit in memory but we can use its factored representation $XX^\\top$ since the kernel is linear)\n\nFor each of these settings, we evaluate the reconstruction error on the test set, which is not contaminated! These errors are summarized in the following table (the column headers denote the number of principal components).\n|           | 50    | 100    | 150    |\n| --------- | ----- | ------ | ------ |\n| Setting 1 (noncontaminated) | 0.011 | 0.0057 | 0.0034 |\n| Setting 2 (contaminated)  | 0.062 | 0.0569 | 0.0520 |\n| Setting 3 (contaminated)  | 0.014 | 0.0107 | 0.0090 |\n| Setting 4 (contaminated)  | 0.014 | 0.0107 | 0.0090 |\n\nComparing Settings 3 and 4 with Setting 2, we observe that our robust formulation better reconstructs the test set when trained on the contaminated data. The fact that Setting 3 and 4 have similar performance is expected from strong duality. Moreover, the top components are less affected by the outliers than the remaining components. This is logical since these latter components explain less “variance” and it becomes more difficult to distinguish noise/outliers from data.\n\nTo further illustrate the strength of the formulation, we can now extract robust features in light of Theorem 3.9. We train a small MLP classifier (1 hidden layer with 20 neurons) on these extracted features (we choose a very simple classifier so the quality of the features becomes more apparent). Additionally, we also consider the following two settings:\n\n5. robust PCA on the noncontaminated data using Algorithm 5\n6. robust kernel PCA on the contaminated data using Algorithm 6 with a RBF kernel ($\\gamma=0.01$) approximated using a Nyström approximation with 500 pivots\n\nWe summarize the test accuracies in the following table (the column headers indicate number of principal components at the intermediate stage)\n\n|                             | 50    | 100   | 150   |\n| --------------------------- | ----- | ----- | ----- |\n| Setting 1 (noncontaminated) | 95.87 | 95.73 | 94.84 |\n| Setting 2 (contaminated)    | 76.57 | 79.93 | 81.09 |\n| Setting 3 (contaminated)    | 86.89 | 87.32 | 85.53 |\n| Setting 4 (contaminated)    | 87.63 | 84.44 | 86.92 |\n| Setting 5 (noncontaminated) | 95.73 | 96.06 | 95.43 |\n| Setting 6 (contaminated)    | 95.21 | 95.43 | 94.34 |\n\nWe observe:\n- The robust formulations always perform better than the non-robust formulation for the contaminated data.\n- Using a robust formulation (Setting 5) on the clean data does not degrade performance with respect to the standard formulation (Setting 1).\n- The RBF kernel trained on the corrupted data performs on par with models trained on the clean data.\n\nFor a comparison with a classical robust PCA formulation (that is not out-of-sample applicable) on a synthetic experiment, see our response to **Reviewer FsFQ**. All the code will be made public later.\n\n> You mention that standard kernel PCA methods suffer from scalability issues. Have you evaluated the performance of your kernelizable DC duals with Nyström approximation or other scalable kernel methods?\n\nSee our answer to the previous question (a Nyström approximation is required for the MNIST dataset since storing the complete kernel matrix of the training data would take 20+ GB of memory).\n\n> You report that formulations (n)-(o) perform best in your experiments. Do you have any theoretical insights into why these perform better than others (e.g., (l)-(m))? Is it due to smoother optimization landscapes or better conditioning?\n\nWe indeed believe that the formulations (n)-(o) have a better condition number than the formulations (l)-(m), leading to their better performance. We’d like to point out an analogy with solving linear least squares where it is known that solving the normal equations is impractical, as its condition number is squared, which is similar to formulations (l)-(m).\n\n> You discuss the link between PCA and linear autoencoders. Could the DC duality framework also inspire new architectures or training methods for nonlinear or deep autoencoders?\n\nAs briefly noted in the introduction, the transformer architecture bears a close connection to kernel PCA, which has already motivated the development of a robust transformer variant in [1]. Similarly, an alternative transformer design was proposed in [2], based on another duality framework. Deep autoencoders derived from PCA-based formulations are also identified in our future work, and related architectures based on alternative formulations have been explored as well, see [3]. We therefore indeed believe that the DC duality framework can inspire new architectures as well as efficient training methods.\n\n---\n\nReferences\n\n[1] Rachel SY Teo and Tan Nguyen. Unveiling the hidden structure of self-attention via kernel principal component analysis. Advances in Neural Information Processing Systems, 37:101393–101427, 2024.\n\n[2] Yingyi Chen, et al. Primal-attention: Self-attention through asymmetric kernel svd in primal representation. Advances in Neural Information Processing Systems, 36:65088–65101, 2023.\n\n[3] Tonin, Francesco, et al. Deep Kernel Principal Component Analysis for multi-level feature learning. Neural Networks 170 (2024): 578-595."}, {"review_id": "KgJk3iWILW", "reviewer": "Reviewer_FsFQ", "summary": "This paper revisits principal component analysis (PCA) through the lens of difference-of-convex (DC) duality, proposing several formulations, interpreting classical iterative methods such as simultaneous iteration as instances of DC, and introducing kernelizable formulations for robust variants of PCA.", "questions": "1. Include convergence guarantees or at least local convergence conditions for each of the proposed DCAs.\n2. Benchmark against methods for sparse/robust PCA on real datasets (e.g., MNIST, Yale faces, motion segmentation).", "limitations": "Yes.", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The paper presents a theoretical reinterpretation of PCA via DC duality and identifies connections with classical linear algebra routines.\n2. The use of Toland duality is methodologically sound and allows formulation of new primal-dual pairs.\n3. The kernelization of DC duals via unitarily invariant primal functions is a useful insight.\n\nWeaknesses:\n1. This paper fails to provide any theoretical convergence guarantees for DC algorithms. No rate of convergence, global optimality conditions, or even local convergence analyses are presented. Given the non-convex nature of the problems, such results are crucial for justifying the use of DCA in practice.\n2. The empirical study does not convincingly demonstrate the superiority of the proposed DC methods over classical PCA solvers. In Table 1, DCA-based algorithms (Algorithms 1–6) often take longer than standard eigensolvers, such as SVDS and KrylovKit, with no clear accuracy benefit. The best-performing formulations are still dominated by classical methods in terms of both speed and simplicity, which weakens the practical motivation for these DC-based approaches.\n3. The numerical experiments are limited to synthetic Gaussian data with a fixed dimensionality and target rank. There is no demonstration on real-world datasets, and no evaluation of approximation quality, numerical stability, or robustness to noise.", "was_revised": true, "final_justification": "Thanks, authors for the detailed the explanation. The referee is generally happy and positive about this paper. However, the referee feels it requires another round of review to check the correctness and proof details of the claims made. And, the convergence to the critical point is not satisfactory, since we know PCA can be polynomially solved to optimality. The authors' explanation of sparse PCA is interesting and should be incorporated into the manuscript. Again, referee is positive about this paper and retains the same score for this paper.", "rebuttal": "We thank the reviewer for their positive feedback and useful comments to improve our work.\n\n---\n\nIn the following we address the stated **Weaknesses**.\n\n> This paper fails to provide any theoretical convergence guarantees for DC algorithms. No rate of convergence, global optimality conditions, or even local convergence analyses are presented. Given the non-convex nature of the problems, such results are crucial for justifying the use of DCA in practice.\n\nWe did not provide any theoretical convergence guarantees for DCA due to space restrictions and the fact that it is a well-established algorithm whose theoretical basis is described in full detail in the literature, see for example [1]. Upon acceptance of the paper, we will make the paper more self-contained and further enhance the theoretical aspects by adding the following paragraph in the DCA section 2.3 regarding the basic convergence properties of DCA in general.\n\nConsider the problem (P) and suppose that the iterates $\\{W^{(k)}\\}$ are generated by (DCA), then\n- $G(W^{(k+1)}) - F(XW^{(k+1)}) \\leq G(W^{(k)}) - F(XW^{(k)})$, i.e., (DCA) is a descent method (without linesearch). This fact can be seen from the fact that (DCA) fits into the majorization-minimization framework [2].\n- If $G(W^{(k+1)}) - F(XW^{(k+1)}) = G(W^{(k)}) - F(XW^{(k)})$, then (DCA) terminates at the $k$th iteration and $W^{(k)}$ is a critical point of (P). Here, a critical point $\\bar{W}$ of (P) is defined as a point satisfying $\\partial G(\\bar{W}) \\cap X^\\top \\partial F(X\\bar{W}) \\neq \\emptyset$, which is a necessary condition for optimality. If $G(W) - F(XW)$ is bounded below, then every limit point of the sequence $\\{W^{(k)}\\}$ is a critical point. We note that the optimal value being bounded below is an assumption that is trivially satisfied for all our PCA formulations since the classical formulation optimizes a continuous function over a compact set, and the optimal values of the other formulations can be related through strong duality or problem transformations.\n- In general, the rate for (DCA) is sublinear (also known as linear in the complexity theory literature). However, if certain growth conditions hold, then linear rates are possible, see for example [3].\n\nMoreover, from the described connection between some of our algorithms and simultaneous iteration (Theorem 3.12), we inherit the linear convergence rate of simultaneous iteration [4, Theorem 28.4]. This insight will be added in the paragraph following Theorem 3.12.\n\n> The empirical study does not convincingly demonstrate the superiority of the proposed DC methods over classical PCA solvers. In Table 1, DCA-based algorithms (Algorithms 1–6) often take longer than standard eigensolvers, such as SVDS and KrylovKit, with no clear accuracy benefit. The best-performing formulations are still dominated by classical methods in terms of both speed and simplicity, which weakens the practical motivation for these DC-based approaches.\n\nIt is true that some of the DC methods do not perform better than the classical PCA solvers. This is to be expected since we discovered that these are related to simultaneous iteration, which is the default eigenvalue solver for finding *complete* eigendecompositions (the QR algorithm or more generally QZ algorithm). Since the setting we consider only cares about the top principal components, the specialized algorithms which extract few components are naturally better. We have merely included these algorithms as a baseline comparison and to show their (unexpected) connection to numerical linear algebra. These methods can be thought of as the “gradient descent” for DC problems and there are better algorithms that can be derived, such as ADMM-type algorithms [5] as well as those inspired by the connection between (DCA) and Frank-Wolfe [6].\n\nWe respectfully disagree regarding the assessment that the best-performing formulations are dominated by the standard eigensolvers in both speed and simplicity. As for speed, both Table 1 and the additional experiments in the Appendix show that ZeroFPR and PG perform better than the classical methods for a certain range of required accuracies. Regarding the simplicity, we wish to emphasize that we compare generic first-order methods against tailored state-of-the-art implementations that include many low-level optimizations. Hence the main takeaway we hope to convey through our preliminary experiments is that for the level of accuracy typically sufficient in machine learning applications, there exist formulations for which faster solvers can be designed.\n\n> The numerical experiments are limited to synthetic Gaussian data with a fixed dimensionality and target rank. There is no demonstration on real-world datasets, and no evaluation of approximation quality, numerical stability, or robustness to noise.\n\nWe would like to emphasize that the main contribution of this paper is its theory, which all the reviewers have acknowledged to be sound. In that sense, we believe that a thorough empirical study should be deferred to a separate paper. Moreover, our experiments show that generic first-order optimization methods can beat state-of-the-art numerical linear algebra routines in the “machine learning regime”, which is an unexpected observation that goes against conventional wisdom. Regarding numerical stability, since our formulations are based on Schatten-norms, their implementations require (small) SVDs, which operate through unitary matrices and are therefore numerically stable. See our response to **Reviewer BYYV** for additional experiments.\n\n---\n\nIn the following we answer the posed **Questions**.\n\n> Include convergence guarantees or at least local convergence conditions for each of the proposed DCAs.\n\nSee our response to the first **Weakness** as well as our response to the second question of **Reviewer jssA**.\n\n> Benchmark against methods for sparse/robust PCA on real datasets (e.g., MNIST, Yale faces, motion segmentation).\n\nWe would like to point out that classical robust PCA formulations [7] are of a different nature than our proposed method. The classical methods often decompose a matrix into a low-rank component and a sparse component, which are well-suited for motion segmentation. However, our robust formulation views PCA from its autoencoder perspective and aims to enforce sparsity on the reconstruction error, which is better aligned with outlier detection (note that this is also different from sparse PCA). Another major difference is that most classical methods are not out-of-sample applicable while ours is through Theorem 3.9. These differences make fair comparisons difficult.\n\nTo further illustrate this fact, consider the synthetic experiment from [7, Section 4.1] (note that they do not present any other numerical results with metrics that can be easily compared). We first generate a low-rank component $L^\\star$ according to the same settings and add a sparse matrix $E$ where $E$ has limited support and the entries of $E$ are independent Bernoulli $\\pm 1$ entries. We consider two settings: \n1. the support of $E$ is chosen uniformly distributed over all its entries.\n2. the support of $E$ consists exclusively of full rows (i.e., an image where complete rows are perturbed). \n\nWe then compare the classical principal component pursuit (PCP) with alternating directions [7, Algorithm 1] and our Algorithm 5. For our experiment, we choose $\\\\mathrm{rank}(L^\\star) = 25$ and $L^\\star\\in\\mathbb{R}^{500\\times 500}$. We assume the size of the support of $E$ is $|E| = 12500$ and compute both the reconstruction error of the low-rank component and the primal cost of Proposition 4.1 (i.e., the $l_1$-norm of the row-wise reconstruction errors) and summarize our result in the following table (all metrics were averaged over $20$ different random problems).\n\n|           | $\\\\\\|L - L^\\\\star\\\\\\|_F / \\\\\\|L^\\\\star\\\\\\|_F$ PCP | $\\\\\\|L - L^\\\\star\\\\\\|_F / \\\\\\|L^\\\\star\\\\\\|_F$ Alg 5 | Primal cost PCP | Primal cost Alg 5 |\n| --------- | ------------------------------------------------- | --------------------------------------------------- | ------------------------ | -------------------------- |\n| Setting 1 | $3 \\\\times 10^{-8}$                               | 0.013                                               | 2488                     | 2364                      |\n| Setting 2 | 0.045                                             | 0.009                                               | 559                      | 545                      |\n\n\nWe observe that the classical algorithm from [7] performs better in the setting where it was proposed while our algorithm is better in the outlier setting. Moreover, we see that the primal cost function is lower for our algorithm in both cases, which is to be expected since the PCP algorithm is not designed to minimize this cost.\n\n---\n\nReferences\n\n[1] Tao, Pham Dinh, and LT Hoai An. Convex analysis approach to DC programming: theory, algorithms and applications. Acta mathematica vietnamica 22(1) (1997): 289-355.\n\n[2] Sun, Ying, Prabhu Babu, and Daniel P. Palomar. Majorization-minimization algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing 65(3) (2016): 794-816.\n\n[3] LT, Hoai An, Van Ngai Huynh, and Pham Dinh Tao. Convergence analysis of difference-of-convex algorithm with subanalytic data. Journal of Optimization Theory and Applications 179(1) (2018): 103-126.\n\n[4] Trefethen, Lloyd N., and David Bau. Numerical linear algebra. Society for Industrial and Applied Mathematics, (2022).\n\n[5] Sun, Tao, et al. Alternating direction method of multipliers with difference of convex functions. Advances in Computational Mathematics 44(3) (2018): 723-744.\n\n[6] Yurtsever, Alp, and Suvrit Sra. CCCP is Frank-Wolfe in disguise. Advances in Neural Information Processing Systems 35 (2022): 35352-35364.\n\n[7] Candès, Emmanuel J., et al. Robust principal component analysis?. Journal of the ACM (JACM) 58(3) (2011): 1-37."}, {"review_id": "a3dKrsiMVx", "reviewer": "Reviewer_jssA", "summary": "This paper revisits Principal Component Analysis (PCA) through the framework of difference-of-convex (DC) duality, offering new theoretical insights and algorithmic developments. The authors establish a novel connection between PCA and DC optimization, showing that applying the Difference-of-Convex Algorithm (DCA) to variance maximization is closely related to simultaneous iteration, a classical numerical linear algebra technique akin to the QR algorithm. This reveals an unexpected link between optimization and traditional matrix computations. Additionally, the paper introduces a kernelizable dual formulation for PCA that supports out-of-sample extensions, making it suitable for modern machine learning applications. A key contribution is a new robust variant of kernel PCA that minimizes the $l_1$-norm deviation of reconstruction errors, along with an associated algorithm connected to iteratively reweighted least squares. Empirical results demonstrate that simple first-order optimization methods can outperform state-of-the-art solvers for certain formulations, depending on the required accuracy. The work opens several promising research directions, including exploring other unitarily invariant objectives, developing deep stacked (kernel) PCA architectures, and designing specialized algorithms for the proposed DC dual pairs. Overall, this study provides a fresh optimization-based perspective on PCA, bridging classical techniques with modern machine learning applications.", "questions": "1. Proposition 3.3 appears to rely on the distinctness of singular values. Could this assumption be relaxed while still preserving the validity of the proposition? Alternatively, does the non-distinct case lead to fundamentally different results?\n2. The paper introduces several novel optimization-based algorithms for PCA. Could you provide a more rigorous discussion of their convergence properties? Specifically: Under what conditions (e.g., step sizes, initialization) do they converge? Are there local vs. global convergence guarantees? How does their convergence rate compare to classical methods? For the robust ℓ₁ variant, does the connection to IRLS imply any inherited convergence guarantees?\n3. While Assumption 3.1 appears natural from a data processing standpoint, it represents a significant restriction from the perspective of numerical linear algebra where such constraints aren't typically required. If we relax this assumption, which of the proposed PCA reformulations remain valid?", "limitations": "yes", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths:\nThe paper provides a fresh perspective on PCA by framing it within difference-of-convex (DC) duality, uncovering new connections between optimization and classical linear algebra. The link between the DC algorithm (DCA) and simultaneous iteration (related to the QR algorithm) is particularly insightful, bridging optimization and numerical methods in an elegant way. The introduction of a kernelizable dual formulation with out-of-sample extensions makes the approach applicable to modern machine learning problems. A notable contribution is the robust kernel PCA variant based on ℓ₁ reconstruction error minimization, which enhances robustness to outliers. The connection to iteratively reweighted least squares (IRLS) provides a well-founded optimization approach for this problem. By connecting PCA, DC optimization, and numerical linear algebra, the paper helps unify perspectives from different domains, making it valuable for researchers in optimization, linear algebra, and machine learning.\n\nWeaknesses:\nWhile the connection between DCA and simultaneous iteration is intriguing, the paper does not provide explicit convergence rates or conditions under which the proposed algorithms are guaranteed to converge to the global minimum. Additionaly, for the robust $l_1$ PCA variant, the relationship to IRLS is mentioned, but formal convergence properties are not rigorously established.", "was_revised": true, "final_justification": "I appreciate the theoretical contributions of the paper and continue to support it. My rating remains unchanged.", "rebuttal": "We thank the reviewer for their high evaluation, insightful questions and useful recommendations to improve our work.\n\n---\n\nIn the following we answer the posed **Questions**.\n\n> Proposition 3.3 appears to rely on the distinctness of singular values. Could this assumption be relaxed while still preserving the validity of the proposition? Alternatively, does the non-distinct case lead to fundamentally different results?\n\nThe distinctness of singular values can be relaxed. For the proof of Proposition 3.3, it suffices to change scalars to scalar matrices, i.e. $\\sigma_i$ becomes $\\sigma_i I_{n_i}$, where $n_i$ denotes the multiplicity. It follows that Proposition 3.3 stays largely the same, with the only difference being that multiplicities should be taken into account and columns of $V$ may be permuted within the space of right singular vectors corresponding to the same singular value. We will remove this assumption in the revised manuscript.\n\n We note that if this assumption does not hold, then the performance of classical methods degrades considerably (recall that the standard power iteration converges linearly with rate $\\lambda_2/\\lambda_1$ where the eigenvalues $\\lambda_i$ of a symmetric matrix are sorted in decreasing order), see for example Table 7 where we consider spectra with small spectral gaps. Empirically, we therefore observe that the first-order methods on these formulations are less sensitive to the spectral properties. However, we would also like to note that this case does not occur in practice.\n\n> The paper introduces several novel optimization-based algorithms for PCA. Could you provide a more rigorous discussion of their convergence properties? Specifically: Under what conditions (e.g., step sizes, initialization) do they converge? Are there local vs. global convergence guarantees? How does their convergence rate compare to classical methods? For the robust ℓ₁ variant, does the connection to IRLS imply any inherited convergence guarantees?\n\nSee our response to **Reviewer FsFQ** for convergence guarantees for the DCA algorithms. Regarding the convergence of PG and zeroFPR, the standard assumptions to minimize objectives of the form $f(x) + g(x)$ using these methods require $f$ to be $L$-smooth, while $g$ is often the indicator of a closed and convex set as in our case. In general, the (global) rates of PG and zeroFPR to stationary points are both sublinear $O(1/k)$ [1]. We note that, strictly speaking, for some formulations our $f$ is nonsmooth. However, this does not pose a problem since the objectives are concave. Therefore, the Euclidean descent lemma [2, Lemma 5.7] holds with $L=0$, and both PG and zeroFPR enjoy the same theoretical guarantees as described above. In fact, we may choose arbitrarily large stepsizes, despite nonsmoothness. This discussion will be added in the revised manuscript.\n\nThe general rate for IRLS algorithms is also sublinear. However, under the restricted isometry property, it is possible to show a linear rate for IRLS [3]. While we have shown a connection between Algorithm 9 (DCA) and Algorithm 10 (IRLS), it is not immediately clear how to concretely map the iterates to each other so a linear rate in this setting is not immediately inherited. We conjecture that such a rate is also possible for our Algorithm 9 under the same property, though we believe that such an endeavor deserves separate work.\n\n> While Assumption 3.1 appears natural from a data processing standpoint, it represents a significant restriction from the perspective of numerical linear algebra where such constraints aren't typically required. If we relax this assumption, which of the proposed PCA reformulations remain valid?\n\nAssumption 3.1 is nowhere used in any of the derivations of alternate formulations, so all formulations remain valid. We have merely included this assumption so that the statistical interpretation of PCA is still valid (i.e. variance maximization). We will remove this assumption and only mention it when speaking about the statistical interpretation to make it clear that this is not a restriction of the paper.\n\n---\n\nReferences:\n\n[1] Themelis, Andreas, Lorenzo Stella, and Panagiotis Patrinos. Forward-backward envelope for the sum of two nonconvex functions: Further properties and nonmonotone linesearch algorithms. SIAM Journal on Optimization 28(3) (2018): 2274-2303.\n\n[2] Beck, Amir. First-order methods in optimization. Society for Industrial and Applied Mathematics, (2017).\n\n[3] Daubechies, Ingrid, et al. Iteratively reweighted least squares minimization for sparse recovery. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 63(1) (2010): 1-38."}, {"review_id": "rRpdUhr6I9", "reviewer": "Reviewer_A5aE", "summary": "The paper aims at proposing a faster algorithm in solving PCA.", "questions": "could the authors comment on practical use  for the proposed algorithms?", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**strengths:**\n- thorough mathematical presentation of the proposed algorithms  \n\n**weaknesses:**\n- paper difficult to follow in particular its purpose, goal and its  contributions e.g.  what are practical uses and potential impact of the proposed algorithm(s) for cited applications (e.g. transformer [64])?\n- lines 108/176/311 missing Fig 1 and/or Fig 2. reference in the sentences ?\n- Appendices are missing (?)\n- Table 1 shows experimental results for computational times but it is difficult to interpret as different units are used (e.g. second, milli seconds), and it is unclear if they are statistically different.\n- Tab1 considers three problem settings with variations only on variable $N$ (number of samples) and $d$ (feature dimension) while $s$ and $\\epsilon$ remains the same. Would plots be  more suited to highlight impacts of $N$ and $d$ on compute  time?", "was_revised": true, "final_justification": "Thanks for authors for their rebuttal.\nI have updated my evaluation.", "rebuttal": "We thank the reviewer for their time and comments.\n\n---\n\nIn the following we address the stated **Weaknesses**.\n\n> paper difficult to follow in particular its purpose, goal and its contributions e.g. what are practical uses and potential impact of the proposed algorithm(s) for cited applications (e.g. transformer [64])?\n\nWe recognize that the paper may not be the easiest to follow since the main contributions of the paper are of theoretical nature. However, we have motivated all our results clearly and made sure to state all the assumptions explicitly, while most mathematical details are deferred to the appendix. The main contributions are also stated in a separate section.\n\nWe reiterate the main goals/contributions of this paper:\n- Developing novel (and useful) formulations for PCA. Note that it is quite remarkable that for such a well-studied and fundamental problem as PCA, there still exist problem formulations that have not been described before.\n- Discovering new connections between PCA, optimization, and numerical linear algebra.\n- Deriving results for a family of PCA like formulations that are kernelizable and out-of-sample applicable.\n- As an example of this last result, we consider a robust PCA formulation that minimizes the $l_1$-deviation of the reconstruction errors, thereby making the principal components less sensitive to outliers.\n\nSome examples of immediate practical applications of this paper are:\n- Faster PCA solvers for the accuracy required in machine learning applications. Note that PCA is still prevalent in many domains, e.g., in bio-informatics.\n- A way to robustly extract features in a kernel setting, with a strong theoretical foundation.\n\nThe most promising research avenues that we have in mind for the future are:\n- Tailored algorithms for solving the novel problem formulations, thereby obtaining even faster solvers for PCA. Moreover, from the connection between kernel PCA and self-attention [1], these algorithms may inspire better algorithms for training LLMs.\n- Novel (deep) architectures inspired by the new formulations (see also our response to **Reviewer BYYV**).\n- Deriving novel theoretical results by utilizing the connection between the various scientific fields. In particular, since simultaneous iteration is known to converge linearly from numerical linear algebra [2, Theorem 28.4], this means that the corresponding difference-of-convex formulation must satisfy some growth condition. By utilizing this connection, we believe it is possible to reverse-engineer this condition, which may subsequently also hold for many other formulations, and so on.\n\n> lines 108/176/311 missing Fig 1 and/or Fig 2. reference in the sentences ?\n\nWhile the tags are properly hyperlinked to the corresponding location, we will make this more clear by giving hyperlinks a distinct color as well as mentioning the explicit Figure numbers.\n\n> Appendices are missing (?)\n\nThe appendices are present in the supplementary material, as is allowed by the NeurIPS formatting guidelines: \"Technical appendices with additional results, figures, graphs and proofs may be submitted with the paper submission before the full submission deadline (see above), or as a separate PDF in the ZIP file below before the supplementary material deadline.\"\n\n> Table 1 shows experimental results for computational times but it is difficult to interpret as different units are used (e.g. second, milli seconds), and it is unclear if they are statistically different.\n\nWe will report all computation times in milliseconds to avoid ambiguity, as has already been done in the supplied Appendix. As for the statistical differences of the timings, we emphasize that we use BenchmarkTools.jl [3], which, to quote: \"is robust in the presence of timer errors, OS jitter and other environmental fluctuations, and is insensitive to the highly nonideal statistics produced by timing measurements\", thereby providing timing data that is more statistically reliable than what is typically encountered in many machine learning papers. In particular, we draw attention to Tables 6 and 7, for which it should be clear that the timings between ZeroFPR/PG and SVDS/KrylovKit are statistically different.\n\n> Tab1 considers three problem settings with variations only on variable $N$ (number of samples) and $d$ (feature dimension) while $s$ and $\\varepsilon$ remains the same. Would plots be more suited to highlight impacts of $N$ and $d$ on compute time?\n\nWe thank the reviewer for this recommendation. We have avoided plots for these results since they are quite cluttered. However, if space permits, we will include visualizations of only the best performing algorithms to better highlight the impacts of $N$ and $d$. See our response to **Reviewer BYYV** for some additional experiments with varying $s$ and see Table 5 in the appendix for an additional experiment with different $\\varepsilon$.\n\n---\n\nRegarding the posed **Question** concerning the practical use, see our response to the first **Weakness**.\n\n---\n\nWe would like to kindly request the reviewer to reconsider their score based on the other reviews and our rebuttal, since it is not clear how the stated weaknesses merit the current score.\n\n---\n\nReferences\n\n[1] Teo, Rachel SY, and Tan Nguyen. Unveiling the hidden structure of self-attention via kernel principal component analysis. Advances in Neural Information Processing Systems 37 (2024): 101393-101427.\n\n[2] Trefethen, Lloyd N., and David Bau. Numerical linear algebra. Society for Industrial and Applied Mathematics, (2022).\n\n[3] Chen, Jiahao, and Jarrett Revels. Robust benchmarking in noisy environments. arXiv preprint arXiv:1608.04295 (2016)."}]

Author rebuttal:

Meta-review: The submission studies the classical PCA problem from the perspective of difference-of-convex (DC) optimization. While DC optimization is known to be generally intractable, some heuristics tend to perform fairly well on practical instances. This is a situation that resembles some numerical lineal algebra problems, such as nonnegative matrix factorization or sparse PCA. Hence, I find this approach natural and insightful.

Reviews agree on the technical merits of this work, but some have concerns about the experimental evaluations. Authors provided further results in this line; while I find this not strictly necessary, it is appreciated that the authors provide more thorough evaluations. 

Based primarily on the technical contributions, this work seems suitable for publication.
Decision: Accept (poster)