{'fraction_non_alphanumeric': 0.04453525521955569, 'fraction_numerical': 0.02043242440002972, 'mean_word_length': 4.496691987258025, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse PCA With Multiple Components

Sparse Principal Component Analysis is a cardinal technique for obtaining combinations of features, or principal components (PCs), that explain the variance of high-dimensional datasets in an interpretable manner. At its heart, this involves solving a sparsity and orthogonality constrained convex maximization problem, which is extremely computationally challenging. Most existing work address sparse PCA via heuristics such as iteratively computing one sparse PC and deflating the covariance matrix, which does not guarantee the orthogonality, let alone the optimality, of the resulting solution. We challenge this status by reformulating the orthogonality conditions as rank constraints and optimizing over the sparsity and rank constraints simultaneously. We design tight semidefinite relaxations and propose tractable second-order cone versions of these relaxations which supply high-quality upper bounds. We also design valid second-order cone inequalities which hold when each PC's individual sparsity is specified, and demonstrate that these inequalities tighten our relaxations significantly. Moreover, we propose exact methods and rounding mechanisms that exploit these relaxations' tightness to obtain solutions with a bound gap on the order of 1%-5% for real-world datasets with p = 100s or 1000s of features and r \in {2, 3} components. We investigate the performance of our methods in spiked covariance settings and demonstrate that simultaneously considering the orthogonality and sparsity constraints leads to improvements in the Area Under the ROC curve of 2%-8% compared to state-of-the-art deflation methods. All in all, our approach solves sparse PCA problems with multiple components to certifiable (near) optimality in a practically tractable fashion.


Introduction
The ever-increasing amount of data generated and collected in the past decades has fueled a growing interest in machine learning (ML) algorithms that scale to ultra-high-dimensional settings.
On the other hand, regulatory requirements such as the right to explanation enshrined in the European Union's General Data Protection Regulation (Goodman and Flaxman 2017) promote algorithms with outputs that are interpretable to humans. However, since most ML algorithms generate uninterpretable learners (Rudin et al. 2022), a key challenge for data scientists is to generate interpretable learners that perform well in high-dimensional settings. A standard technique for achieving this is dimensionality reduction, wherein we reduce a p-dimensional dataset to an r-dimensional one (with r p) with minimal information loss.
Perhaps the most popular dimensionality reduction technique is Principal Component Analysis, or PCA, initially proposed by Pearson (1901) (see also Hotelling 1933, Eckart andYoung 1936). Given a normalized and centered data matrix A ∈ R n×p , and its sample covariance matrix Σ := 1 n−1 AA , practitioners identify the top r principal components of Σ by solving: and subsequently project the data matrix A onto the principal components U . PCA is now a cardinal unsupervised learning paradigm that is practically useful across a range of fields, including pattern recognition (Naikal et al. 2011) and sequence classification (Tan et al. 2014).
Unfortunately, while modern implementations of PCA are extremely time-efficient (Udell et al. 2016, Tropp et al. 2017, PCA generates components that are dense, error-prone, and statistically meaningless in high-dimensional settings due to the curse of dimensionality (Amini and Wainwright 2008, Berthet and Rigollet 2013). Accordingly, several authors (such as Jolliffe et al. 2003, d'Aspremont et al. 2007) have proposed sparse PCA, namely augmenting Problem (1) with a sparsity constraint. When r = 1, sparse PCA can be formulated as the following optimization problem: max u∈R p Σ, uu s.t. u 2 2 = 1, u 0 ≤ k, where u 0 denotes the cardinality of u or the size of its support: u 0 = | supp(u)| = |{j : u j = 0}|.
Problem (2) can be formulated as a mixed-integer semidefinite optimization problem and solved via global optimization techniques like branch-and-cut or convex optimization techniques like relaxand-round (Bertsimas et al. 2022, Li andXie 2020). High-quality solutions can also be obtained by greedy heuristics (d'Aspremont et al. 2008) or solving its 1 relaxation (Zou et al. 2006, d'Aspremont et al. 2007, Dey et al. 2022a. Unfortunately, no consensual formulation extends Problem (2)'s formulation to r > 1, and there are no practically relevant algorithms with optimality guarantees that successfully address this extension. This is a significant issue for two reasons.
First, from a statistical recovery perspective, sparse PCA performs provably better than PCA in some rigorously-identified settings (see Gamarnik et al. 2021, for further details). Empirically, it is believed to generally outperform PCA as well (d'Aspremont et al. 2007, Gamarnik 2021). Furthermore, among sparse PCA methods, there exists an information-theoretic performance gap between certifiably optimal methods and heuristic approaches (Berthet andRigollet 2013, Gamarnik 2021).
Second, certifiably optimal methods for sparse PCA are also valuable from a practical perspective. Indeed, for r = 1, sparse PCA methods which do not possess optimality guarantees are typically and often significantly less accurate than exact methods. For example, Berk and Bertsimas (2019) compared the performance of seven popular sparse PCA heuristics across four UCI datasets (k ∈ {5, 10}) and observed that no heuristic found an optimal solution across all datasets, with some methods being routinely 20% suboptimal or more. On the other hand, recent exact methods for r = 1, including branch-and-bound (Berk and Bertsimas 2019) and branch-and-cut  found the optimal solution in all cases on the same benchmark.
In this paper, we propose a generic optimization formulation that extends Problem (2) to r > 1, reformulate it as a sparsity and rank constrained optimization problem, and design certifiably (near) optimal techniques to solve it. As in the r = 1 case, we demonstrate numerically that exact methods obtain substantially higher quality solutions than existing state-of-the-art heuristics.

A Generic Formulation for Sparse PCA with Multiple PCs
Perhaps the most natural extension of sparse PCA to multiple principal components, and the one which we advocate in this paper, is to augment Problem (1) with a constraint on the number of non-zero entries in the matrix U , This gives a formulation which enforces two desirable properties on the matrix U : orthogonality (U U = I), as present in the prototypical formulation of PCA; and sparsity, to address interpretability and accuracy concerns.
Formally, introducing a binary matrix Z to encode the support of U , we consider the problem: where pr > k > r in order that the problem is well-posed and non-trivial.
Alternatively, instead of imposing an overall sparsity budget (k) on the entire matrix U , we could restrict the size of the support of each column of U separately, i.e., impose U t 0 ≤ k t ∀t ∈ [r] and adapt (3) accordingly. Indeed, optimizing this formulation over all integer combinations k t which sum to k is equivalent to solving (3). We address both modeling options in this paper. From a generative model perspective, Problem (3) is consistent with a spiked covariance model (see, e.g., Amini andWainwright 2008, d'Aspremont et al. 2008), where the true covariance matrix Σ can be decomposed as the sum of a sparse and low-rank term plus some noise: While sparse PCA with one PC has received much attention over the past 20 years, sparse PCA with multiple PCs has received substantially less attention, arguably because of the significant technical difficulty of simultaneously optimizing over orthogonality and sparsity constraints. Indeed, in the single PC case, various authors including Gally and Pfetsch (2016), , Li and Xie (2020), Kim et al. (2021), Dey et al. (2022a) have shown that Problem (2) can be recast as a mixed-integer semidefinite optimization (MISDO) problem and have derived both high-quality solutions and valid dual bounds. For the r > 1 case however, we are not aware of any equivalent results-except when very strong assumptions are made about the sparsity pattern, as discussed in Section 1.2. Moreover, as established in our numerical results (Section 4), existing heuristics routinely terminate with highly suboptimal solutions. For example, on the pitprops dataset, with 6 PCs, we explain 81% of the variance with an overall sparsity of 24 while previous studies could only explain less than 70% of the variance with twice times as many variables (see Section 4.2).

Literature Review
To identify the extent to which the state-of-the-art for sparse PCA could be improved, we now review methods that have been proposed to approximately solve (3) or its special cases.
Methods for generic sparse PCA: It is well known that an optimal solution to Problem (1) can be obtained via a greedy procedure where, at each iteration, the leading PC of Σ, u, is computed (e.g., by solving (1) with r = 1) and then the matrix Σ is updated (or deflated ) to annihilate u. This so-called deflation method provides an optimal solution to the PCA problem with multiple PCs (see, e.g., Mackey 2008). Consequently, Mackey (2008), building upon the work of d 'Aspremont et al. (2007), propose a sparse extension of this deflation scheme where, at each iteration, a sparse PC is computed (e.g., by solving (2) or a relaxation) and Σ is updated by projecting out the eigenspace modeled by u: Σ new = (I − uu )Σ(I − uu ). Empirically, this method often performs reasonably well (see Berk and Bertsimas 2019, Section 5.3), particularly when Problem (2) is solved to global optimality; see also Hein and Bühler (2010), Bühler (2014) for a related deflation-based scheme. However, unlike in the traditional case, deflation need not return an optimal solution to (3). Actually, deflation need not even return a feasible solution, since the orthogonality constraint is not explicitly imposed by the method, and is therefore often violated in practice. Indeed, because of their iterative nature, deflation-based procedures can easily control the sparsity of each component but usually struggle to enforce the orthogonality constraint.
A second approach for solving Problem (3) is to apply a heuristic which (approximately) optimizes all r PCs simultaneously, rather than sequentially. Among others, Zou et al. (2006) propose an alternating minimization scheme for an 1 relaxation of Problem ( (3) approximately. Unfortunately, these approaches are often suboptimal if r = 1, and only provide candidate solutions, with no indication on the optimality gap. Indeed, none of these approaches explicitly control both the sparsity and orthogonality constraints, and therefore none of the methods reviewed here are actually guaranteed to return feasible solutions to (3).
Row-Sparsity: Motivated by tractability concerns, another line of work studies a special case of (3), namely row-sparse principal component analysis or principal component analysis with global support. In this formulation, we replace the sparsity constraint on U with one requiring U has at most k non-zero rows, as advocated by Boutsidis et al. (2011), Probel andTropp (2011), . This implies we can rewrite sparse PCA as performing a top-r SVD on a subset of k rows of Σ, i.e., solving Problem (5) is a special case of Problem (3) where each PC has the same support, which corresponds to constraining Z such that z i,t = z j,t ∀t ∈ [r] in (3). This restriction is advantageous from a computational perspective, but disadvantageous from a statistical one. Indeed, because of the global support assumption, Problem (5) can be reformulated as a mixed-integer semidefinite problem (see Bertsimas et al. 2022, Li andXie 2021, Bertsimas and Kitane 2022, for derivations) and solved via relax-and-round (Dey et al. 2022b, Li and Xie 2021) or branch-and-bound (Del Pia 2022 strategies. However, from a generative model perspective (4), it is equivalent to making the very strong assumption that all leading eigenvectors u t have the same sparsity pattern. Indeed, if one PC u t is sparser than k ( u t 0 ≤ k t < k) or two PCs are partially disjoint then Problem (5) will nonetheless estimate that the support of each u t is of size k. Therefore, Problem (5) Note that (6) is a special case of (3) where we additionally require that t∈[r] Z i,t ≤ 1 ∀i ∈ [p], i.e., that each PC is supported by non-overlapping features. Interestingly, this restriction allows (6) to be recast as a MISDO and solved as such (c.f. Kitane 2022).
An obvious criticism of (6)'s formulation is that it only allows each feature to be assigned to one PC, while in practice we may wish to assign a feature to multiple PCs. Therefore, (6) is best thought of as a special case of sparse PCA -and of the generative model (4). Indeed, if the true generative model involves vectors u t with partially overlapping supports, then (6) could not recover them. Nonetheless, as we explore in our numerical experiments (Section 4), disjoint solutions often perform well when k is very small relative to p. This is because, in high-dimensional settings, there are often several disjoint submatrices which are near-optimal in the rank-one case, and selecting the leading PCs from each of them is often a reasonable approach in practice.

Contributions and Structure
To our knowledge, no existing algorithm can solve sparse PCA problems with multiple components and obtain certificates of (near) optimality, except in the very special cases of row-sparsity or disjoint support. Accordingly, in this paper, we undertake a detailed study of Problem (3).
The main contributions of the paper are twofold. First, we reformulate Problem (3) as a mixedinteger low-rank problem. Based on this reformulation, we derive tight, yet tractable, semidefinite and second-order cone relaxations of this formulation, and propose valid inequalities which strengthen the relaxation. Second, from these convex relaxations, we propose rounding techniques and an alternating minimization scheme to obtain high-quality solutions. By combining strong relaxations and good solutions, we obtain bound gaps of 5% on real-world instances of sparse PCA with multiple PCs where p = 100s or 1000s. All in all, we provide the first reformulation of the sparse PCA problem (3) that is amenable to certifiably optimal algorithms and admits tractable relaxations, while demonstrating its benefit compared with state-of-the-art heuristics.
The structure and detailed contributions of the paper are as follows: • In Section 2.1, we prove that rank constraints successfully model orthogonality constraints. Leveraging this result, we reformulate Problem (3) as a sparsity and rank constrained optimization problem, and derive its semidefinite relaxation. We propose valid inequalities which strengthen the relaxation in Section 2.2, and additional valid inequalities which hold if we restrict the support of each PC separately in Section 2.3. For high-dimensional instances, we propose a tractable second-order cone relaxation in Section 2.4.
• In Section 3, we develop numerical strategies to obtain high-quality feasible solutions to Problem (3). In Section 3.1, we develop a greedy rounding strategy which rounds an optimal solution to one of Section 2's relaxations into r sparse PCs with disjoint support (hence, orthogonal). In Section 3.2, we develop an exact formulation of Problem (3) in the original space of decision variables, which can be addressed via global branch-and-bound solvers, and propose valid inequalities which improve the tightness of the formulation. Finally, in Section 3.3, we consider a Lagrangean relaxation of Problem (3) and design an alternating minimization strategy to obtain near-orthogonal solutions by iteratively solving a sequence of sparse PCA problems with r = 1. Compared with existing deflation-based techniques, we explicitly penalize the orthogonality violation at each iteration, with an increasing penalty parameter, thus converging towards an orthogonal and typically high-quality solution.
• In Section 4, we investigate the quality of our semidefinite relaxations and feasible methods and observe empirically that they collectively give small bound gaps in practice. First, we invoke our convex relaxations to obtain high-quality upper bounds on UCI datasets in Section 4.1. We also compare our approximate methods with existing methods from the literature, including the methods of Zou et al. (2006), Hein and Bühler (2010), Berk and Bertsimas (2019). On both synthetic and UCI datasets, we demonstrate in Sections 4.2-4.3 that our approaches, particularly the alternating minimization scheme developed in Section 3.3, explain more of the variance in the data and exhibit a lower false discovery/higher true discovery rate for a given sparsity budget.
Finally, in Section 4.4, we explore the relationship between the symmetry in the size of each PCs sparsity budget and the amount of variance explained.

Preliminaries and Notation
We let nonbold face characters such as u denote scalars, lowercase bold faced characters such as u denote vectors, uppercase bold faced characters such as U denote matrices, and calligraphic uppercase characters such as Z denote sets. If U is a matrix then U t denotes the tth column vector of U , and U i,t denotes the (i, t)th entry of U . We let [p] denote the set of running indices {1, ..., p}.
We let e denote a vector of all 1's, 0 denote a vector of all 0's, and I denote the identity matrix, with dimension implied by the context.
We also use an assortment of matrix operators. We let ·, · denote the Euclidean inner product between two matrices, · F denote the Frobenius norm, and S p + denote the p × p positive semidefinite cone; see Horn and Johnson (1990) for a general theory of matrix operators.
We repeatedly use some basic properties of orthogonal projection matrices, often without proof.
Let Y n := {Y ∈ S n : Y 2 = Y } denote the set of n × n orthogonal projection matrices and Y k n := {Y ∈ Y n : tr(Y ) ≤ k} denote projection matrices with rank at most k. Among others, the convex hulls of Y n and Y k n are well-studied, as we now remind the reader: is given by The above result is sometimes known as a perspective reformulation, since we strengthen the Lemma 3. (Theorem 1 of Atamtürk and Gomez 2019) The convex closure of the set is given by Lemma 3 is extremely useful when a single continuous variable depends upon multiple indicator variables, as occurs in certain substructures of our reformulations of Problem (3).

Exact Formulations and Their Relaxations
In this section, we reformulate Problem (3) as a mixed-integer low-rank problem, study its semidefinite and second-order cone relaxations, and propose valid inequalities for strengthening these relaxations. In conjunction with the methods for generating high-quality feasible solutions we develop in the next section, this allows us to, for the first time, solve instances of Problem (3) in a manner that maintains optimality guarantees at scale.

An Extended Formulation and Its Semidefinite Relaxation
Our sparse PCA formulation (3) is a mixed-integer quadratic optimization problem that exhibits three primary sources of difficulty. First, as typical in PCA problems, Problem (3) maximizes a convex function in the decision variable U . Second, there is a sparsity constraint, which is computationally challenging to model, although recent evidence suggests that sparsity constraints need not imply intractability (see, e.g., . Finally, and more consequentially, there is an orthogonality constraint. To the best of our knowledge, existing generic non-convex solvers such as Gurobi cannot optimize over such orthogonality constraints at scale. To address the three aforementioned difficulties, we now derive an orthogonality-free reformulation in five steps. First, we introduce the matrix Y = U U I, thus linearizing the non-convex objective. Second, we introduce rank-one matrices Y t to model the outer products of each column of U , U t , with itself, U t U t . Third, we reassign the indicator variable Z i,t to model whether Y t i,j , rather than U i,t , is non-zero. Fourth, by letting Y = r t=1 Y t , we observe that we can omit the matrix U (and the constraints involving U ) without altering the set of feasible Y 's. Finally, we use the fact that Y t i,j is only supported on indices i, t where Z i,t > 0 to strengthen the constraint Y I to Y Diag (min (e, t Z t )). Formally, we have: Theorem 1. Problem (3) attains the same optimal objective value as the problem: tr Remark 1. We do not explicitly require that Y 0, since Y is the sum of positive semidefinite matrices Y t . Indeed, omitting Y 0 substantially improves the tractability of (7)'s relaxations.
The proof of Theorem 1 requires an intermediate result (Proposition 1). Proposition 1 shows that by imposing rank-one constraints on each Y t , the condition that the Y t 's are mutually orthogonal can be reforumlated as a linear semidefinite constraint. Independently, Proposition 1 is crucial for designing our main algorithm in Section 3.3 (proof deferred to Section EC.1).
Proposition 1. Consider r matrices, Y t ∈ S p + , such that tr(Y t ) = 1 and Rank(Y t ) = 1. Then, Proof of Theorem 1 It suffices to show that for any feasible solution to (3), we can construct a feasible solution to Problem (7) with an equal or greater payoff, and vice versa.
• Let (U , Z) be a solution to Problem (3). Then, since U i,t can only be non-zero if Z i,t = 1 and U U I, it follows that U U Diag (min (e, t Z t )). Therefore, (Y : is a feasible solution to (7) with an equal cost.
• Let (Y , Y t , Z) denote a feasible solution to Problem (7). Then, since each Y t is symmetric and rank-one, we can decompose Y t as Y t = U t U t for a vector U t such that U i,t = 0 if Z i,t = 0, and concatenate these vectors U t into a matrix U such that (U , Z) has the same cost in (3) as (Y , Y t , Z) does in (7). Therefore, it remains to show that U U = I. To see this, observe that Theorem 1 provides a formulation which is less compact than (3), but contains rank constraints rather than orthogonality constraints. Therefore, it is amenable to exact approaches for addressing sparsity (Bertsimas et al. 2021b) and rank (Bertsimas et al. 2021a) constraints.
Furthermore, this formulation already provides valid upper bounds on (3)'s objective by relaxing the rank and sparsity constraints: where M i,i = 1 and M i,j = 1/2 if i = j is an upper bound on |Y t i,j |, since Y t was, before relaxing the rank constraint, a rank one matrix (c.f. ).

Valid Inequalities for Strengthening the Extended Formulation
In this section, we propose valid inequalities which allow us to improve the quality of the convex relaxation (8) introduced in the previous section. Formally, we have the following result: Theorem 2. Let P strengthened denote the optimal value of the following problem: where w i models min(1, r t=1 Z i,t ). Further, let P Relax denote the optimal value of Problem (8), and P Exact denote the optimal objective value of Problem (7). Then, for any covariance matrix Σ and any sparsity and rank (k, r): Proof of Theorem 2 The first inequality holds by verifying that any feasible solution to (7) is also a solution to (9) with equal cost.
follows by aggregating the 2 × 2 minors of Diag(w) − Y and invoking Lemma 3. We provide a full derivation of these valid inequalities in Section EC.2.
The second inequality holds by observing that (8) is a relaxation of (9).

Strong Inequalities with a Per-Component Sparsity Budget
In this section, we consider specifying a sparsity budget k t for each component Y t , in addition to an overall sparsity budget k (with k = r t=1 k t ), as originally proposed in Section 1.1. The primary motivation for this approach is that it allows us to derive further valid inequalities which strengthen Problem (9), sometimes substantially. Indeed, existing works which generate sparse and sometimes orthogonal PCs (c.f. Hein and Bühler 2010, Berk and Bertsimas 2019) and we benchmark against in our numerical experiments in Section 4 require that each k t is specified separately.
Formally, we have the following result (proof deferred to Section EC.3): (7). Then, the following inequalities hold: Interestingly, as we observe in Section 4.1, combining Problem (9) with Constraints (10)- (11) often yields much tighter upper bounds than (9) alone, even if we take the worst-case upper bound over all feasible splits k t 's which sum to k.
To further tighten the semidefinite relaxation, we now leverage additional problem structure.
For each component t ∈ [r], the feasible set of all k t -sparse components {u ∈ R p : u 0 ≤ k t } is permutation and sign invariant, i.e., for any feasible vector u, any vector obtained by permuting or changing the sign of the coordinates of u is also feasible. Based on this observation, Kim et al. (2021) recently propose a lifted formulation for sparse PCA with a single PC which contains additional valid inequalities. To the best of our knowledge, their formulation leads to the strongest existing relaxation for sparse PCA with r = 1, although it is less tractable than the relaxations developed in , Li and Xie (2020) and cannot scale beyond p = 100. Therefore, we now extend the valid inequalities in Kim et al. (2021)'s "T-relaxation" to sparse PCA problems with r > 1. To do so, for each t ∈ [r], we introduce an additional variable F t to capture the entry-wise absolute value of Y t , and an additional matrix G t which contains a sorted version of F t . We obtain: The additional variables r t,D and t t,D are introduced to enforce coupling constraints between the diagonal entries of F t and G t (Kim et al. 2021, eq. 44), while T t allows to couple F t with the binary variables Z (Kim et al. 2021, eq. 50). In contrast to Kim et al. (2021), we explicitly require that each F t is positive semidefinite (rather than that its 2 × 2 minors are), in order to obtain a stronger relaxation; we consider the 2 × 2 minors when developing a more tractable relaxation in the next section of the paper.
We now offer several remarks on the strength and generality of this formulation compared to the ones proposed in the previous section: Remark 2. Problem (12)'s relaxation requires that we have knowledge of k t , the sparsity of each individual PC, in addition to k = t∈[r] k t alone. Indeed, the constraint G t i,j = 0 if i > k t or j > k t is only linear if k t is known. Thus, (12)'s relaxation is less general than (9)'s.
Remark 3. Problem (12)'s relaxation dominates Problem (9)'s, even when (9) is strengthened with the inequalities proposed in (10)-(11). Indeed, Kim et al. (2021, Theorem 13) can be extended to show that (10)-(11) are redundant in (12)'s formulation. However, we did include in (12) the inequalities we introduced in Section 2.2, which capture notions of orthogonality between the Y t 's, since the lifted formulation based on Kim et al. (2021) applies to each Y t separately.
In our numerical experiments (see Table 1 in Section 4), we observe that (12) is strictly tighter than (9), at the price of a very modest amount of tractability. Therefore, after an initial comparison, we favor (12) throughout our numerical results. Unfortunately, (12) cannot scale beyond p = 100, at least with current technology, due to the presence of multiple semidefinite matrices and constraints.

Second-Order Cone Relaxations for High-Dimensional Settings
We now develop a more tractable, albeit less tight, version of (12)'s relaxation which scales to p > 100 features. Namely, we replace all semidefinite constraints of the form X ∈ S p + with the non-negativity of their 2 × 2 minors, , as proposed by Bertsimas and Cory-Wright (2020) and references therein. This gives a second-order cone relaxation of (12), (EC.5), which we lay out in Section EC.4 for brevity.
Remark 4. We can improve the second-order cone relaxation (EC.5) without compromising its tractability by iteratively solving (EC.5) and imposing linear cuts of the form X, vv ≥ 0, for each matrix X which is positive semidefinite in (12), where v is a trailing eigenvector of X in the most recent solution to the relaxation, as presented in Bertsimas and Cory-Wright (2020). In our numerical experiments (Section 4), we consider a version of this scheme where we iteratively impose one such cut corresponding to the most negative eigenvector in the matrices Y t , G t , F t , and Diag(w) − Y , and perform up to 50 iterations of this scheme, stopping early if the most negative eigenvalue in all constraints is −10 −4 or larger.

Algorithmic Strategies
In this section, we propose three numerical strategies which provide high-quality solutions to Problem (3). First, in Section 3.1, we propose a rounding mechanism which converts a solution to one of Section 2's convex relaxations into a certifiably near-optimal solution to Problem (3). Second, in Section 3.2, we propose a formulation of Problem (3) which is amenable to existing spatial branch-and-bound codes, and can be warm-started using the rounding mechanism proposed in Section 3.1. Finally, in Section 3.3, we propose an iterative deflation heuristic that exploits the relative maturity of sparse PCA technology in the rank-one case to obtain high-quality PCs in the rank-r case for r > 1. All in all, we provide a suite of algorithms that allow practitioners to address instances of Problem (3) in an accurate yet practically efficient fashion. Moreover, when combined with the upper bounds developed in Section 2, these strategies allow us to, for the first time, solve large-scale instances of Problem (3) to provable near-optimality (see Section 4).

Feasible Solutions from the Relaxation via Greedy Disjoint Rounding
In this section, we develop a rounding mechanism that converts an optimal solutions to a convex relaxation (see Section 2) into a high-quality feasible solution. Historically, a numerically useful strategy for similar integer optimization problems has been to (a) solve a convex relaxation in (Z, Y ), (b) greedily round Z , the solution to the relaxation, to obtain a feasible binary matrixẐ that is close to Z , and (c) resolve for U under the constraints U i,t = 0 ifẐ i,t = 0. For example, in our prior work , we established that a similar rounding strategy performs well for the single component case. In addition, such strategies are typically an important building block for exact approaches, such as spatial branch-and-bound, to obtain good upper and lower bounds for each subproblem, especially at the root node. Therefore, a natural question is whether rounding remains numerically useful in the presence of the greater technical difficulties induced by optimizing over multiple components simultaneously.
To answer this, observe that the "relax" and "round" steps are tractable for any r, since they respectively involve solving a convex relaxation and sorting. Therefore, the numerical tractability of greedy rounding hinges on the tractability of solving for U withẐ fixed. This is equivalent to: Unfortunately, (13) is very expensive to solve in high-dimensional settings, due to the orthogonality constraint U U = I. Indeed, in preliminary experiments on the UCI pitprops dataset (p = 13), spatial branch-and-bound routinely required more than 7200s to solve (13), and the solution found within the first 7200 seconds was usually significantly worse than solutions identified by other methods developed in this paper (even when the same support was identified). Indeed, the orthogonality constraints currently constitute the main bottleneck for spatial branch-and-bound codes like Gurobi, to the extent that, on the same pitprops instances, we observed that solving (13) was not significantly faster than jointly optimizing for Z and U . In order for our rounding scheme to be scalable and potentially useful for accelerating exact approaches, we avoid solving (13) directly.
To do so, we modify the rounding step (b) of the relax-round-and-resolve strategy described above. In particular, we greedily round Z to the closest binary solutionẐ encoding for disjoint . As a result, any matrix U with supportẐ automatically satisfies the orthogonality constraint. We can then obtain a solution of (13) by solving for each PC independently. For each t ∈ [r], we consider the submatrix of Σ over the indices {i :Ẑ i,t = 1}, extract its leading eigenvector via SVD, and pad it with zeros to construct U t . Although the restriction to disjoint support is not without loss of optimality, we will observe numerically in Section 4 that disjoint solutions are not particularly suboptimal for Problem (3) when k and r are small relative to p-an observation already made by Asteris et al. (2015).
We formalize our approach in Algorithm 1. In the rounding step, we require that each PC is supported on at least one feature, since otherwise the constraint U U = I cannot be satisfied.
Interestingly, Algorithm 2 of  can be seen as a special case of Algorithm 1, where there is only one PC and thus the solution is trivially disjoint.

Exact Non-Convex Formulation with Warmstart and Presolving
Problem (3) is a non-convex mixed-integer quadratically constrained problem which is addressable via non-convex MIQCP solvers. Since non-convex quadratic constraints cannot be satisfied exactly, we relax the orthogonality constraint U U = I to require that it is satisfied to within an elementwise s.t.
where δ i,j = 1 if i = j, and 0 otherwise. Problem (14) is a non-convex quadratically constrained mixed-integer problem with pr continuous variables, pr binaries, and r 2 quadratic constraints.
However, as previously discussed, spatial branch-and-bound technology can only solve problems with pr < 50. Therefore, we now propose techniques to accelerate the convergence of spatial branchand-bound by leveraging information from the convex relaxations derived in Section 2.
First, from the solution of a convex relaxation, Algorithm 1 generates a feasible solution (Z, U ) to Problem (14). We propose to use this solution as a high-quality warm-start for branch-andbound. Indeed, it is well documented in the mixed-integer optimization literature that warm-starts can dramatically accelerate the convergence of branch-and-bound (see, e.g., Bertsimas et al. 2016).
Second, when searching for a high-quality solution, we use information from the convex relaxation to pre-solve Problem (14) and restrict the search space. Namely, for any index i for which t∈[r] Z i,t < 10 −4 in the solution of the convex relaxation Z , we set t∈[r] Z i,t = 0 in (14). Recent work (c.f. Fischetti and Monaci 2013, Bertsimas and Digalakis 2022) similarly and successfully used information from convex relaxations to reduce the number of binary variables in other contexts. As we observe in our numerical results in Section 4, this numerical strategy accelerates branch-and-bound's convergence at the price of almost no loss in solution quality.
In addition, note that we strengthen Problem (14)'s formulation with valid inequalities derived from sparse PCA's 1 relaxation, as explored by Dey et al. (2022b,a) and in Section EC.2. Indeed, if each PC's sparsity k t is specified a priori, we have the valid inequalities Moreover, if k is specified but k t is not, we instead impose the second-order cone inequalities which allows us to model k t = p j=1 Z j,t in a more expensive yet tractable fashion. In practice, these warm-starting and presolving strategies obtained via Section 2's relaxations allow MIQCP solvers to obtain high-quality solutions in a reasonable amount of time for pr ≤ 50 (Section 4.2).
As spatial branch-and-bound technology improves over time, we believe that it should be possible to solve Problem (14) exactly at larger problem sizes. Indeed, recent works like Dong and Luo (2018), Gupta et al. (2022) solve some quadratically constrained problems with up to 50 variables to optimality using custom branch-and-bound solvers, and Gupta et al. (2022)

further reports that
Gurobi's off-the-shelf QCQP solver has achieved a machine independent speedup factor of 67.5 in less than two years, which suggests that larger instances of (14) (or (13)) may soon be within reach.
Unfortunately, we observe in our experiments in Section 4.1 that, due to numerical instability, the upper bounds obtained are currently not reliable for high values of (r, k) even when p = 13. Therefore, we use (14) primarily to obtain high-quality feasible solutions in this paper.

An Iterative Deflation Heuristic
In this section, we propose a local improvement technique which identifies near-optimal and nearfeasible solutions to Problem (3). The technique is based on the theory of Lagrangean relaxations (see Geoffrion 1974), which argues that if a non-convex problem is decomposable as a sum of easier (but still non-convex) subproblems with a coupling constraint, a good strategy is often to penalize the coupling constraint in the objective and iteratively solve the non-convex subproblems with different penalty multipliers on the coupling constraint. For example, Lu and Zhang (2012) propose an augmented Lagrangian method for solving an 1 relaxation of the sparse PCA, where the sparsity constraints are replaced by an 1 penalty in the objective.
To apply this perspective to our sparse PCA problem with multiple PCs, let us first rewrite (7) without the variable Y = t∈[r] Y t as: Hence, Problem (17) is the sum of r rank-1 sparse PCA problems coupled via linear constraints on their respective supports Z t and the semidefinite constraint t∈[r] Y t I. Therefore, if the coupling constraint can be handled appropriately then this problem can be addressed via scalable methods for rank-one sparse PCA, as reviewed in the introduction. To this end, we assume that a sparsity budget k t is imposed for each component as in Section 2.3, i.e., i∈[p] Z i,t ≤ k t , ∀t ∈ [r], and invoke Proposition 1 to replace the orthogonality constraint with r(r − 1)/2 bilinear scalar constraints Y t , Y t = 0, ∀t = t. Finally, since X, W ≥ 0 for any positive semidefinite matrices X, W of the same size, we replace Y t , Y t = 0, ∀t = t with Y t , Y t ≤ 0, ∀t = t without loss of generality. Therefore, for any penalty λ > 0 we have the following valid Lagrangean relaxation: Given an index t and a sparsity budget k t , optimizing for Y t (with all other Y t , t = t, and λ fixed) is equivalent to finding the leading k t -sparse PC of the matrix Σ − λ t =t Y t . Since tr(Y t ) = 1, we can add a term λ offset tr(Y t ) to the objective without impacting the optimal solution. Accordingly, in practice, we consider the matrix Σ−λ t =t Y t +λ offset I, where the constant λ offset > 0 is sufficiently large that the entire cost matrix is positive semidefinite, as required by most sparse PCA algorithms for r = 1. In our implementation, we solve these subproblems via the relax-and-round strategy from Section 3.1 in the special case r = 1 (analogous to ). Algorithm 2 proceeds by optimizing for each Y t sequentially and then iteratively increasing the penalty parameter λ > 0.
By increasing the penalty parameter λ, we can improve the orthogonality of the resulting PCs. Figure 1 illustrates the convergence of Algorithm 2 with the update rule λ = on the pitprops dataset with r = 3 and k = (10, 10, 10) as explored in in Section 4. Namely, Algorithm 2 initially identifies a set of high-quality but not orthogonal PCs, and converts this into a feasible solution by iteratively increasing the cost of violating the orthogonality constraint. In Section 4, for correlation matrices Σ, we use the update rule λ = for ≤ 50 and λ ← 50 + 5( − 50) for > 50, since, as Algorithm 2 Lagrangean Alternating Minimization for Problem (3) Require: Matrix Σ, rank parameter r, sparsity parameters k 1 , . . . , k r , number of iterations L Figure 1, using λ = throughout the method sometimes results in slow convergence.
For covariance matrices, we scale λ by a factor U U ,Σ p U U −I 1 , where U is a solution obtained after one iteration with λ = 0, to account for the different magnitude of the objective values. For futher scalability improvements, we also consider running Algorithm 2 only on indices i where t∈[r] Z i,t > in Z , a solution to the relaxation solved by Algorithm 1, as discussed in the previous section. We observe in our numerical results in Section 4 that this numerical strategy dramatically accelerates Algorithm 2 convergence at the price of a small loss in solution quality.

Numerical Results
In this section, we evaluate the algorithmic strategies derived in the previous two sections, imple- For a correlation matrix, this corresponds to dividing by p, the number of features.

Benchmarking on Pitprops Data
We first compare the bounds generated by each method-in terms of proportion of correlation explained-on the UCI pitprops dataset (p = 13) as we vary r ∈ {2, 3} and k ∈ {2r, 3r, 4r, 5r}, in Table 1. We consider both imposing an overall sparsity budget alone (denoted by "k, −") and imposing a seperate budget for each PC, denoted by "k, (k 1 , . . . , k t )".
We impose a relative optimality tolerance of 10 −6 for the exact formulation, set the parameters FuncPieceError and FuncPieceLength to 10 −6 and 10 −5 respectively (their minimum possible values), maintain an elementwise feasibility tolerance of 10 −8 for the orthogonality constraint, and report the upper bound at the time limit if branch-and-bound does not terminate.
On the instances presented in Table 1, we observe that spatial branch-and-bound on our nonconvex formulation (14) terminates within minutes and provides the tightest (i.e., smallest) upper bound. In comparison, the conic relaxations terminate in less than a second, while providing upper bounds that are only weaker at the third decimal. Unfortunately however, for r = 3 and r > 9, we started experiencing instability issues with spatial branch-and-bound. For example, for r = 3 and k = 12, it returned a global upper bound of 0  Table 1 Performance of upper bounds on the pitprops dataset (p = 13), as we vary the overall sparsity (k), the number of PCs (r) and the allocation of a sparsity budget to the different PCs. We denote the best performing solution (least upper bound) in bold. Note that all results are normalized by dividing by the trace of Σ, i.e., p, the number of features, to report results in terms of the proportion of variance explained.
by solving (14) may not be reliable as k or r increases or for p > 13. Accordingly, in the rest of this paper, we only use our convex relaxations to certify a solution's quality.
Regarding the two semidefinite relaxations "Disjoint-Ineq" and "Perm-Ineq", we should remind the reader that Perm-Ineq cannot compute a bound if we specify an overall sparsity budget k only.
Disjoint-Ineq's bound for a given k is much weaker than the worst-case bound over all possible allocations of k t 's that sum to k. This can be explained by the strong second-order cone inequalities (10)-(11) that can be added in the latter case. Therefore, time permitting, we recommend computing the lower bound by solving the relaxations for all possible allocations of k t and taking the worst-case bound. Furthermore, when individual sparsity budgets k t are given, we observe that Perm-Ineq provides uniformly and sometimes significantly tighter bounds than Disjoint-Ineq. Ac-cordingly, in the rest of the paper, we only consider instances of Problem (3) where we know both k and k t and consider Perm-Ineq and its second-order cone relaxations, but not Disjoint-Ineq.
For concision, we will mainly report the results for symmetric sparsity allocations, i.e., k t = k/r, only. In Table 1, we observe that for each value of (r, k), the worst-case upper bound is often achieved with asymmetric sparsity allocation, which suggests that sparse PCs with asymmetric sparsity allocation could explain more variance under a fixed sparsity budget. We investigate this matter in detail in Section 4.4.
Benchmarking on Larger-Scale Datasets We now investigate the scalability of the relaxation Perm-Ineq and its second-order cone relaxations on larger UCI datasets.  (12) with (10) We observe that the PSD relaxation can be solved within a few minutes for p ≈ 50 but quickly requires a prohibitive amount of memory and time at higher dimensions. The SOC-Cuts relaxation scales up to p ≈ 300 and provides a high-quality upper bounds, within 1 − 2% of PSD. Without the additional cuts, the SOC relaxation alone is weak when (k, r) are large relative to p. For example, for the pitprops dataset, when k t = (10, 10, 10), SOC returns an upper bound of 1.007 while the proportion of correlation explain can trivially not exceed 1. However, SOC is still preferable in high-dimensional settings, where the semidefinite formulation cannot be solved via an interior point method due to excessive memory requirements, and the formulation SOC-Cuts provides a small improvement in the upper bound at the price of significantly more runtime.

Performance of Feasible Methods
In this section, we numerically evaluate the quality of the three methods developed in Section 3 in terms of their ability to recover approximately orthogonal high-quality principal components on real-world datasets. We first validate that the variable fixing strategy proposed in Section 3.2 effectively reduces the number of variable without damaging the quality of the obtained solutions.
We then compare our algorithms with state-of-the-art techniques on seven UCI datasets. Since Algorithm 2 and all three benchmarked algorithms from the literature require that the sparsity of each PC is specified separately, we consider this formulation in this section (and fix k t = k/r).
Each of the algorithms proposed in Section 3 involves solving one of Section 2's convex relaxations. Based on the scalability results presented in the previous section, we use different convex relaxations depending on the dimensionality of the problem. Namely, for Algorithm 1, we use the full semidefinite relaxation (12) if p ≤ 50, the SOC relaxation (EC.5) with 50 PSD cuts if p ≤ 200, and the SOC relaxation (EC.5) otherwise. The exact method from Section 3.2 uses the solution from Algorithm 1 as a warm-start. We impose a two hour time limit for branch-and-bound and report the best solution found at the time limit. For Algorithm 2, we use the full semidefinite relaxation (12) when rounding each subproblem if p ≤ 50, and the SOC relaxation (EC.5) otherwise.
We also impose a limit of 100 iterations for p ≤ 200 or 200 iterations for p > 200.

Effectiveness of Variable Fixing:
We first apply our three methods on four small and mediumsized UCI datasets (Pitprops, Wine, Ionosphere, and Geography) as we vary r ∈ {2, 3} and k ∈ {5, 10, 20} : k ≤ p. Summary results are reported in Table 2 (see also  Table 2 Average performance of methods across the UCI pitprops, wine, ionosphere and geographical datasets with k ∈ {5, 10, 20}, r ∈ {2, 3} : k ≤ p. Note that the objective value is 1 p U U , Σ , the constraint violation is U U − I 1, and the objective value minus the constraint violation is 1 p U U , Σ − U U − I 1. For Algorithm 2, the variable fixing heuristic reduces computational time by nearly 10% and surprisingly improves the average objective from 0.391 to 0.410, thus halving the average relative optimality gap. By inspecting the instance-wise results in Table EC.3-EC.4, we observe that this gain happens mostly for the Geography dataset with p = 68 features. Recalling that we use a PSD relaxation for p ≤ 50 and a SOC relaxation otherwise, it is clear that the variable fixing heuristic allows the effective dimension of the Geography dataset to fall below 50, enabling "Algorithm 2 (fixed)" to use PSD relaxations, which are tighter than the SOC relaxation used by "Algorithm 2" for the same data. For the branch-and-bound approach, the scalability gains translate into a higher number of nodes explored (+25%) within the same time limit, although most instances are still not solved to optimality.
Comparing against the semidefinite bounds computed in the previous section, we observe that all methods return solutions which are at most 15% suboptimal on average and that Algorithm 2 (fixed) provides the best average performance (≤ 5% optimality gap on average).
Benchmarking on Pitprops Data: Lu and Zhang (2012) extensively compared the performance of existing sparse PCA methods on the pitprops dataset with r = 6 PCs. All six methods they benchmarked required an overall sparsity of around 45-60 to explain less than 70% of the variance.
The best performing method could explain 69.55% of the variance with an overall sparsity of 46 (Lu and Zhang 2012, Table 11). They concluded that "there do not exist six highly sparse, nearly orthogonal and uncorrelated PCs while explaining most of variance". The three methods we developed in Section 3 strongly challenge their conclusions. As reported in Table EC.2, with r = 6 PCs and k t = 2, hence an overall sparsity of 12, solutions returned by any of our methods explain 73%-75% of the variance. Algorithm 2 even provides a solution that explains 81% of the variance with an overall sparsity of 6 × 4 = 24. That is to say, what was previously considered by the community to be impossible can be done with the techniques in this paper, in seconds.

Benchmarking on Larger-Scale Datasets
We now investigate the performance of Algorithm 1, Algorithm 2 (fixed) and branch-and-bound (fixed) on seven UCI datasets, whose dimension range from p = 13 (pitprops) to p = 1300 (micromass). We compare them with three state-of-the-art methods from the literature. Namely, • The branch-and-bound method of Berk and Bertsimas (2019) for optimally computing one sparse PC, combined with the deflation scheme of Mackey (2008) to obtain multiple PCs, implemented in Julia and made available at github.com/lauren897/Optimal-SPCA. According to Berk and Bertsimas (2019), this method outperformed four others across three UCI datasets (r = 3, k = 5).
• The deflation method of Hein and Bühler (2010), using the custom deflation method developed in Bühler (2014), implemented in Matlab and made publicly available at github.com/tbuehler/sparsePCA, using default parameters. This approach was found by Berk and Bertsimas (2019 , Table 9) to be second-best of the methods in their comparison.
• The Lasso-inspired method of Zou et al. (2006), using the spca function in the elasticnet package version 1.3, using default parameters. This approach is perhaps the most commonly used one in practice, since it is distributed via the ubiquitous elasticnet package.
We report summary results in Table 3 (see also Tables EC.4 In terms of the average proportion of correlation explained minus the average orthogonality constraint violation, we observe that Algorithm 2 (fixed) and branch-and-bound (fixed) perform the best overall, with scores of 0.236 and 0.241 respectively, which compares favorably with the average upper bound computed by Algorithm 1 using the conic relaxations developed in Section 2 (0.276). Of the remaining methods, Algorithm 1 and the method of Berk and Bertsimas (2019) perform the next best, with respective scores of 0.232 and 0.221. Finally, the methods of Hein and Bühler (2010) and Zou et al. (2006) have the lowest respective scores, 0.200 and −1.530. One could also score the methods using the relative optimality gap alone; we avoid doing so here since some methods, e.g. the method of Berk and Bertsimas (2019) or Zou et al. (2006), often violate the orthogonality constraint significantly, thus potentially achieving artificially low optimality gaps.
We remark that no one method performs best on every instance. The method of Berk and Bertsimas (2019) performs best on instances where k is very small and where disjoint solutions are nearly optimal, in which case the orthogonality constraint can essentially be ignored (avg. relative gap of 3.38% and avg. constraint violation of 0.00875 over instances where p ≥ 101). Algorithm 2 (fixed), on the other hands, performs best on instances where k is large relative to p and the orthogonality constraint is important to account for (avg. relative gap of 2.55% and avg. constraint violation of 0.00614 over instances where p ≤ 34). These results suggest that both k and the amount of overlap between the optimal PCs impact the performance of each method, and motivate a comparison on synthetic data, where we control the ground truth, in the next section.
Finally, we note that while branch-and-bound (fixed) performs well on instances where p ≤ 34, Table EC.4 reveals that, for larger instances it does not improve upon the warm-start provided by Algorithm 1 within 7200s. Therefore, we do not consider the use of branch-and-bound any further in this paper. This situation should be revisited as global non-convex solvers improve.

Statistical Recovery Properties of Method
To evaluate the support recovery ability of each method, we now generate synthetic data according to a spiked covariance matrix model and compare the performance of Algorithms 1 and 2 (fixed) against the same three methods from the literature. We use an experimental setup inspired by d' Aspremont et al. (2008). Namely, we generate a test matrix Σ ∈ S p + , where Σ = U U + σx 1 x 1 + σx 2 x 2 , U ∈ [0, 1] p×p is a noisy matrix with i.i.d. uniform entries, σ = 2 is the signal-to-noise ratio, p = 50, and x 1 , x 2 are k true -sparse binary vectors (we fix k true = 20). We also control the proportion of overlap between the supports of x 1 and x 2 , q := x 1 x 2 0 x 1 0 , where denotes elementwise vector multiplication. Note that q = 0 and q = 1 correspond to the special cases of sparse PCA with disjoint support and row-sparsity discussed in Section 1.2. We impose a limit of 200 iterations for Algorithm 2. Figure 2 depicts the ROC curves (true positive rate vs. false positive rate for recovering the support of x 1 and x 2 ) as we vary k, averaged over 20 random instances, for low and high overlap (q ∈ { 1 20 , 19 20 }). Note that, since Algorithm 1 cannot allow any overlap between x 1 and x 2 , it converges (as k → p) towards setting Z i,t = 1 ∀i, t in the relaxation and randomly setting each Z i,t = 0 with probability 1/2 in the rounding step, leading the ROC curve to end at (1/2, 1/2). First, we observe high that overlap proportion negatively impacts the support recovery ability of all methods, since all ROC curves are noticeably closer to the diagonal line (random classifier) when the overlap proportion q increases. For low overlap (q = 1 20 ), Algorithm 2 appears to strictly dominate all other methods. We also observe that the methods of Hein and Bühler (2010), Berk and Bertsimas (2019) and Algorithm 1 behave almost identically for k ≤ 10, which aligns with our finding in Section 4.2 that these methods are comparable when k is small.
We confirm these findings with Table 4 that reports the Areas Under the ROC Curve (AUC) and the true positive and false positive rates (TPR, FPR) at the operator point k = k true = 20, for each method. Note that there is no clear correspondence between U 1 , U 2 and x 1 , x 2 , since any permutation of the columns of U gives an equivalent solution. Therefore, we compute the better true/false positive rates over both permutations. Moreover, to compute the AUC for Algorithm 1, we augment its ROC curve with the point (1, 1).  Table 4 Area under the Receiver Operator Curve (AUC) depicted in Figure 2 and True Positive Rate (TPR), False Positive Rate (FPR) at the operator point k = ktrue = 20 for each method, as we vary q. For each q, we denote the best performing method in terms of TPR/FPR in bold.
Across the three degrees of overlap, we observe that the method of Zou et al. (2006) is the least accurate, with a substantially larger false detection rate than the remaining methods. hence is outperformed by (resp. outperforms) Hein and Bühler (2010) for q = 1 20 (resp. q = 1 2 ). In terms of (TPR, FPR), Algorithm 1 actually dominates both Berk and Bertsimas (2019) and Hein and Bühler (2010), for q = 1 20 , 1 2 . For q = 19 20 (high overlap), Algorithm 1 and the method of Berk and Bertsimas (2019) obtain the best overall AUC. Algorithms 1 and 2 dominate the other methods in terms of (TPR, FDR).
However, as already observed in Figure 2, all AUCs are in the 0.5 − 0.6 range and the differences between methods are less acute than in the low-to medium-overlap cases.
Overall, Algorithm 2 dominates current state-of-the-art methods, in terms of (FDR, TPR) across all overlap regimes, and in terms of AUC in two out of three cases. Although it only returns PCs with disjoint support, Algorithm 1 performs is one of the best performing methods, achieving similar AUC as Berk and Bertsimas (2019).

Specifying the Sparsity Pattern: The Benefits of Asymmetry
While Problem (3) only requires a bound on the total sparsity, thus allowing flexibility on how this budget is allocated across PCs, the worst-case semidefinite upper bound over all sparsity patterns (k 1 , . . . , k r ) : t∈[r] k t = k is often significantly tighter than Problem (3)'s semidefinite relaxation with a sparsity budget of k alone, as demonstrated in Section 4.1. Moreover, Algorithm 2, which as demonstrated in Sections 4.2-4.3 is currently the best performing method for obtaining feasible solutions to Problem (3), requires that (k 1 , . . . , k t ) are individually specified. Collectively, these observations suggest that it may be necessary to enumerate all allocations of the sparsity budget k, which could be expensive. In our experiments, as is often done in practice, we restricted our search to symmetric allocations. In this section, we revisit the symmetry assumption, investigate when it is justified, and study the relative benefits of asymmetric sparsity budget allocations in terms of obtaining equally sparse sets of PCs that explain more variance.
We consider the pitprops, ionosphere, geographical and communities UCI datasets with a fixed number of PCs r = 3 and a given overall sparsity budget k ∈ {15, 30}. Accordingly, in Figures 3 and EC.1 of Section EC.5.5, we depict the relationship between the proportion of correlation explained in the data for each possible allocation of the sparsity budget (k 1 , k 2 , k 3 ) : k 1 + k 2 + k 3 = k, p ≥ k 1 ≥ k 2 ≥ k 3 ≥ 1, as computed by Algorithm 2 with a limit of 200 iterations and the same setup as in Section 4.2 (and the corresponding upper bound computed by Algorithm 1), against the relative asymmetry in the sparsity budget, as measured by KL((k 1 , k 2 , k 3 )||(k/3, k/3, k/3)) max p≥k 1 ≥k 2 ≥k 3 :k 1 +k 2 +k 3 =k KL((k 1 , k 2 , k 3 )||(k/3, k/3, k/3)) , where KL(p, q) := i p i log(p i /q i ) denotes the KL divergence. We observe a general trend that more symmetric sparsity budget allocations tend to explain more of the correlation in the data (both in terms of actual correlation explained by a solution from Algorithm 2 and in terms of the upper bound). This suggests that, when time is a concern, requiring that all PCs are equally sparse is a reasonable approach. Table 5  the proportion of correlation explained obtained in each case. We observe that in several instances a perfectly symmetric allocation of the sparsity budget yields the highest quality solution, and in all instances a perfectly symmetric allocation is within 10% in the worst-case (and within 1.65% in average) of the best solution. In the enumerated case, we also compute the optimality gap between the worst-case upper bound over all sparsity budget allocations, and the best solution found, and observe that on average it is less than 1% over the instances considered. Note that this is a different gap to the one reported in Section 4.2, where the upper bound is computed after assuming that all PCs are equally sparse.  Table 5 Comparison of symmetric and enumerated solution computed by Algorithm 2 for a given sparsity budget k total and given number of PCs (r = 3). We report the largest upper bound over all sparsity budget allocations as our enumerated upper bound, and report the relative optimality gap between the best solution and the worst-case bound. On average, considering asymmetric sparsity budget allocations improves the proportion of correlation explained by 1.65%.
All in all, more symmetric allocations of the sparsity budget tend to perform better on average.
Therefore, for a given sparsity budget k, a reasonable strategy could be to (a) run Algorithm 2 to compute a perfectly symmetric allocation, (b) compute Problem (12)'s upper bound across all possible allocations, and (c) run Algorithm 2 only on the asymmetric allocations for which the upper bound from (b) allows for a significant potential improvement upon the symmetric solution.
As a robustness check, we also report the quality of the solution and bound obtained by Algorithm 1 where we specify k but not k t in Section EC.5.6. We observe that the best solution obtained by Algorithm 1 without specifying k t is significantly worse than the solution obtained by Algorithm 2 with k 1 = k 2 = k 3 , and moreover there is little correlation between the k t 's identified by Algorithm 1 and the best k t 's obtained by Algorithm 2, which validates the idea that running Algorithm 2 with a perfectly symmetric sparsity allocation is a good idea in practice.

Summary and Guidelines From Numerical Experiments
In summary, our main findings from our numerical experiments are as follows: • The sparsity constraint on U can either be imposed via an overall sparsity budget of U 0 ≤ k or with component-specific sparsity constraints U t 0 ≤ k t . As reflected in Section 4.1, constraining each column separately yields substantially tighter conic relaxations. Moreover, some of the best performing algorithms for obtaining feasible solutions, including our Algorithm 2, explicitly require a constraint on the non-zero entries in each column. Therefore, constraining the sparsity of each column separately should be preferred wherever feasible.
• If practitioners have an overall sparsity budget but are agnostic about the sparsity of each column, a reasonable strategy is to require that all columns are equally sparse (i.e., set k t = k/r), as shown in Section 4.4. In our experiments, considering asymmetrically sparse sets of PCs increases the amount of variance explained by around 2% on average, at the price of increasing the total runtime by an order of magnitude.
• Under a separate sparsity constraint on each column of U , our proposed combination of solving a conic relaxation and either rounding to obtain a disjoint solution or running Algorithm 2 yields certifiably near optimal solutions within minutes (resp. hours) for problems with hundreds (resp. thousands) of features. Moreover, it substantially outperforms existing methods for sparse PCA with multiple PCs in terms of obtaining higher-quality solutions (Section 4.2) with a lower false discovery rate (Section 4.3). Therefore, it should be considered as a viable and more accurate alternative for sparse PCA problems with multiple PCs.

Conclusion
In this paper, we studied the problem of selecting a set of mutually orthogonal sparse principal components and proposed techniques which, for the first time, allow this problem to be solved to certifiable (near) optimality with 100s or 1000s of features in minutes or hours. In particular, we proposed a strong semidefinite relaxation (Section 2) which provides high-quality upper bounds on the amount of variance explainable by any set of sparse and mutually orthogonal components, and a suite of numerically efficient algorithms (Section 3) which, as demonstrated in Section 4, recover components that nearly match these upper bounds. This contributes towards an ever-growing body of work demonstrating that computationally challenging non-convex optimization problems can often be solved to provable (near) optimality in practice.  Moreover, since Y i,j = r t=1 Y t i,j is a rank-one quadratic under logical constraints Y t i,j = 0 if Z i,t = 0, invoking Lemma 3 reveals that this quadratic constraint's convex closure under these logical constraints is given by the strengthened inequality:

References
Second, in any feasible solution we have: Let us denote by k t the sparsity of the tth column of U , U t . Then, it is well known that U t 1 ≤ √ k t . Therefore: Next, squaring both sides and invoking the Cauchy-Schwarz inequality reveals that Finally, noting that the expression where δ i,j = 1{i = j} is an indicator denoting whether i = j. Summing these constraints over all indices i = j and using k − r + 1 as an upper bound on i∈ Finally, we recognize the right hand side as a sum of rank-one quadratic terms ( r t=1 Y t i,j ) 2 under logical constraints Y t i,j = 0 if Z j,t = 0 and invoke Lemma 3 to obtain the convex closure, giving:

EC.3. Proof of Proposition 2
Proof of Proposition 2 First, let us observe that if i∈[p] Z i,t ≤ k and Y t = U t U t is a rank-one matrix such that U 2 = 1 then we have U 1 ≤ √ k t by norm equivalence. Therefore Squaring both sides of this inequality then yields and combining Lemma 2 with this inequality yields (10).
Second, in the same spirit, since U t U t is only supported on indices where Z t is non-zero, we have that Y t Diag(Z t ). This constraint implies the following 2 × 2 minors are non-negative where δ i,j = 1 if i = j and 0 otherwise. Summing these inequalities over indices i = j and setting Finally, using Lemma 2 to take the convex closure of this inequality under the logical constraints Y t i,j = 0 if Z j,t = 0 gives Equation (11).

EC.4. A Second-Order Cone Relaxation
We now supply the full second-order cone relaxation we alluded to in Section 2:

EC.5. Supplementary Numerical Results
In this section, we provide supplementary results which support the numerical experiments performed in Section 4.   i.e., the number of features, to report in terms of the proportion of correlation explained. The notation "( * )" denotes that an instance could not be solved with 16 GB RAM and was instead solved using a high-performance computing environment with 256 GB RAM. The notation "−" denotes that neither environment could solve the instance. Next, we provide an instance-by-instance account of the results summarized in Table 2, in Table   EC.3. Note that we do not report on Algorithm 1, Algorithm 2 (fixed) or branch-and-bound (fixed) here, as these results are already reported in Table EC

EC.5.4. Instance-Wise Results on Larger UCI Datasets
Next, we provide an instance-by-instance account of the results summarized in Table 3, in Tables EC.4-EC.5.      Solution reported by Algorithm 1 for a given sparsity budget k total and number of PCs (r = 3), with kt not specified. We report the kt identified by Algorithm 1's rounding mechanism. The sparsity of the solutions returned by Algorithm 1 on the pitprops dataset can be explained by the fact that t kt ≤ p in Algorithm 1's solutions, since it returns disjoint solutions.