Iterative Algorithms for Sparse PCA

Why Use Iterative Algorithms?

Ordinary principal component analysis can find a leading direction by computing an eigenvector of the covariance matrix. Sparse PCA adds a constraint: the direction should have only a limited number of nonzero entries, so it can be interpreted using a small subset of variables. That constraint makes the optimization nonconvex, and a single eigenvector calculation no longer solves the problem in general. Iterative thresholding methods provide a practical way to seek a useful sparse direction without searching all possible supports.

Let \(S\) be a symmetric positive semidefinite \(p\times p\) covariance or sample covariance matrix. For a vector \(v\), write \(\|v\|_2\) for its Euclidean norm and \(\|v\|_0\) for the number of its nonzero coordinates. The cardinality-constrained sparse PCA objective is

where \(s\) is the chosen support-size limit. The value \(v^\top S v\) is the variance in direction \(v\). A global solution has the largest variance among all unit vectors with at most \(s\) nonzero coordinates. Iterative algorithms generally seek a good solution, not a certificate that this global maximum has been found.

Hard Thresholding and the Truncated Power Update

The update first applies the usual power-method multiplication \(v\mapsto Sv\), which amplifies directions associated with large eigenvalues. Hard thresholding then restores the desired sparsity before the vector is normalized. Normalization is essential because the optimization is over unit vectors; it does not change which coordinates thresholding keeps.

Hard thresholding has a precise projection interpretation, but that interpretation concerns Euclidean distance—not the sparse PCA objective itself.

To see this, fix a candidate support \(J\). Among vectors supported on \(J\), the squared distance \(\|z-u\|_2^2\) is minimized by setting \(u_j=z_j\) for \(j\in J\), giving residual \(\sum_{j\notin J}z_j^2\). Minimizing this residual over supports of size at most \(s\) means retaining the coordinates with the largest squared magnitudes. Thus hard thresholding is an exact Euclidean projection onto the sparse vectors, up to ties. It does not follow that it maximizes \(v^\top S v\).

Convergence: What Can Be Guaranteed?

Thresholded iterations are nonlinear because the selected support can change from one step to the next. Consequently, the familiar power-method convergence argument cannot be applied without qualification to the full algorithm. A clean local guarantee is available when the support has stabilized.

This is a conditional result: it describes convergence once the support stays fixed. It does not assert that the support will stabilize, that it will be the globally best support, or that the algorithm reaches the global sparse PCA solution.

Proof. Since \(v^{(T)}\) is supported on \(J\), and the update retains that support, multiplication and thresholding on subsequent steps act as multiplication by \(S_{JJ}\) followed by normalization. Thus, for \(k\geq 1\), the iterate \(v^{(T+k)}\), embedded in \(\mathbb R^p\), is \(S_{JJ}^k v^{(T)}\) normalized to unit length. Expand \(v^{(T)}\) in an orthonormal eigenbasis \(u_1,\ldots,u_{|J|}\): $$ v^{(T)}=a_1u_1+\sum_{j=2}^{|J|}a_ju_j, \qquad a_1=u_1^\top v^{(T)}\neq 0. $$ Then $$ S_{JJ}^kv^{(T)} =\lambda_1^k\left(a_1u_1+ \sum_{j=2}^{|J|}a_j\left(\frac{\lambda_j}{\lambda_1}\right)^ku_j\right). $$ For every \(j\geq2\), \(0\leq\lambda_j/\lambda_1<1\), so each term in the sum tends to zero as \(k\to\infty\). After normalization the vector therefore tends to \(\operatorname{sign}(a_1)u_1\). The sign is immaterial to the principal direction. \(\square\)

The theorem explains why power iteration can work well when the selected variables settle into a stable support and there is an eigenvalue gap on that support. When the support keeps changing, or the leading eigenvalue on a support is nearly tied with the next one, this argument gives no comparable guarantee. Multiple initializations are therefore useful in practice: they can reveal sensitivity to starting values, though they do not prove global optimality.

Initialization, Tuning, and Related Thresholding Methods

An initialization should be a unit vector, and its choice can affect which support the algorithm visits. Common choices include the leading eigenvector of \(S\), a vector formed from selected large-variance coordinates, and random unit vectors. The ordinary leading eigenvector uses all coordinates, so it can be a useful direction from which to start, but it may be poorly aligned with the sparse solution. Random starts can expose this sensitivity. In a careful analysis, run several starts and compare their final objective values and supports.

The sparsity parameter \(s\) controls the maximum support size, not necessarily the number of nonzero entries in every iterate or the true support size of a scientifically meaningful signal. A small \(s\) enforces interpretability but may omit useful variables; a large \(s\) reduces the restriction and can produce less parsimonious directions. Researchers can compare choices using stability across resamples, scientific constraints, or a separate model-selection procedure. The best choice depends on the goal and the data; the iteration itself does not select \(s\) automatically.

A related option is soft thresholding. For a threshold \(\tau\geq0\), define the coordinatewise operator $$ \bigl[\mathcal S_\tau(z)\bigr]_j =\operatorname{sign}(z_j)\max\{|z_j|-\tau,0\}. $$ A soft-thresholded power update replaces \(H_s(Sv)\) with \(\mathcal S_\tau(Sv)\), then normalizes. Larger \(\tau\) tends to remove more coordinates and shrink the surviving entries; unlike hard thresholding, it does not directly impose an exact support-size cap. If every coordinate is removed, the result is zero and normalization is undefined. An implementation must detect that case and, for example, change the threshold or restart.

Both hard and soft thresholding involve tuning choices: \(s\) for the former and \(\tau\) for the latter. Neither update is a general global optimizer for the cardinality-constrained variance objective. Their value is computational: each creates a sparse candidate direction using inexpensive coordinate operations.

Computational Cost and a Common Pitfall

For a dense \(p\times p\) matrix \(S\), multiplication by \(S\) costs \(O(p^2)\) operations per iteration. Finding the largest \(s\) absolute entries costs \(O(p\log p)\) with sorting, or can be done in expected \(O(p)\) time with selection algorithms. Thus the matrix multiplication usually dominates for a dense covariance matrix. If \(S=\mathbf X^\top\mathbf X/n\), where \(\mathbf X\) is a centered \(n\times p\) data matrix, an iteration can compute $$ Sv=\mathbf X^\top(\mathbf Xv)/n $$ without forming \(S\). This takes \(O(np)\) operations for dense data, in addition to the thresholding cost. Sparse data or fast matrix-vector products can reduce the multiplication cost. For \(K\) iterations and \(R\) starting values, the total cost is approximately \(K R\) times the per-iteration cost.

Common mistake: “Hard thresholding is the closest sparse approximation to \(Sv\), so one truncated power step must solve the sparse PCA optimization problem.” The first claim is correct in Euclidean distance; the conclusion is not. For example, take $$ S=\begin{pmatrix}2&1&0\\1&2&0\\0&0&2.5\end{pmatrix}, \qquad v=(0.8,0,0.6)^\top,\qquad s=1. $$ The vector is unit length, and \(Sv=(1.6,0.8,1.5)^\top\), so hard thresholding retains coordinate 1 and produces \((1,0,0)^\top\). Its sparse PCA objective is 2. But the feasible vector \((0,0,1)^\top\) has objective 2.5. The correction is to treat thresholding as a computational update, not as a global optimizer. Compare starts and assess the resulting directions rather than interpreting a single thresholding step as proof of optimality.

Iterative thresholding is most useful when a sparse direction is needed and repeated matrix-vector products are affordable. Its output depends on the sparsity tuning, initialization, and support changes. A stable support with a clear restricted eigenvalue gap offers a simple local convergence explanation; outside those conditions, convergence to a desired global solution should not be assumed.