Multi-Rank Sparse and Functional PCA: Manifold Optimization and Iterative Deflation Techniques

TL;DR

This paper introduces advanced methods for estimating multiple sparse and functional principal components using manifold optimization and iterative deflation, improving efficiency and performance over existing techniques.

Contribution

It extends SFPCA to multiple components with manifold optimization and develops efficient deflation algorithms for better estimation.

Findings

Manifold optimization effectively estimates multiple components.

Iterative deflation improves computational efficiency.

Both methods outperform traditional approaches in non-orthogonal scenarios.

Abstract

We consider the problem of estimating multiple principal components using the recently-proposed Sparse and Functional Principal Components Analysis (SFPCA) estimator. We first propose an extension of SFPCA which estimates several principal components simultaneously using manifold optimization techniques to enforce orthogonality constraints. While effective, this approach is computationally burdensome so we also consider iterative deflation approaches which take advantage of existing fast algorithms for rank-one SFPCA. We show that alternative deflation schemes can more efficiently extract signal from the data, in turn improving estimation of subsequent components. Finally, we compare the performance of our manifold optimization and deflation techniques in a scenario where orthogonality does not hold and find that they still lead to significantly improved performance.

Tables

Table 1. Table 1 : Properties of Hotelling’s Deflation (HD), Projection Deflation (PD), and Schur Complement Deflation (SD). Only SD captures all of the individual signal of each principal component without re-introducing signal at later iterations. Additionally, only SD allows for the non-unit-norm PCs estimated by SFPCA to be used without rescaling.

Method	Two-Way Orthogonality	One-Way Orthogonality	Subsequent Orthogonality ($\forall s\ge 0$)	Robust to
	${𝒖}_t^T{𝑿}_t{𝒗}_t=0$	${𝒖}_t^T{𝑿}_t,{𝑿}_t{𝒗}_t=\mathrm{𝟎}$	${𝒖}_t^T{𝑿}_{t+s},{𝑿}_{t+s}{𝒗}_t=\mathrm{𝟎}$	Scale of ${𝒖}_t,{𝒗}_t$
Hotelling’s Deflation (HD)	✓	✗	✗	✗
Projection Deflation (PD)	✓	✓	✗	✗
Schur Deflation (SD)	✓	✓	✓	✓

Table 2. Table 2 : Comparison of Manifold-ADMM, Manifold Proximal Gradient, and Alternating Manifold Proximal Gradient approaches for Manifold SFPCA ( 2 ). MADMM is consistently more efficient and obtains better solutions, but both MADMM and A-ManPG perform well in terms of signal recovery, as measured by relative subspace recovery error ( rSS-Error = ‖ 𝑼 ^ ​ 𝑼 ^ T − 𝑼 ∗ ​ ( 𝑼 ∗ ) T ‖ / ‖ 𝑼 ^ SVD ​ 𝑼 ^ SVD T − 𝑼 ∗ ​ ( 𝑼 ∗ ) T ‖ rSS-Error norm ^ 𝑼 superscript ^ 𝑼 𝑇 superscript 𝑼 superscript superscript 𝑼 𝑇 norm subscript ^ 𝑼 SVD superscript subscript ^ 𝑼 SVD 𝑇 superscript 𝑼 superscript superscript 𝑼 𝑇 \text{rSS-Error}=\|\hat{\bm{U}}\hat{\bm{U}}^{T}-\bm{U}^{*}(\bm{U}^{*})^{T}\|/\|\hat{\bm{U}}_{\text{SVD}}\hat{\bm{U}}_{\text{SVD}}^{T}-\bm{U}^{*}(\bm{U}^{*})^{T}\| ), true positive rate (TPR) and false positive rate (FPR).

		MADMM	ManPG	A-ManPG
Time (s)		19.46	217.40	175.99
Suboptimality		0	4.58	12.23
rSS-Error	$𝑼$	68.66%	74.03%	64.54%
	$𝑽$	36.85%	50.69%	43.17%
TPR	$𝑼$	87.75%	69.72%	89.99%
	$𝑽$	94.75%	70.30%	89.87%
FPR	$𝑼$	12.25%	30.28%	10.01%
	$𝑽$	5.25%	29.70%	10.13%

Table 3. Table 3 : Cumulative Proportion of Variance Explained (CPVE) of Rank-One SFPCA with (normalized) Hotelling, Projection, and Schur Complement Deflation and of (order 3) Manifold SFPCA. SD gives the best CVPE of the iterative approaches and appears to be more robust to violations of orthogonality (Scenario 2). Manifold SFPCA outperforms the iterative methods in both scenarios.

Scenario 1: ${𝐔}^{\ast }$ and ${𝐕}^{\ast }$ Orthogonal
		HD	PD	SD	ManSFPCA
CPVE	PC1	15.92%	21.05%	21.87%	37.12%
	PC2	22.21%	29.42%	30.59%
	PC3	26.80%	35.57%	37.09%
Scenario 2: ${𝐔}^{\ast }$ and ${𝐕}^{\ast }$ Not Orthogonal
CPVE	PC1	8.85%	19.74%	29.80%	50.85%
	PC2	13.03%	28.30%	39.87%
	PC3	16.16%	34.22%	46.48%

Table 4. Table 4 : Extended Version of Table 3 , comparing Cumulative Proportion of Variance Explained (CVPE), relative subspace recovery error ( rSS-Error = ‖ 𝑼 ^ ​ 𝑼 ^ T − 𝑼 ∗ ​ ( 𝑼 ∗ ) T ‖ / ‖ 𝑼 ^ SVD ​ 𝑼 ^ SVD T − 𝑼 ∗ ​ ( 𝑼 ∗ ) T ‖ rSS-Error norm ^ 𝑼 superscript ^ 𝑼 𝑇 superscript 𝑼 superscript superscript 𝑼 𝑇 norm subscript ^ 𝑼 SVD superscript subscript ^ 𝑼 SVD 𝑇 superscript 𝑼 superscript superscript 𝑼 𝑇 \text{rSS-Error}=\|\hat{\bm{U}}\hat{\bm{U}}^{T}-\bm{U}^{*}(\bm{U}^{*})^{T}\|/\|\hat{\bm{U}}_{\text{SVD}}\hat{\bm{U}}_{\text{SVD}}^{T}-\bm{U}^{*}(\bm{U}^{*})^{T}\| ), true positive rate (TPR) and false positive rate (FPR) of Manifold SFPCA with Rank-One SFPCA using (normalized) Hotelling, Projection and Schur Complement Deflation.

		Scenario 1: ${𝐔}^{\ast }$ and ${𝐕}^{\ast }$ Orthogonal				Scenario 2: ${𝐔}^{\ast }$ and ${𝐕}^{\ast }$ Not Orthogonal
		HD	PD	SD	ManSFPCA	HD	PD	SD	ManSFPCA
CPVE	PC1	15.92%	21.05%	21.87%	37.12%	8.85%	19.74%	29.80%	50.85%
	PC2	22.21%	29.42%	30.59%		13.03%	28.30%	39.87%
	PC3	26.80%	35.57%	37.09%		16.16%	34.22%	46.48%
rSS-Error	$𝑼$	129.54%	129.55%	128.35%	69.32%	215.73%	206.30%	205.74%	97.77%
	$𝑽$	143.01%	143.72%	141.15%	36.98%	211.15%	207.77%	204.38%	78.26%
TPR	$𝑼$	54.57%	54.29%	54.63%	87.72%	11.31%	13.11%	15.28%	88.91%
	$𝑽$	59.20%	58.84%	59.41%	94.72%	12.89%	15.33%	18.50%	80.08%
FPR	$𝑼$	0.92%	0.91%	0.92%	12.28%	1.24%	1.05%	1.18%	11.09%
	$𝑽$	0.61%	0.60%	0.60%	5.28%	2.80%	2.72%	2.79%	19.92%

Equations

$$\underset{\mathbf{u}\in {\overline{\mathbb{B}}}_{{\mathbf{S}}_{\mathbf{u}}}^n,\mathbf{v}\in {\overline{\mathbb{B}}}_{{\mathbf{S}}_{\mathbf{v}}}^p}{\mathrm{argmax}}{\mathbf{u}}^T\mathbf{X}\mathbf{v}-{\lambda }_{\mathbf{u}}{P}_{\mathbf{u}}(\mathbf{u})-{\lambda }_{\mathbf{v}}{P}_{\mathbf{v}}(\mathbf{v})$$

$$\underset{\mathbf{U}\in {\mathcal{V}}_{n\times k}^{{\mathbf{S}}_{\mathbf{u}}},\mathbf{V}\in {\mathcal{V}}_{p\times k}^{{\mathbf{S}}_{\mathbf{v}}}}{\mathrm{argmax}}\mathrm{Tr}({\mathbf{U}}^T\mathbf{X}\mathbf{V})-{\lambda }_{\mathbf{U}}{P}_{\mathbf{U}}(\mathbf{U})-{\lambda }_{\mathbf{V}}{P}_{\mathbf{V}}(\mathbf{V})$$

$$\mathbf{U}\in {\mathcal{V}}_{n\times k}^{{\mathbf{S}}_{\mathbf{u}}}⟺\mathbf{U}\in {\mathbb{R}}^{n\times k}\text{ and }{\mathbf{U}}^T{\mathbf{S}}_{\mathbf{u}}\mathbf{U}={\mathbf{I}}_k.$$

$$\hat{\mathbf{U}}=\underset{\mathbf{U}\in {\mathcal{V}}_{n\times k}^{{\mathbf{S}}_{\mathbf{u}}}}{\mathrm{argmin}}-\mathrm{Tr}({\mathbf{U}}^T\mathbf{X}\hat{\mathbf{V}})+{\lambda }_{\mathbf{U}}{P}_{\mathbf{U}}(\mathbf{U})$$

$$\hat{\mathbf{V}}=\underset{\mathbf{V}\in {\mathcal{V}}_{n\times p}^{{\mathbf{S}}_{\mathbf{v}}}}{\mathrm{argmin}}-\mathrm{Tr}({\hat{\mathbf{U}}}^T\mathbf{X}\mathbf{V})+{\lambda }_{\mathbf{V}}{P}_{\mathbf{V}}(\mathbf{V})$$

$$\hat{\mathbf{D}}=\underset{\begin{array}{c}\mathbf{D}\in {\mathbb{R}}^{n\times k}\end{array}}{\mathrm{argmin}}-⟨\mathbf{X}\hat{\mathbf{V}},\mathbf{D}{⟩}_F+{\lambda }_{\mathbf{U}}{P}_{\mathbf{U}}({\mathbf{U}}^{(k)}+\mathbf{D})$$

$${\mathbf{U}}^{(k+1)}$$

$${\mathbf{W}}^{(k+1)}$$

$$={\mathrm{prox}}_{{\lambda }_{\mathbf{U}}/\rho P_{\mathbf{U}}(\cdot )}\left({\mathbf{U}}^{(k+1)}+{\mathbf{Z}}^{(k)}\right)$$

$${\mathbf{Z}}^{(k+1)}$$

$${\mathbf{X}}_t:={\mathbf{X}}_{t-1}-d_t{\mathbf{u}}_t{\mathbf{v}}_t^T\text{ where }{d}_t={\mathbf{u}}_t^T{\mathbf{X}}_{t-1}{\mathbf{v}}_t.$$

$${\mathbf{X}}_t:=({\mathbf{I}}_n-{\mathbf{u}}_t{\mathbf{u}}_t^T){\mathbf{X}}_{t-1}({\mathbf{I}}_p-{\mathbf{v}}_t{\mathbf{v}}_t^T).$$

$${\mathbf{X}}_t:={\mathbf{X}}_{t-1}-\frac{{\mathbf{X}}_{t-1}{\mathbf{v}}_t{\mathbf{u}}_t^T{\mathbf{X}}_{t-1}}{{\mathbf{u}}_t^T{\mathbf{X}}_{t-1}{\mathbf{v}}_t}.$$

$${\mathbf{X}}_t^{\text{HD}}$$

$${\mathbf{X}}_t^{\text{PD}}$$

$${\mathbf{X}}_t^{\text{SD}}$$

$${\mathbf{u}}_t^T{\mathbf{X}}_t{\mathbf{v}}_t$$

$$={\mathbf{u}}_t^T{\mathbf{X}}_{t-1}{\mathbf{v}}_t-d_t{\mathbf{u}}^T{\mathbf{u}}_t{\mathbf{v}}_t^T{\mathbf{v}}_t$$

$$=d_t-d_t\Vert {\mathbf{u}}_t{\Vert }^2\Vert {\mathbf{v}}_t{\Vert }^2$$

$$\mathbf{X}=\left(\begin{array}{cc}2& -4/3\\ 2& 2/3\\ 1& 4/3\end{array}\right)=\underset{=\mathbf{U}}{\underbrace{\left(\begin{array}{ccc}2/3& -2/3& 1/3\\ 2/3& 1/3& -2/3\\ 1/3& 2/3& 2/3\end{array}\right)}}\underset{=\mathbf{D}}{\underbrace{\left(\begin{array}{ccc}3& 0& 0\\ 0& 2& 0\\ 0& 0& 1\end{array}\right)}}\underset{={\mathbf{V}}^T}{\underbrace{{\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right)}^T}}.$$

$${\mathbf{X}}_1=\left(\begin{array}{cc}0& -4/3\\ 0& 2/3\\ 1& 4/3\end{array}\right).$$

$${\mathbf{u}}_t^T{\mathbf{X}}_t$$

$$=({\mathbf{u}}_t^T-\Vert {\mathbf{u}}_t{\Vert }^2{\mathbf{u}}_t^T){\mathbf{X}}_{t-1}({\mathbf{I}}_p-{\mathbf{v}}_t{\mathbf{v}}_t^T)$$

$$\mathbf{X}=\left(\begin{array}{ccc}[r]-2& -3/2& 1\\ 8/3& 1/6& 1/3\\ 0& 5/2& 1\\ 2/3& 7/6& 7/3\end{array}\right)=\underset{=\mathbf{U}}{\underbrace{\left(\begin{array}{cccc}[r]-1/2& 1/2& 1/2& 1/2\\ 1/2& -1/2& 1/2& 1/2\\ 1/2& 1/2& -1/2& 1/2\\ 1/2& 1/2& 1/2& -1/2\end{array}\right)}}\underset{=\mathbf{D}}{\underbrace{\left(\begin{array}{cccc}[r]4& 0& 0& 0\\ 0& 3& 0& 0\\ 0& 0& 2& 0\\ 0& 0& 0& 1\end{array}\right)}}\underset{={\mathbf{V}}^T}{\underbrace{\left(\begin{array}{ccc}[r]2/3& 2/3& 1/3\\ -2/3& 1/3& 2/3\\ 1/3& -2/3& 2/3\\ 0& 0& 0\end{array}\right)}}.$$

$${\mathbf{X}}_1=\left(\begin{array}{ccc}[r]-1/8& 1/8& -1/6\\ 11/8& -11/8& -5/6\\ -9/8& 9/8& -1/6\\ -1/8& 1/8& 7/6\end{array}\right)$$

$${\mathbf{X}}_2\approx \left(\begin{array}{ccc}[r]0.0417& 0.2917& 0\\ 2.2083& -0.5417& 0\\ 0.2750& 0.9650& 0\\ -0.3667& -1.2867& 0\end{array}\right)$$

$${\mathbf{u}}_t^T{\mathbf{X}}_t$$

$$={\mathbf{u}}_t^T{\mathbf{X}}_{t-1}-\frac{({\mathbf{u}}_t^T{\mathbf{X}}_{t-1}{\mathbf{v}}_t){\mathbf{u}}_t^T{\mathbf{X}}_{t-1}}{{\mathbf{u}}_t^T{\mathbf{X}}_{t-1}{\mathbf{v}}_t}$$

$$=0$$

$${\mathbf{u}}_t^T{\mathbf{X}}_{t+s+1}$$

$$={\mathbf{u}}_t^T{\mathbf{X}}_{t+s}-\frac{({\mathbf{u}}_t^T{\mathbf{X}}_{t+s}){\mathbf{v}}_{t+s+1}{\mathbf{u}}_{t+s+1}{\mathbf{X}}_{t+s}}{{\mathbf{u}}_{t+s+1}^T{\mathbf{X}}_{t+s}{\mathbf{v}}_{t+s+1}}$$

$$=0-\frac{0{\mathbf{v}}_{t+s+1}{\mathbf{u}}_{t+s+1}{\mathbf{X}}_{t+s}}{{\mathbf{u}}_{t+s+1}^T{\mathbf{X}}_{t+s}{\mathbf{v}}_{t+s+1}}$$

$$=0$$

$${\mathbf{X}}_t^{\text{HD}}$$

$${\mathbf{X}}_t^{\text{PD}}$$

$${\mathbf{X}}_t^{\text{SD}}$$

$${\mathbf{U}}_t^T{\mathbf{X}}_t^{\text{HD}}{\mathbf{V}}_t$$

$$={\mathbf{U}}_t^T{\mathbf{X}}_{t-1}{\mathbf{V}}_t-{\mathbf{U}}_t^T{\mathbf{U}}_t({\mathbf{U}}_t^T{\mathbf{U}}_t{)}^{-1}{\mathbf{U}}_t{\mathbf{X}}_{t-1}{\mathbf{V}}_t({\mathbf{V}}_t^T{\mathbf{V}}_t{)}^{-1}{\mathbf{V}}_t^T{\mathbf{V}}_t$$

Full text

Department of Statistics

Rice University

Houston, TX 77005

michael.weylandt@rice.edu

1 Introduction

Principal Components Analysis (PCA, [42]) is a widely-used approach to finding low-dimensional patterns in complex data, enabling visualization, dimension reduction (compression), and predictive modeling. While PCA performs well in a wide range of low-dimensional settings, its performance degrades rapidly in high-dimensions, necessitating the use of regularized variants. Recently, [41] proposed Sparse and Functional PCA (SFPCA), a unified regularization scheme that allows for simultaneous smooth (functional) and sparse estimation of both row and column principal components (PCs). The rank-one SFPCA estimator is given by

[TABLE]

where $P_{\bm{u}}(\cdot)$ is a regularizer inducing sparsity in the row PCs, with strength controlled by $\lambda_{\bm{u}}$; $\bm{\Omega}_{\bm{u}}$ is a positive semi-definite penalty matrix, typically a second- or fourth-order difference matrix; $\bm{S}_{\bm{u}}=\bm{I}+\alpha_{\bm{u}}\bm{\Omega}_{\bm{u}}$ is a smoothing matrix for the row PCs, with strength controlled by $\alpha_{\bm{u}}$; and $\overline{\mathbb{B}}^{n}_{\bm{S}_{\bm{u}}}$ is the unit ellipse of the $\bm{S}_{\bm{u}}$-norm, i.e., $\overline{\mathbb{B}}^{n}_{\bm{S}_{\bm{u}}}=\{\bm{u}\in\mathbb{R}^{n}:\bm{u}^{T}\bm{S}_{\bm{u}}\bm{u}\leq 1\}$. (Respectively, $P_{\bm{v}}(\cdot)$, $\lambda_{\bm{v}}$, $\bm{\Omega}_{\bm{v}}$, $\alpha_{\bm{v}}$, $\bm{S}_{\bm{v}}$, and $\overline{\mathbb{B}}^{p}_{\bm{S}_{\bm{v}}}$ for the column PCs.)

[41] show that SFPCA unifies much of the existing regularized PCA literature [3, 4, 5, 6, 7, 8] into a single framework, avoiding many pathologies of other approaches. Finally, they propose an efficient alternating maximization scheme with guaranteed global convergence to solve the bi-concave SFPCA problem (1). The SFPCA estimator only allows for a single pair of PCs to be estimated for a given data matrix $\bm{X}$. [41] suggest applying SFPCA repeatedly to the deflated data matrix if multiple PCs are desired. While this approach performs acceptably in practice, it loses the interpretable orthogonality properties of standard PCA. In particular, the estimated PCs are no longer guaranteed to be orthogonal to each other, hindering the common interpretation of PCs as statistically independent sources of variance, or to the deflated data matrix, suggesting that additional signal remains uncaptured.

We extend the work of [41] to address these shortcomings: first, in Section 2, we modify the SFPCA estimator to simultaneously estimate several sparse PCs subject to orthogonality constraints. The resulting estimator is constrained to a product of generalized Stiefel mainfolds and we propose three efficient algorithms to solve the resulting manifold optimization problem. Next, in Section 3, we propose improved deflation schemes which provably remove all of the signal from the data matrix, allowing for more accurate iterative estimation of sparse PCs. Finally, we demonstrate the improved performance of our manifold estimators and deflation schemes in Section 4. Supplemental materials for this paper, including proofs, counter-examples, and additional algorithmic details, are available online at https://arxiv.org/abs/1907.12012.

2 Manifold Optimization for SFPCA

One of the most attractive properties of PCA is the factors it extracts are orthogonal ($\bm{u}_{t}^{T}\bm{u}_{s}=\bm{v}_{t}^{T}\bm{v}_{s}=1$ if $t=s$ and [math] otherwise). Because of this, PCs can be interpreted as separate sources of variance and, under an additional Gaussianity assumption, statistically independent. While this follows directly from the properties of eigendecompositions for standard PCA, it is much more difficult to obtain similar results for sparse PCA. Some authors have suggested that there exists a fundamental tension between orthogonality and sparsity, with [44] calling the goal of sparse and orthogonal estimation “questionable.” Indeed ex post orthogonalization of sparse PCs, e.g., using a Gram-Schmidt step, destroys any sparsity in the estimated PCs. To avoid this, it is necessary to impose orthogonality directly in the estimation step, rather than trying to impose it afterwards.

We modify the SFPCA estimator to simultaneously estimate multiple PCs subject to an orthognality constraint:

[TABLE]

where $\mathcal{V}^{\bm{S}_{\bm{u}}}_{n\times k}$ is the generalized Stiefel manifold of order $k$ over $\mathbb{R}^{n}$, i.e.,

[TABLE]

The generalized Stiefel manifold constraint ensures orthogonality of the estimated PCs, while still allowing us to capture most of the variability in the data.111When estimating orthogonal factors, it is common to re-express the problem using a (generalized) Grassmanian manifold constraint to avoid identifiability issues. We cannot use the Grassmanian approach here as sparse estimation (implicitly) fixes a single coordinate system. We note that, because we use a generalized Stiefel constraint, the estimated PCs will be orthogonal with respect to $\bm{S}_{\bm{u}}$, i.e., $\bm{u}_{t}^{T}\bm{S}_{\bm{u}}\bm{u}_{s}=0$ for $t\neq s$, rather than orthogonal in the standard sense. This is commonly observed for functional PCA variants [3, 5, 6] and can be interpreted as orthogonality under the inner product generating the $\bm{S}_{\bm{u}}$-norm. If no roughness penalty is imposed ($\alpha_{\bm{u}}=0$ or $\alpha_{\bm{v}}=0$), then our method gives orthogonality in the standard sense.

To solve the Manifold SFPCA problem (2), we employ an alternating maximization scheme, first holding $\bm{V}$ fixed while we update $\bm{U}$ and vice versa, as described in Algorithm 1. Even with one parameter held fixed, the resulting sub-problems are still difficult manifold optimization problems, which require iterative approaches to obtain a solution [62, 63, 60, 38, 28, 29, 45, 17].

[41] developed a custom projected + proximal gradient algorithm to solve the $\bm{u}$- and $\bm{v}$ subproblems of the rank-one SFPCA estimator. Assuming $P_{\bm{U}}$ and $P_{\bm{V}}$ are positive homogeneous, (e.g. $P(\cdot)=\|\bm{A}\cdot\|_{p}$ for arbitrary $p\geq 1$ and $\bm{A}$), they establish convergence to a stationary point. In order to extend this idea to the multi-rank (manifold) case, we use the recently-proposed Manifold Proximal Gradient (ManPG) scheme of [28], detailed in Algorithm 2. ManPG proceeds in two-steps: first, we solve for a descent direction $\bm{D}$ of the objective along the tangent space of the generalized Stiefel manifold, subject to the tangency constraint of $\bm{D}^{T}\bm{S}_{\bm{u}}\bm{U}^{(k)}$ being skew-symmetric; secondly, back-tracking line search is used to determine a step-size $\alpha$, after which the estimate is projected back onto the generalized Steifel manifold using a retraction anchored at the previous $\bm{U}^{(k)}$. The retraction, which plays the same role as the projection step in the original SFPCA algorithm [41, Algorithm 1], can be computed using a Cholesky factorization [17, Algorithm 3.1]. [29] showed that a single step of ManPG is sufficient to ensure convergence despite the bi-concave objective: we refer to their approach as Alternating ManPG (A-ManPG).

Note that in general manifold proximal gradient schemes [28, 29] impose a maximum step-size to ensure that linearization of the smooth portion of the objective actually leads to descent: because the smooth portion of our objective function in linear in $\bm{U}$ and $\bm{V}$, we can omit this term from Algorithm 2.

While efficient when tuned properly, we have found the performance of ManPG on the $\bm{U}$- and $\bm{V}$-subproblems quite sensitive to infeasibility in the descent direction. A more robust scheme can be derived by using the Manifold ADMM (MADMM) scheme of [38] to solve the subproblems. Like standard ADMM schemes, MADMM allows us to split a problem into two parts, each of which can be solved more easily than the global problem. When applied to the $\bm{U}$- and $\bm{V}$-subproblems, the MADMM allows us to separate the manifold constraint from the sparsity inducing regularizer, thereby side-stepping the orthogonality / sparsity tension at the heart of this paper. After this splitting, the smooth update can be shown to equivalent to the unbalanced Procrustes problem [26, 27] with a closed-form update: $\bm{U}^{(k+1)}=\bm{S}_{\bm{u}}^{-1/2}\bm{A}\bm{B}^{T}$ where $\bm{A}\bm{\Delta}\bm{B}^{T}$ is the SVD of $\bm{S}_{\bm{u}}^{-1/2}\bm{X}\hat{\bm{V}}+\rho\bm{S}_{\bm{u}}^{1/2}(\bm{W}^{(k)}-\bm{Z}^{(k)})$. The sparse update is simply the proximal operator of $P_{\bm{U}}(\cdot)$, typically a threhsolding step [20, Chapter 6]. To the best of our knowledge, the convergence of MADMM has not yet been established, but we have not observed significant non-convergence problems in our experiments.

3 Iterative Deflation for SFPCA

We next consider the use of iterative deflation schemes for multi-rank SFPCA. As discussed by [43], the attractive orthogonality properties of standard (Hotelling’s) deflation depend critically on the estimated PCs being exact eigenvectors of the covariance matrix. Because the PCs estimated by sparse PCA schemes are almost surely not eigenvectors, [43] proposes several alternate deflation schemes which retain some of the attractive properties of Hotelling’s deflation even when non-eigenvectors are used. We extend these to the low-rank model and allow for deflation by several PCs, possibly with non-unit norm, simultaneously, e.g., as produced by ManSFPCA. To ease exposition, we first work in the vector setting and consider the general case at the end of this section. The properties of our proposed deflation schemes are summarized in Table 1 above.

The simplest deflation scheme is essentially that used by [42], extended to the low-rank model:

[TABLE]

For two-way sparse PCA variants [7, 8], this deflation gives a deflated matrix which is “two-way” orthogonal to the estimated PCs, i.e., $\bm{u}_{t}^{T}\bm{X}_{t}\bm{v}_{t}=0$. We may interpret this as Hotelling’s deflation (HD) capturing all the signal jointly associated with the pair $(\bm{u}_{t},\bm{v}_{t})$.

We may also ask if HD captures all of the signal associated with $\bm{u}_{t}$ or only the signal which is also associated with $\bm{v}_{t}$. If HD captures all of the signal associated with $\bm{u}_{t}$, then we would expect $\bm{u}_{t}^{T}\bm{X}_{t}\tilde{\bm{v}}=0$ for all $\tilde{\bm{v}}\in\mathbb{R}^{p}$, or equivalently, $\bm{u}_{t}^{T}\bm{X}_{t}=\mathbf{0}_{p}$. Interestingly, HD does not have this left-orthogonality property, suggesting that it leaves additional $\bm{u}_{t}$-signal in the deflated matrix $\bm{X}_{t}$.

HD fails to yield left- and right-orthogonality because it is not based on a projection operator. To address this in the covariance model, [43] proposed a deflation scheme which projects the covariance matrix onto the orthogonal complement of the estimated principal component. We extend this idea to the low-rank model by projecting the column- and row-space of the data matrix into the orthogonal complement of the left- and right-PCs respectively, giving two-way projection deflation (PD):

[TABLE]

Unlike HD, PD captures all of the linear signal associated with $\bm{u}_{t}$ and $\bm{v}_{t}$ individually ($\bm{u}_{t}\bm{X}_{t}\tilde{\bm{v}}=\tilde{\bm{u}}\bm{X}_{t}\bm{v}_{t}=0,\forall\tilde{\bm{u}}\in\mathbb{R}^{n},\tilde{\bm{v}}\in\mathbb{R}^{p}$).

If we use PD repeatedly, however, the multiply-deflated matrix will not continue to be orthogonal to the PCs: that is, $\bm{u}_{t}^{T}\bm{X}_{t+s}\neq\mathbf{0}$ for $s\geq 1$. This suggests that repeated application of PD can reintroduce signal in the direction of the PCs that we previously removed. This occurs because PD works by sequentially projecting the data matrix, but in general the compositition of two orthogonal projections is not another orthogonal projection without additional assumptions. To address this, [43] proposed a Schur complement deflation (SD) technique, which we now extend to the low-rank (two-way) model:

[TABLE]

While [43] motivates this approach using conditional distributions and a Gaussianity assumption on $\bm{X}$, it can also be understood as an alternate projection construct which is more robust to scaling and non-orthogonality.

So far, we have only considered the behavior of the proposed deflation schemes for two-way sparse PCA. If we consider SFPCA in generality, however, $\bm{u}_{t}$ and $\bm{v}_{t}$ are unit vectors under the $\bm{S}_{\bm{u}}$- and $\bm{S}_{\bm{v}}$-norms, not under the Euclidean norm. Consequently, the projections used by PD may not be actual projections and a PD deflated matrix may fail to be two- or one-way orthogonal. Normalizing the estimated PCs before deflation addresses this problem and is recommended in practice: conversely, because its deflation term is invariant under rescalings of $\bm{u}_{t}$ and $\bm{v}_{t}$, SD works without renormalization.

The extension of these techniques to the multi-rank case is straightforward. We give the normalized variants here:

[TABLE]

As in the covariance model, if $\bm{u}_{t}$ and $\bm{v}_{t}$ are true singular vectors, all three deflation schemes are equivalent.

4 Simulation Studies

In this section, we compare the performance of Manifold SFPCA and the iterative rank-one deflation schemes proposed above in illustrative simulation studies. Manifold SFPCA using Manifold ADMM (Algorithm 3) to solve the subproblems achieves better solutions than Manifold Proximal Gradient or A-ManPG (Algorithm 2) in less time. Furthermore, despite the additional flexibility of the iterative rank-one variants, Manifold SFPCA achieves both better signal recovery and a higher proportion of variance explained, even when the orthogonality assumptions are violated.

We first consider the relative performance of the three algorithms proposed for solving the Manifold SFPCA problem (2). We generate data in $\bm{X}^{*}=\bm{U}^{*}\bm{D}^{*}\bm{V}^{*}\in\mathbb{R}^{250\times 100}$ with three distinct PCs: the left PCs ($\bm{U}^{*}$) are localized sinusoids of varying frequency; the right PCs ($\bm{V}^{*}$) are non-overlapping sawtooth waves. (See Figure 1.) We add independent standard Gaussian noise ($\bm{E}$) to give a signal-to-noise ratio (SNR) of $\|\bm{X}^{*}\|/\|\bm{E}\|\approx 1.2$. We fix $\lambda_{\bm{u}}=\lambda_{\bm{v}}=1$ and $\alpha_{\bm{u}}=\alpha_{\bm{v}}=3$ which is near optimal for all three schemes. This is a favorable setting for Manifold SFPCA as the underlying signals are orthogonal, sparse, smooth, and of comparable magnitude.

Table 2 shows the performance of our three Manifold SFPCA algorithms on several metrics, averaged over 100 replicates. Overall, the Manifold ADMM (MADMM [38]) and Alternating Manifold Proximal Gradient (A-ManPG [29]) variants perform best, handily beating the Manifold Proximal Gradient scheme [28] on all measures. MADMM achieved the best objective value on every replicate.

In terms of signal recovery, MADMM achieves slightly better performance than A-ManPG on the right singular vectors, while A-ManPG is slightly better on the left singular vectors. Computationally, MADMM dominates both proximal gradient variants even though it requires many more matrix decompositions. The descent direction subproblems of ManPG and A-ManPG are rather expensive to solve repeatedly and their performance is very sensitive to the solver used. Overall, the MADMM variant of Manifold SFPCA achieves the best optimization and statistical performance in far less time than the proximal gradient-based variants.

Next, we compare Manifold SFPCA with the iterative deflation schemes proposed in Section 3 in two different scenarios: the favorable scenario used above and a less-favorable scenario where the true PCs are shifted and no longer orthogonal. ($n=p=100$ and $\|(\bm{U}^{*})^{T}(\bm{U}^{*})-\bm{I}_{3}\|,\|(\bm{V}^{*})^{T}\bm{V}^{*}-\bm{I}_{3}\|\approx 0.37$.) We compare the proportion of variance explained using Manifold SFPCA with iterative rank-one SFPCA using the normalized Hotelling, Projection, and Schur Complement deflation strategies. As can be seen in Table 3, PD and SD consistently dominate HD. Because PD and SD fully remove the signal associated with estimated PCs, the subsequent PCs are able to capture different signals and explain a larger fraction of variance. By ensuring that the signal is never re-introduced, SD does even better than PD as we consider higher ranks. Interestingly, while PD and SD perform about as well when the underlying signals are orthogonal, SD performs much better in the non-orthogonal scenario. By estimating all three PCs simultaneously, Manifold SFPCA is able to find a better set of PCs than any of the greedy deflation methods.

5 Discussion

We have introduced two practical extensions to Sparse and Functional PCA: first, we presented a multi-rank scheme which estimates multiple PCs simultaneously and proposed algorithms to solve the resulting manifold optimization problems. The resulting estimator inherits many of the attractive properties of rank-one SFPCA and is, to the best of our knowledge, the first multi-rank PCA scheme for the low-rank model. ManSFPCA combines both the flexiblity and superior statistical performance of rank-one SFPCA with the superior interpretability of orthogonal (non-regularized) PCA.

Secondly, we re-considered the use of Hotelling’s deflation, and developed two additional deflation schemes which have attractive theoretical properties and emprical performance. Our schemes extend the results of [43] in several ways: they allow for deflation by multiple PCs in a single step, they are applicable to the low-rank model and the covariance model, and they are robust to non-orthogonality and non-unit-scaling common to Functional PCA variants. While developed for SFPCA, these deflation schemes are useful for any regularized PCA model.

We note here that our results can be extended naturally to other multivariate analysis techniques which can be expressed in a regularized SVD framework (e.g., PLS, CCA, etc.). Our deflation approaches can be also extended to the higher-order / multi-way context and may be particularly useful in the context of regularized CP decompositions [22, 23]. In the tensor setting, ManSFPCA is a sparse and smooth version of a Tucker decomposition, which suggests several interesting extensions we leave for future work [24, 25].

6 References