Penalized Principal Component Analysis Using Smoothing

In order to allow for the addition of a penalty to the eigenvalue or eigenvector problem, we first formulate the computation of the leading eigenvector as an optimization problem (Golub and van Loan,, 1996). Precisely, letting $Q\in\mathbb{R}^{p\times p}$ be a symmetric and positive semidefinite matrix for a fixed dimension $p\in\mathbb{N}$, and $C\in\mathbb{R}^{p\times p}$ be a symmetric and strictly positive definite matrix, it can be shown that the solution of the optimization problem

is the leading eigenvector of the matrix $C^{-1}Q$. Following this observation, the penalized eigenvalue problem of Gaynanova et al., (2017) is defined by adding an $L_{1}$ penalty to the objective function being maximized in eq. (1), leading to the Lasso-type (Tibshirani,, 1996) objective function

where $C$ is assumed to be the identity matrix, and $\lambda>0$ is a penalty parameter. Note that in eq. (2), the objective is maximized, thus the penalty (which is non-negative) is being subtracted.

We employ the unified framework of Chen and Mangasarian, (1995); Nesterov, (2005); Trendafilov and Gallo, (2021) to smooth the objective function in eq. (2). We consider the smoothing of a piecewise affine and convex function $f:\mathbb{R}^{q}\rightarrow\mathbb{R}$, where $q\in\mathbb{N}$. As noted in Chen and Mangasarian, (1995); Nesterov, (2005); Trendafilov and Gallo, (2021), piecewise affine functions can be written as

where $z\in\mathbb{R}^{q}$ and $[z,1]\in\mathbb{R}^{q+1}$ denotes the concatenation of $z$ and the scalar $1$. The function in eq. (3) is called piecewise affine as it is composed of $k\in\mathbb{N}$ linear pieces, parameterized by the rows of $A\in\mathbb{R}^{k\times(q+1)}$. To be precise, each row in $A$ contains the coefficients of one of the $k$ linear pieces, with the constant coefficients being in column $q+1$.

A smooth surrogate for the function $f$ of eq. (3) is given by

In eq. (4), the optimization is carried out over the unit simplex in $k$ dimensions, given by

The function $\rho$ which appears in eq. (4) is called the prox-function. The prox-function must be nonnegative, continuously differentiable, and strongly convex. A specific choice of the prox-function is given below. The parameter $\mu>0$ is called the smoothing parameter. The choice $\mu=0$ recovers $f^{0}=f$.

The advantage of the smoothed surrogate of eq. (4) consists of the fact that it is both smooth for any choice $\mu>0$ and uniformly close to $f$, that is the approximation error is uniformly bounded as

see (Nesterov,, 2005, Theorem 1). As can be seen from eq. (4) and eq. (5), larger values of $\mu$ result in a stronger smoothing effect and a higher approximation error, while smaller values of $\mu$ result in a weaker smoothing effect and a higher degree of similarity between $f^{\mu}$ and $f$.

The choice of the prox-function $\rho$ is not unique, and in Chen and Mangasarian, (1995); Nesterov, (2005); Trendafilov and Gallo, (2021), several choices for the prox-function can be found. In the remainder of the article, we focus on one particular choice called the entropy prox-function as it allows for simple closed-form expressions of the smoothed objective for the penalized eigenvalue problem of eq. (2). The entropy prox-function is given by

and as shown in (Nesterov,, 2005, Section 4.1) and Hahn et al., 2021b ; Hahn et al., (2020), it satisfies the uniform bound

on the approximation error between $f^{\mu}$ and $f$.

Equipped with the results from Section 2.2, we now state the proposed smoothed version of the penalized eigenvalue problem of eq. (2). To this end, it is important to note that in eq. (2), the quadratic form $v^{\top}Qv$ is differentiable, while the $L_{1}$ penalty is not. We thus apply the smoothing to the $L_{1}$ penalty of the objective function only. Moreover, since $\|v\|_{1}$ in eq. (2) for $v\in\mathbb{R}^{p}$ can be written as $\|v\|_{1}=\sum_{i=1}^{p}|v_{i}|$, it suffices to smooth the non-differentiable absolute value applied to each entry of the vector $v$ separately.

The absolute value can be expressed in the form of eq. (3) with $k=2$ pieces, that is the function $f(z)=\max\{-z,z\}=\max_{i=1,2}\left(A[z,1]^{\top}\right)_{i}$ with

where $z\in\mathbb{R}$ is a scalar. We aim to replace the absolute value by the smooth surrogate with the entropy prox-function of eq. (6). Simplifying eq. (6) for the specific choice of $A$ results in

which is the smooth surrogate for the absolute value $|\cdot|$ we use. Leaving the quadratic form $v^{\top}Qv$ in eq. (2) unchanged, and substituting the $L_{1}$ penalty $\|v\|_{1}=\sum_{i=1}^{p}|v_{i}|$ by its smoothed version $\sum_{i=1}^{p}f_{e}^{\mu}(v_{i})$, then results in the smoothed penalized eigenvalue problem

Due to its simple form, the first derivative of $f_{e}^{\mu}$ can be stated explicitly as

Therefore, denoting the objective function in eq. (9) as

its gradient is given by

which is useful for the maximization of eq. (9) via gradient descent/ ascent solvers.

Solving the unsmoothed penalized eigenvalue problem of eq. (2), or the smoothed version of eq. (9), yields the leading eigenvector. Naturally, we are interested in further eigenvectors, in particular in Section 3 we will present population stratification plots of the first two leading eigenvectors, and compute a polygenic risk score using the first $10$ leading eigenvectors.

Deflation is a standard technique to extract higher order eigenvectors from a matrix $X\in\mathbb{R}^{m\times n}$, see Mackey, (2009); Trendafilov and Gallo, (2021). The idea is based on the Eckart-Young theorem (Eckart and Young,, 1936) for the SVD representation $X=U\Sigma V$, where $U$ is an $m\times m$ unitary matrix, $\Sigma$ is an $m\times n$ matrix with off-diagonal elements (that is, $\Sigma_{ij}$ for $i\neq j$) being zero, and $V\in\mathbb{R}^{n\times n}$ is a unitary matrix. The representation can be rewritten as $X=\sum_{i}\sigma_{i}u_{i}\otimes v_{i}$, where $\otimes$ denotes the outer product between two vectors and $\sigma_{i}$ is the $i$th entry on the diagonal of $\Sigma$. Using this representation allows one to subtract the leading eigenvector from $X$ once it is known.

We therefore first solve the unsmoothed penalized eigenvalue problem of eq. (2), or the smoothed version of eq. (9), and obtain the leading eigenvector of a matrix $X$, denoted as $v_{1}$. We compute the corresponding eigenvalue using the Rayleigh coefficient as $\alpha_{1}=v_{1}^{\top}Xv_{1}$. We then compute $X_{1}=X-\alpha_{1}v_{1}\otimes v_{1}$, the matrix with leading eigencomponent subtracted. As seen from the Eckart–Young theorem, the representation $X-\alpha_{1}v_{1}\otimes v_{1}=\sum_{i>1}\alpha_{i}v_{i}\otimes v_{i}$ still holds true, now with the second eigencomponent being the leading one. Applying the (unsmoothed or smoothed) penalized eigenvalue problem to $X-\alpha_{1}v_{1}\otimes v_{1}$ thus yields the second penalized eigenvector.

Using the smoothed version of the penalized eigenvalue problem of eq. (9) and its explicit gradient $\nabla F_{\lambda}$ given in Section 2.3, the minimization of eq. (9) is carried out using the quasi-Newton method BFGS (Broyden–Fletcher–Goldfarb–Shanno) which is implemented in the function optim in R (R Core Team,, 2014).

To improve the numerical accuracy of optimizing eq. (9), an iterative approach can be used. Given a target value of $\mu>0$ and an additional parameter $s\in\mathbb{N}$ (the number of steps), we start solving eq. (9) using a random initial starting value and the smoothing parameter $2^{s}\mu$. After the optimization is complete, we use the obtained maximizer as the initial value for a new run, this time with smoothing parameter $2^{s-1}\mu$. We continue seeding each new optimization with the solution of the previous run while lowering the smoothing parameter further until the target smoothing parameter $2^{0}\mu$ is reached. The target value of $\mu$ should be chosen as the machine precision or the square root of the machine precision. This is to ensure that the obtained estimates with smoothed PEP of eq. (9) closely resemble the ones obtained if the original PEP of eq. (2) had been solved.

Solving the smoothed PEP objective of eq. (9) yields an eigenvector which is not guaranteed to be sparse. This is due to the lack of an $L_{1}$ penalty in the formulation of eq. (9). Nevertheless, we can enforce sparsity again via thresholding. To be precise, let $v=v_{\lambda}^{\mu}$ be the vector obtained by maximizing eq. (9). Denoting the vector as $v=(v_{1},\dots,v_{p})$, we first choose a desired sparsity level $\pi\in[0,1]$. We then compute the threshold $\tau$ as the empirical $\pi$-quantile of the values $\{|v_{1}|,\ldots,|v_{p}|\}$, thus ensuring that a proportion $\pi$ of $v$ is smaller than $\tau$ in absolute value. Finally, we threshold all entries in $v$ against $\tau$, meaning we compute a new thresholded vector $\overline{v}$ with entries $\overline{v}_{i}=v_{i}\mathbb{I}(|v_{i}|>\tau)$ for all $i\in\{1,\ldots,p\}$.

In all of our experiments, whenever we refer to the smoothed version of PEP, we use the iterative solving approach outlined in Section 2.5 to solve eq. (9) with a target smoothing parameter of $\mu=0.1$. The smoothed version of PEP is used with $s=5$ steps, that is starting from $2^{5}\mu$. The choice of the parameter $\mu$ is investigated further in Section 3.6.

Moreover, in all of our experiments, sparsity is enforced via thresholding as described in Section 2.6. We use a target sparsity of $0.05$ (5% percent zeros).

All experiments are carried out using the statistics software $R$ (R Core Team,, 2014). The maximization of unsmoothed PEP in eq. (2) is conducted with the FISTA algorithm of Beck and Teboulle, (2009). The smoothed PEP of eq. (9) is carried out using the optim function in $R$. The parameter $\lambda$ of eq. (2) and eq. (9) is given individually in each simulation subsection.

In the comparison study of Section 3.9, we evaluate our proposed smoothed PEP to seven state-of-the-art sparse PCA algorithms. The details of those algorithms, including references to code we used, is given in Section 3.9.

The first dataset under consideration is the 1000 Genomes Project (The 1000 Genomes Project Consortium,, 2015, 2023). The goal of the 1000 Genomes Project was to discover genetic variants with frequencies of at least $1$ percent in the populations studied. In order to prepare the raw 1000 Genomes Project data and select rare variants, we used PLINK2 (Purcell and Chang,, 2019) with the cutoff value set to 0.01 for the –max-maf option. We then used the LD pruning method with the parameteters –indep-pairwise set to $2000$ $10$ $0.01$. These results focus on the European super population (subpopulations: CEU, FIN, GBR, IBS, and TSI), with our dataset containing 503 subjects and 5 million rare variants in total.

Next, we compute a similarity measure on the $503$ genomes included in the dataset. We employ the genomic relationship matrix (GRM) of Yang et al., (2011). We then compute two sets of leading eigenvectors, precisely the first two eigenvectors of the unsmoothed (see eq. (2)) and the smoothed (see eq. (9)) penalized eigenvector problems, respectively.

The second dataset under consideration consists of 1000 nucleotide sequences of the SARS-CoV-2 virus, which were extracted from patients diagnosed with Covid-19. The sequences were downloaded from the GISAID database (Elbe and Buckland-Merrett,, 2017; Shu and McCauley,, 2017). The sequences have accession numbers in the range of EPI_ISL_404227 to EPI_ISL_766048. All nucleotide sequences were aligned with MAFFT (Katoh et al.,, 2002) to the official SARS-CoV-2 reference sequence (available on GISAID under the accession number EPI_ISL_402124) using the keeplength option, and with all other parameters set to their default values. This procedure yields aligned sequences of length $29891$ base pairs.

After alignment, we compare each aligned nucleotide sequence against the reference sequence, resulting in a binary vector with $1$ indicating a mismatch to the reference genome, and $0$ otherwise. As before, we compute the genomic relationship matrix (GRM) of Yang et al., (2011) on the binary matrix obtained in this way (with each row corresponding to one sample), resulting in a $1000\times 1000$ similarity matrix for all subjects. The eigenvectors of the unsmoothed and smoothed PEP are computed on this similarity matrix (with $\lambda=0.1$ for smoothed PEP). Additional to the SARS-CoV-2 nucleotide sequence for each patient, we also obtain patient status meta information from GISAID, consisting of age (numeric), sex (male, female), geographic region (WHO geographic region encoded as WPRO, EURO, AFRO, PAHO, EMRO, SEARO), and binary mortality (patient deceased which is encoded with 1, or alive which is encoded with 0).

The third dataset under consideration is the Iris dataset of UC Irvine Machine Learning Repository (UC Irvine Machine Learning Repository,, 1988), a simple but widely used benchmark for clustering. The dataset contains $150$ observations (rows), where each observation is a plant. The true labels is the plant species (setosa, versicolor, virginica). The covariates are the plant’s sepal length, sepal width, petal length, and petal width (all measured as floating point numbers). The task is to recover the species from the four covariates measured for each plant.

Figure 1 displays the results of unsmoothed PEP in eq. (2) using the first two eigenvectors, colored by subpopulation (CEU, FIN, GBR, IBS, and TSI) of the 1000 Genomes Project dataset. We evaluated unsmoothed PEP for $\lambda$ equal to 0, 1, 10, and 100. We can see a good stratification for $\lambda$ equal to 0, 1, and 10, indicated by the relative discernibility of colored population clusters (notably becoming less discernable as $\lambda$ increases). We also observe that some entries of the principal component vectors are zero, leading to points which are located on the x-axis or y-axis. The amount of entries shrunk to zero seems to increase with $\lambda$, as expected. However, when $\lambda$ increases to 100, the populations are no longer discernible.

In Appendix A we present two additional figures showing population stratification plots for the dataset of the 1000 Genomes Project. Those two additional figures are computed with (a) the classic eigenvectors, as well as (b) the first two principal components as in eq. (2), however with an $L_{2}$ norm instead of the $L_{1}$ norm.

Similarly to Figure 1, Figure 2 shows the results for the smoothed version of PEP of eq. (9) using the first two eigenvectors, colored by subpopulation (CEU, FIN, GBR, IBS, and TSI). As before, we evaluated smoothed PEP for $\lambda$ equal to 0, 1, 10, and 100. However, in contrast to Figure 1, we see good stratification across all penalty values, again indicated by the relative discernibility of colored population clusters. While the clusters still become less discernable as $\lambda$ increases, for $\lambda=100$ we see a plot that is much more clearly clustered and conspicuous than the equivalent one in Figure 1. As before, the thresholding operation shrinks some of the entries in the first and second principal components to zero, leading to points which are located on the x-axis or y-axis. The amount of entries shrunk to zero seems to increase with $\lambda$, as expected.

The population stratification plots of Section 3.3 and Section 3.4 exhibit slightly different clusters by population. Naturally, our aim is to achieve a good separation between the clusters. In order to quantify the quality of the clustering achieved, we compute three metrics, the within- and between-cluster sums of squares as well as the silhouette score. As in Section 3.1, $n$ denotes the number of subjects, which is also the dimension of the matrix $Q\in\mathbb{R}^{n\times n}$ of Section 2. We denote the number of clusters with $k$, and the set of points contained in the $i$-th cluster with $C_{i}$, where $\cap_{i=1}^{k}C_{i}=\emptyset$ and $\cup_{i=1}^{k}C_{i}=\{1,\ldots,n\}$. Moreover, we denote each point in the population stratification plots with $v_{i}\in\mathbb{R}^{2}$ and the cluster this point belongs to with $c(i)\in\{1,\ldots,k\}$, where $i\in\{1,\ldots,n\}$. The three aforementioned metrics are defined as follows:

1. 1.

   The within-cluster sum of squares is defined as

   $\displaystyle\text{SS}_{\text{within}}=\sum_{i=1}^{n}\|v_{i}-\mu_{c(i)}\|_{2}^{2},$

   (10)

   where $\mu_{j}=\frac{1}{|C_{j}|}\sum_{r\in C_{j}}v_{r}$ denotes the mean of all the points in the $j$-th cluster $C_{j}$, where $j\in\{1,\ldots,k\}$, and $|C_{j}|$ denotes the number of elements in the set $C_{j}$.

2. 2.

   The between-cluster sum of squares is defined as

   $\displaystyle\text{SS}_{\text{between}}=\sum_{i=1}^{k}|C_{i}|\cdot\|\mu_{i}-\mu\|_{2}^{2},$

   (11)

   where $\mu=\frac{1}{n}\sum_{i=1}^{n}v_{i}$ is the mean of all the points $v_{i}$, where $i\in\{1,\ldots,n\}$.

3. 3.

   The silhouette score for point $v_{i}$ is defined as

   $\displaystyle\text{Silhouette}_{i}=\frac{b_{i}-a_{i}}{\max\{a_{i},b_{i}\}},$

   (12)

   where $a_{i}=\frac{1}{|C_{c(i)}|-1}\sum_{r\in C_{c(i)},r\neq i}\|v_{i}-v_{r}\|_{2}$ is the average distance of $v_{i}$ to all the other data points in the same cluster $C_{c(i)}$ as $v_{i}$, and the quantity $b_{i}=\min_{j\neq c(i)}\frac{1}{|C_{j}|}\sum_{r\in C_{j}}\|v_{i}-v_{r}\|_{2}$ is the minimum among all the mean distances of $v_{i}$ to the other clusters $j\neq c(i)$ that $v_{i}$ does not belong to. We compute the silhouette score for each point with the help of the function silhouette of the R-package cluster (Maechler et al.,, 2022). We report the mean over all silhouette scores as a summary metric.

By evaluating $\text{SS}_{\text{within}}$, $\text{SS}_{\text{between}}$, and $\text{Silhouette}_{i}$ for the projected data points, we assess the effectiveness of the penalized eigenvalue problem solutions in separating clusters in meaningful ways. These metrics provide a robust way to compare the results obtained from the unsmoothed PEP of Section 3.3 and the smoothed PEP of Section 3.4.

Tables 1, 2, and 3 show the within-cluster sum of squares, between-cluster sum of squares, and silhouette scores for the clusters displayed in Figure 1 and Figure 2 (again for $\lambda\in\{0,1,10,100\}$). For the within-cluster sum of squares and the silhouette score (Tables (1) and (3)), a smaller value indicates better clustering, and for the between-cluster sum of squares metric (Table (2)), a larger value indicates a better clustering. Across all tables, for both the smoothed and unsmoothed PEP, we can see that as we increase the Lasso penalty $\lambda$, the quality of the clusters decreases across the board. We also note that as the $\lambda$ values increase, the smoothed models achieve better clustering in comparison to the unsmoothed models across all three metrics. Therefore, we posit that this smoothing method enhances the performance of these models via more discernable clustering.

We investigate the effect of the smoothing parameter $\mu$ on the obtained clustering. For this, we compute the figures of Section 3.4 with different choices of $\mu$. We note that if the Lasso penalty $\lambda$ in eq. (2) is negligible, then so will be the $L_{1}$ norm penalty for the objective function in eq. (2). Therefore, exchanging the $L_{1}$ norm for the smoothed counterpart will likely have a minor effect only. Therefore, we choose a relatively large Lasso penalty of $\lambda=1$ in order to see the influence of the $L_{1}$ penalty and, accordingly, of its smoothed approximation of Section 2.3.

Figure 3 shows results for the fixed choice $\lambda=1$ and varying values of $\mu\in\{0.01,0.1,1,10\}$. In contrast to Figure 1 and Figure 2, the points in Figure 3 are not colored by population label in the 1000 Genomes Project data but unison by value of $\mu$, precisely black for $\mu=0.01$, blue for $\mu=0.1$, green for $\mu=1$, and red for $\mu=10$. We observe that as $\mu$ increases, the plots become more and more compressed in one spot, whereas smaller values of $\mu$ nicely stratify the data. We explain this finding with the fact that for a fixed $\lambda$, larger values of $\mu$ will excessively smooth out the objective function, thus rendering the maximum delocalized. In contrast, smaller values of $\mu$ cause the smoothed objective of eq. (9) to be pointwise close to the unsmoothed one of eq. (2), thus recovering the original principal components.

We aim to compare the impact of using either the unsmoothed or the smoothed PEP eigenvectors to compute a polygenic risk score to predict SARS-CoV-2 mortality. To be precise, using the dataset described in Section 3.1, we fit a linear regression model to the outcome (mortality as binary vector) as a function of age, sex, geographic region, and $10$ principal components computed on the GRM (genomic relationship matrix) similarity measure of the SARS-CoV-2 nucleotide sequences.

While keeping the regression model unchanged otherwise, we once use $10$ principal components computed with the unsmoothed PEP of eq. (2) and smoothed PEP of eq. (9). For smoothed PEP we employ $\lambda=0.1$. We train the model on a proportion of $\pi\in\{0.1,\ldots,0.9\}$ of the individuals, and evaluate the prediction against the true outcome on the withheld proportion of $1-\pi$ individuals. As the response is binary, we employ the AUC (area under the curve) metric to assess the accuracy of the prediction.

Figure 4 shows the results of this experiment as a function of the proportion $\pi$ used for training, both for the model using $10$ principal components from the unsmoothed (blue) and the smoothed (red) PEP. With an AUC above $0.8$, we observe that indeed, the covariates age, sex, geographic region, and the $10$ principal components of the nucleotide data are good predictors of mortality. Moreover, we observe that the prediction accuracy increases with an increasing training proportion, which is to be expected. Importantly, using the principal components from the smoothed PEP seems to give a slight edge in prediction in comparison to using the unsmoothed PEP eigenvectors.

As a second real data example, we aim to visually assess the clustering obtained on the Iris benchmark dataset of the UC Irvine Machine Learning Repository (UC Irvine Machine Learning Repository,, 1988). This simple but widely used benchmark for clustering contains data on $150$ plants, each belonging to one of three species (setosa, versicolor, virginica). The covariates are the plant’s sepal length, sepal width, petal length, and petal width.

As described in Section 3.1, we first compute a similarity measure on the dataset to obtain a similarity matrix of dimensions $150\times 150$, where each matrix entry is a measure of similarity between any pair of plants. As a similarity measure we use the Jaccard similarity matrix computed on $X$ (Jaccard,, 1901; Prokopenko et al.,, 2016; Hahn et al., 2021a, ), where $X\in\mathbb{R}^{150\times 4}$ is the Iris dataset (150 observations and 4 covariates). Afterwards, we compute two principal components with both unsmoothed and smoothed PEP (for $\lambda=0.075$) and plot them. The results are shown in Figure 5, where each point is colored based on its true label. We observe that the principal components computed with smoothed PEP lead to visibly clearer clustering by species.

We compare our smoothed PEP approach of eq. (9) with seven other state-of-the-art algorithms on simulated data. Those approaches are the power method for sparse principal component analysis of Journée et al., (2008), the semidefinite programming formulation of d’Aspremont et al., (2007), the penalized orthogonal iteration of Jung et al., (2019), the penalized formulation of Gaynanova et al., (2017), minorization-maximization algorithm for sparse PCA of Song et al., (2015), the first-order algorithm of Zhang et al., (2010), and the sPCA-rSVD algorithm Shen and Huang, (2008).

Specifically, we implement Algorithm 6 of Journée et al., (2008) (referred to as Journee6), Algorithm 2 of d’Aspremont et al., (2007) (referred to as Aspremont2), Algorithm 2 of Jung et al., (2019) (referred to as Jung2), Algorithm 1 of Gaynanova et al., (2017) (referred to as Gaynanova1), Algorithm 4 of Song et al., (2015) (referred to as Song4), Algorithm 1 of Zhang et al., (2010) (referred to as Zhang1), and Algorithm 1 of Shen and Huang, (2008) (referred to as Shen1). Our codes also use functionality of the R-packages RSpectra (Qiu et al.,, 2022), alabama (Varadhan,, 2023), irlba (Baglama et al.,, 2022), and epca (Chen,, 2023).

In order to be able to know the true sparse eigenvector, we employ the following simulation scenario. We generate square matrices $U$ of dimensions $n\in\{10,20,50,100\}$ with a planted first eigenvector in its first column. The sparsity $\rho$ (proportion of zeros) of the first eigenvector is varied in $\rho\in\{0.1,0.2,0.3,0.4,0.5\}$. The nonzero entries of the first eigenvector are drawn independently from a uniform distribution in $[0,1]$. The remaining columns are also randomly drawn from a uniform distribution in $[0,1]$. We also generate a diagonal matrix $D$ of dimension $n$ with randomly drawn eigenvalues, sorted in decreasing order, on its diagonal. We use both $U$ and $D$ as the inputs of an eigendecomposition, and thus generate a test matrix as $X=UDU^{-1}$.

Using $X$ constructed in the aforementioned way, we apply our smoothed PEP algorithm of Section 2.5 (with $\lambda=0.1$) and the aforementioned seven state-of-the-art algorithms to recover the first sparse eigenvector of $X$. As the methods of Journée et al., (2008) and Zhang et al., (2010) return a full matrix, we only consider its first column. We assess all algorithms with respect to the accuracy of the obtained eigenvectors, their support recovery, and their runtime. This is done on the same data, that is we apply all algorithms and then return four metrics for each result. To be precise, we assess the accuracy of the returned sparse eigenvector by calculating the cosine similarity measure between the returned sparse eigenvector and the true (planted) eigenvector (ranging from $0$ to $1$, where $0$ encodes maximual disagreement of the two vectors, and $1$ encodes maximal agreement). We also consider the support recovery, meaning the ability to recover the true nonzero entries. We measure this with the true positive rate, meaning the proportion of truly recovered nonzero entries (ranging from $0$ to $1$, where $0$ is worst and $1$ is best). At the same time, we also report the false positive rate, meaning the proportion of falsely proclaimed nonzero entries (ranging from $0$ to $1$, where $0$ is best and $1$ is worst). Finally, we report the runtime of all algorithms in seconds.

The accuracy results are given in Table 4. We observe that, as expected, the accuracy of the solution decreases with $n$. Interestingly, we do not observe a strong dependence of the accuracy on the sparsity $\rho$. All algorithms achieve a reasonable cosine similarity with the planted true eigenvector. However, the algorithms of Gaynanova et al., (2017), Shen and Huang, (2008), and our smoothed PEP perform best.

Next, we look at the support recovery, specifically the ability to recover the true nonzero entries measured with the True Positive Rate in Table 5. We observe that all methods, possibly apart from the one of Shen and Huang, (2008) as $n$ increases, accurately recover the support in this experiment. It is noteworthy that the algorithms of Jung et al., (2019), Gaynanova et al., (2017), Zhang et al., (2010), and our smoothed PEP perform particularly well in this task. At the same time, we also assess the False Positive Rate, meaning the proportion of falsely identified nonzero entries in Table 6. Table 6 again shows that all methods, possibly apart from the ones of Song et al., (2015) and Shen and Huang, (2008), achieve a low false positive rate.

Finally, we assess the runtime (in seconds) of all algorithms in Table 7. We observe most methods are able to compute the leading sparse eigenvector in a fraction of a second for all matrices under consideration. Notably, the algorithms of Journée et al., (2008) and d’Aspremont et al., (2007) are much slower than the others.

The authors have no competing interests to declare that are relevant to the content of this article.

Funding for this research was provided through Cure Alzheimer’s Fund, the National Institutes of Health [1R01 AI 154470-01; 2U01 HG 008685; R21 HD 095228 008976; U01 HL 089856; U01 HL 089897; P01 HL 120839; P01 HL 132825; 2U01 HG 008685; R21 HD 095228, P01 HL132825], the National Science Foundation [NSF PHY 2033046; NSF GRFP 1745302], and a NIH Center grant [P30-ES002109].

The authors thankfully acknowledge support of an ERC training grant (T42 OH008416). The authors gratefully acknowledge the contributors of the 1000 Genomes Project (The 1000 Genomes Project Consortium,, 2023).

The authors gratefully acknowledge the contributors, originating and submitting laboratories of the sequences from GISAID’s EpiCoV™ Database (Elbe and Buckland-Merrett,, 2017; Shu and McCauley,, 2017) on which this research is based. The accession numbers of all nucleotide sequences used in this work are given in the supplementary material.

The datasets generated and/or analysed during the current study are available in the International Genome Sample Resource repository at the address https://www.internationalgenome.org/.

The data that support the findings of this study are publicly available in the GISAID database (Elbe and Buckland-Merrett,, 2017; Shu and McCauley,, 2017), see https://gisaid.org/. The supplementary material of this manuscript contains a list of all EPI_ISL identifiers of the nucleotide sequences which were included in the experiments.