Statistical Analysis of Karcher Means for Random Restricted PSD Matrices

TL;DR

This paper provides a non-asymptotic statistical analysis of the Karcher mean on the manifold of restricted positive semi-definite matrices, with applications to distributed PCA algorithms and strong numerical validation.

Contribution

It offers the first non-asymptotic analysis of the Karcher mean on restricted PSD matrices and applies it to distributed PCA performance guarantees.

Findings

Deterministic error bounds for the Karcher mean under noise.

Distributed PCA matches full sample PCA performance.

Numerical experiments support theoretical results.

Abstract

Non-asymptotic statistical analysis is often missing for modern geometry-aware machine learning algorithms due to the possibly intricate non-linear manifold structure. This paper studies an intrinsic mean model on the manifold of restricted positive semi-definite matrices and provides a non-asymptotic statistical analysis of the Karcher mean. We also consider a general extrinsic signal-plus-noise model, under which a deterministic error bound of the Karcher mean is provided. As an application, we show that the distributed principal component analysis algorithm, LRC-dPCA, achieves the same performance as the full sample PCA algorithm. Numerical experiments lend strong support to our theories.

Equations

$$\tilde{\mathbf{z}}={\mathrm{argmin}}_{\mathbf{z}\in \mathcal{M}}\sum _m{d}^2(\mathbf{z},{\mathbf{z}}_m).$$

$$\mathfrak{h}:S^{\ast }(p,K)\mapsto {\mathcal{L}}^{\ast }(p,K),\mathbf{A}\mapsto \mathbf{N},$$

$$\mathfrak{g}:{\mathcal{L}}^{\ast }(p,K)\mapsto \mathcal{L}(p,K),\mathbf{N}\mapsto {\mathbf{N}}^{\prime },$$

$$\tilde{\mathbf{A}}={\mathfrak{h}}^{-1}\circ {\mathfrak{g}}^{-1}(\frac{1}{M}\sum _{m=1}^M\mathfrak{g}\circ \mathfrak{h}({\mathbf{A}}^m)).$$

$${\mathbf{A}}^m={\mathfrak{h}}^{-1}\circ {\mathfrak{g}}^{-1}({\mathbf{N}}^{\prime }+{\mathbf{E}}^m),m=1,\dots ,M,$$

$$\tilde{\mathbf{A}}={\mathfrak{h}}^{-1}\circ {\mathfrak{g}}^{-1}({\mathbf{N}}^{\prime }+\frac{1}{M}\sum _{m=1}^M{\mathbf{E}}^m).$$

$$\Vert \tilde{\mathbf{A}}-\mathbf{A}{\Vert }_{\mathrm{F}}\le \sqrt{\frac{c_{2}pK{\sigma }^2}{M}}$$

$${\mathbf{A}}^m=(\mathbf{N}+{\mathbf{E}}^m)(\mathbf{N}+{\mathbf{E}}^m{)}^{\top },$$

$$=\mathbf{Q}+f_{\mathbf{R}}(\mathbf{E}{\mathbf{Q}}^{\top })\mathbf{Q}+{\mathcal{O}}_{\max }({ϵ}_0^2),$$

$$=\mathbf{R}+\mathbf{E}{\mathbf{Q}}^{\top }-\mathbf{R}{f}_{\mathbf{R}}(\mathbf{E}{\mathbf{Q}}^{\top })+{\mathcal{O}}_{\max }({ϵ}_0^2),$$

$$f_{\mathbf{R}}(\mathbf{E})=\mathcal{U}({\mathbf{R}}^{-1}\mathbf{E})-(\mathcal{U}({\mathbf{R}}^{-1}\mathbf{E}){)}^{\top },$$

$$\mathcal{U}(\mathbf{P}{)}_{ij}={\mathbf{P}}_{ij},i<j,\mathcal{U}(\mathbf{P}{)}_{ij}=0,\text{otherwise}.$$

$$\tilde{\mathbf{N}}$$

$$+{\mathcal{O}}_{\max }({ϵ}_0^2),$$

$$\Vert \tilde{\mathbf{N}}-\mathbf{N}{\Vert }_F\le \mathcal{O}(\Vert \frac{1}{M}\sum _{m=1}^M{\mathbf{E}}^m{\Vert }_{\text{F}})+\mathcal{O}(pK{ϵ}_0^2).$$

$$\mathbb{E}[e^{\lambda ⟨\mathbf{u},\mathbf{z}⟩}]\le e^{\frac{{\lambda }^2{\sigma }^2}{2}},\forall \lambda \in \mathbb{R},\forall \mathbf{u}\in {\mathbb{R}}^p,\Vert \mathbf{u}{\Vert }_2=1.$$

$${\hat{\mathbf{H}}}^m={\mathrm{argmin}}_{\mathbf{O}\in {\mathcal{O}}_{K\times K}}\Vert {\hat{\mathbf{V}}}^m\mathbf{O}-\mathbf{V}{\Vert }_{\text{F}}$$

$${\hat{\mathbf{V}}}^m{\hat{\mathbf{H}}}^m=\mathbf{V}+g({\mathcal{E}}^m\mathbf{V})+{\mathcal{O}}_{\text{F}}({ϵ}^2),$$

$${\hat{\mathbf{E}}}^m={\mathcal{E}}^m\mathbf{V}{\mathbf{Q}}^{\ast }+\mathbf{\Sigma }g({\mathcal{E}}^m\mathbf{V}){\mathbf{Q}}^{\ast }+{\mathcal{O}}_F({ϵ}^2).$$

$$\Vert \frac{1}{M}\sum _{m=1}^M{\hat{\mathbf{E}}}^m{\Vert }_{\text{F}}\le C\Vert \frac{1}{M}\sum _{m=1}^M{\mathcal{E}}^m{\Vert }_2+\mathcal{O}({ϵ}^2)$$

$$\underset{m}{\max }\Vert {\hat{\mathbf{E}}}^m{\Vert }_{\max }\le C_3\sqrt{\frac{\log (p)}{n}}+{\delta }_1+{\delta }_2,$$

$$\underset{m}{\max }\Vert {\hat{\mathbf{E}}}^m{\Vert }_{\max }\le C\sqrt{\frac{\log (pM)}{n}},$$

$$\Vert \tilde{\mathbf{N}}-\mathbf{N}{\Vert }_{\text{F}}\le \mathcal{O}\left(\Vert \frac{1}{M}\sum _{m=1}^M{\mathcal{E}}^m{\Vert }_2\right)+\mathcal{O}({ϵ}^2)+\mathcal{O}(\sqrt{p}{ϵ}_0^2)$$

$$\Vert \tilde{\mathbf{N}}-\mathbf{N}{\Vert }_{\text{F}}\le \mathcal{O}\left(\frac{\log (p)\sqrt{r}}{\sqrt{Mn}}\right)+\mathcal{O}\left(\frac{{\tilde{r}}_1}{n}\right)+\mathcal{O}\left(\frac{{\tilde{r}}_2}{n^2}\right).$$

$$\Vert \tilde{\mathbf{N}}-\mathbf{N}{\Vert }_{\text{F}}\le \mathcal{O}\left(\frac{\log (p)\sqrt{r}}{\sqrt{Mn}}\right).$$

$$\Vert \frac{1}{M}\sum _{m=1}^M{\mathbf{E}}^m{\Vert }_{\text{F}}^2\le \frac{pK{\sigma }^2}{M}(1+t),$$

$$\underset{i=1,\dots ,K}{\max }|\frac{1}{M}\sum _{m=1}^M{\mathbf{E}}_{ii}^m{|}^2\le \frac{{\sigma }^2}{M}(1+t)$$

$$|\exp (\frac{1}{M}\sum _{m=1}^M{\mathbf{E}}_{ii}^m)-1|\le 2|\frac{1}{M}\sum _{m=1}^M{\mathbf{E}}_{ii}^m|,$$

$$\Vert \tilde{\mathbf{A}}-\mathbf{A}{\Vert }_F\le \sqrt{\frac{c_{2}pK{\sigma }^2}{M}},$$

$${\check{\mathbf{Q}}}_{ii}^{ij}={\check{\mathbf{Q}}}_{jj}^{ij}=\cos ({\theta }^{ij}),{\check{\mathbf{Q}}}_{ij}^{ij}=-{\check{\mathbf{Q}}}_{ji}^{ij}=\sin ({\theta }^{ij}),{\check{\mathbf{Q}}}_{kk}^{ij}=1,\forall k\ne i,j,{\check{\mathbf{Q}}}_{kl}^{ij}=0,\text{otherwise},$$

$${\check{\mathbf{R}}}^{(1)}=(\mathbf{R}+\mathbf{E})({\mathbf{I}}_K+{\check{\mathbf{Q}}}^{(1{)}^{\top }}-{\mathbf{I}}_K)=\mathbf{R}+\mathbf{E}+\mathbf{R}({\check{\mathbf{Q}}}^{(1{)}^{\top }}-{\mathbf{I}}_K)+{\mathcal{O}}_{\max }({ϵ}_0^2).$$

$${\check{\mathbf{Q}}}^{ij}={\mathbf{I}}_K+{\check{\mathbf{P}}}^{ij}+{\mathcal{O}}_{\max }({ϵ}_0^2),$$

Full text

1 Introduction

Positive semi-definite (PSD) matrices arise in a wide range of applications, such as covariance matrices in statistics (Wainwright,, 2019), kernel matrices in machine learning (Hastie et al.,, 2009), diffusion tensor images in medical imaging (Dryden et al.,, 2009), semi-definite programming (Journée et al.,, 2010), and covariance descriptors in image set classification (Wang et al.,, 2012), to name a few. From a geometric perspective, the cone of PSD matrices is not a vector space, since linear combinations of multiple PSD matrices are not necessarily PSD matrices. Instead, the set of (restricted) PSD matrices of fixed rank has been endowed with different metrics such that it forms a Riemannian manifold (Bonnabel and Sepulchre,, 2010; Vandereycken et al.,, 2013; Massart and Absil,, 2020; Neuman et al.,, 2021). By utilizing the geometric structures, researchers have developed many powerful statistical or computational methods (Faraki et al.,, 2016; Cornea et al.,, 2017; Patrangenaru and Ellingson,, 2016).

One important concept in Riemannian geometry or more generally metric spaces is the Karcher mean (Karcher,, 1977). The Karcher mean is often referred to as the Fr$\rm\acute{e}$chet mean or the barycenter of mass. Given $M$ points $\{\bm{z}_{m}\}_{m=1}^{M}$ on a metric space $(\mathcal{M},d)$ with distance function $d(\cdot,\cdot)$, the Karcher mean $\widetilde{\bm{z}}$ of these points is given by

[TABLE]

When the underlying space is Euclidean, the Karcher mean is reducesd to the arithmetic mean. In general, the existence and computation of the Karcher mean is already complicated due to the possibly intricate non-Euclidean structure (Karcher,, 1977; Bini and Iannazzo,, 2013). As a result, most works focus on the computation and applications of the Karcher mean, while few provide statistical guarantees. Statistically, Bhattacharya and Patrangenaru, (2003, 2005) establish a large sample theory of the Karcher mean on manifolds with applications to spheres and projective spaces. Bigot and Gendre, (2013) shows the minimax optimality of the Karcher mean of discretely sampled curves. In this paper, we consider the manifold of restricted PSD matrices by Neuman et al., (2021). In particular, we first study an intrinsic mean model, inspired by the geometric structure of the restricted PSD manifold. A non-asymptotic statistical analysis of the Karcher mean is provided under this intrinsic model. We further consider a general extrinsic signal-plus-noise model, which does not necessarily coincide with the manifold geometry by Neuman et al., (2021). For this general model, we give a deterministic error bound for the Karcher mean, which is then used to provide an error bound for a distributed principal component analysis algorithm.

The Karcher mean is closely related to distributed learning problems, especially the divide-and-conquer (DC) framework (Mackey et al.,, 2011). In distributed learning problems, massive datasets are scattered across distant servers and directly fusing these datasets is challenging due to concerns on communication cost, privacy, data security, and ownership, among others. A commonly used distributed framework is the DC framework which first computes local estimators locally and then aggregate them on the central server, where the last step is often equivalent to computing the Karcher mean on certain manifolds. For example, the divide-and-conquer principal component analysis (PCA) algorithms (Fan et al.,, 2019; Bhaskara and Wijewardena,, 2019; Neuman et al.,, 2021) essentially compute the Karcher means on the Grasssmann manifold, Euclidean space, or the manifold of restricted PSD matrices, respectively. Motivated by this observation, we give theoretical guarantees of the DC PCA algorithm, LRC-dPCA, proposed in Neuman et al., (2021) by applying our non-asymptotic statistical analysis of the Karcher mean on the restricted PSD manifold. Specifically, we show that given sufficienly large local sample size, LRC-dPCA achieves the same performance as the full sample PCA algorithm, which outputs the top eigenvectors of the covariance matrix based on full data.

Our contributions are three-fold. First, we provide a non-asymptotic statistical analysis of the Karcher mean on the restricted PSD manifold under an intrinsic model. Second, for a generic signal-plus-noise model, we give a deterministic characterization of the Karcher mean and then obtain a deterministic error bound. Third, as an application, we show that LRC-dPCA and full sample PCA share the same performance given sufficiently large local sample size. Numerical experiments are carried out to support our theories.

The rest of this paper proceeds as follows. We conclude this section with a discussion on related works. Section 2 reviews the geometry for restricted PSD matrices proposed in Neuman et al., (2021). Then in Section 3, we provide the theoretical analysis of the Karcher mean on the restricted PSD manifolds. Applications to distributed PCA algorithms are given in Section 4. Numerical experiments are carried out in Section 5 and we give concluding remarks in Section 6. Proofs are left to the Appendix.

1.1 Related work

Manifolds of PSD matrices The cone of symmetric positive definite (SPD) matrices is not a vector space. It can be viewed as different Riemannian manifolds when endowed with different metrics, such as the affine-invariant metric (Moakher,, 2005) and the Log-Euclidean metric (Arsigny et al.,, 2007). It is, however, non-trivial to generalize these metrics to the rank-deficient (PSD) case. To this end, Bonnabel and Sepulchre, (2010) treated a PSD matrix of rank $K$ in a quotient space as a $K$-dimensional subspace coupled with a $K$-by-$K$ SPD matrix and then endowed the manifold of PSD matrices with a weighted product metric. Using this geometry, Bonnabel et al., (2013) developed a rank-preserving geometric mean of PSD matrices. Later, Vandereycken et al., (2013) viewed a PSD manifold as a homogeneous space and Massart and Absil, (2020) analyzed a quotient geometry on the manifold of PSD matrices. However, it is hard to give a statistical model on these manifolds. More recently, Neuman et al., (2021) proposed a geometry for restricted PSD matrices which has closed-form solutions for many geometric concepts including the Karcher mean. Our paper provides statistical analysis of the Karcher mean corresponding to this geometry.

Distributed PCA To estimate the leading eigenvector, Garber et al., (2017) proposed a sign-fixing averaging approach. To estimate the top $K$ eigenspace, Fan et al., (2019) proposed a projector averaging approach and Charisopoulos et al., (2021) proposed to average local eigenvector matrices after carefully rotating them. Disregarding the information of eigenvalues, both Fan et al., (2019)’s and Charisopoulos et al., (2021)’s methods require the knowledge of the precise location $K$ of a large eigen gap. To alleviate this issue, Bhaskara and Wijewardena, (2019) proposed to average the local rank-$K$ approximation matrices and then conduct PCA on the aggregated matrix. Neuman et al., (2021) utilized the same methods as Bhaskara and Wijewardena, (2019) except that the average of local rank-$K$ approximation matrices is taken on the manifold of restricted PSD matrices. Neuman et al., (2021) did not provide statistical analysis for their proposed method, while our paper fixes this gap as an application of the main results. Another branch of research turns PCA into the problem of solving a linear system and then solves distributed PCA by some multi-round algorithms. Among them, some make use of the shift-and-invert framework (Garber et al.,, 2017; Chen et al.,, 2021), while some use incremental update schemes (Gang et al.,, 2019; Grammenos et al.,, 2020; Li et al.,, 2021).

Notation. By convention, we use regular letters for scalars and bold letters for both vectors and matrices. Given a vector $\bm{u}\in\mathbb{R}^{p}$, denote by $\|\bm{u}\|_{2}$ its $\ell_{2}$ norm. Given a matrix $\bm{A}\in\mathbb{R}^{n\times p}$, we use $\|\bm{A}\|_{\textrm{F}}$, $\|\bm{A}\|_{2}$ and $\|\bm{A}\|_{\max}=\max_{i,j}|\bm{A}_{ij}|$ to denote its Frobenius norm, $\ell_{2}$ norm and max norm, respectively. We use $\textnormal{span}(\bm{A})$ to represent the subspace spanned by the columns of $\bm{A}$. For a symmetric matrix $\bm{A}$, denote by $\lambda_{j}(\bm{A})$ its $j$th largest eigenvalue. For two sequences of real numbers $\left\{a_{n}\right\}_{n\geq 1}$ and $\left\{b_{n}\right\}_{n\geq 1}$, we write $a_{n}\lesssim b_{n}$ (or $a_{n}\gtrsim b_{n}$) if $a_{n}\leq Cb_{n}$ (or $a_{n}\geq Cb_{n}$) for some constant $C>0$ independent of $n$. For an infinitesimal number $\epsilon$, we denote a matrix whose Frobenius norm or max norm is $\mathcal{O}(\epsilon)$ (i.e., $\lesssim\epsilon$) by $\mathcal{O}_{\textrm{F}}(\epsilon)$ or $\mathcal{O}_{\max}(\epsilon)$, respectively. Given a random variable $x\in\mathbb{R}$, we define $\|x\|_{\psi_{2}}=\sup_{p\geq 1}(\mathbb{E}|x|^{p})^{1/p}/\sqrt{p}$ and $\|x\|_{\psi_{1}}=\sup_{p\geq 1}(\mathbb{E}|x|^{p})^{1/p}/p$. Given two integers $p\geq K>0$, we denote by $\mathcal{O}_{p\times K}$ the set of matrices in $\mathbb{R}^{p\times K}$ whose columns are orthonormal. Denote by $S(p,K)$ the set of all $p\times p$ PSD matrices of rank $K$. Denote by $a\vee b=\max\{a,b\}$.

2 The Manifold of Restricted PSD Matrices

In this section, we briefly recap the geometry for restricted PSD matrices (Neuman et al.,, 2021). To start with, any PSD matrix $\bm{A}\in S(p,K)$ has a unique Cholesky decomposition $\bm{A}=\bm{L}\bm{L}^{\top}$ such that $\bm{L}\in\mathbb{R}^{p\times p}$ is a lower triangular matrix and has precisely $K$ positive diagonal elements and $p-K$ zero columns. The $j$th column of $\bm{L}$ is zero if and only if the $j$th column of $\bm{A}$ is linearly dependent on the previous $j-1$ columns of $\bm{A}$. Thus, we can rewrite $\bm{A}=\bm{N}\bm{N}^{\top}$, where $\bm{N}\in\mathbb{R}^{p\times K}$ consists of $K$ non-zero columns of $\bm{L}$ without changing the order. Note that $\bm{N}$ is mock lower triangular, i.e., $\bm{N}_{ij}=0$ if $i<j$. We refer to $\bm{N}$ as the reduced Cholesky factor of $\bm{A}$. To further develop a geometric structure, Neuman et al., (2021) consider the restricted subset $S^{*}(p,K)$ of $S(p,K)$ such that the first $K$ columns of $\bm{A}\in S^{*}(p,K)$ are linearly independent. The set of all reduced Cholesky factors of matrices in $S^{*}(p,K)$ is denoted by $\mathcal{L}^{*}(p,K)$, which is equivalent to the set of all mock lower triangular matrices in $\mathbb{R}^{p\times K}$ with positive diagonal elements. Neuman et al., (2021) impose a Riemannian structure on $S^{*}(p,K)$ and $\mathcal{L}^{*}(p,K)$ such that the following mappings are isometric,

[TABLE]

where $\bm{N}=\mathfrak{h}(\bm{A})$ is the reduced Cholesky factor of $\bm{A}$, $\mathcal{L}(p,K)=\{\bm{N}^{\prime}\in\mathbb{R}^{p\times K}:\bm{N}^{\prime}_{ij}=0,i<j\}$ is endowed with a Euclidean structure, and $\bm{N}^{\prime}=\mathfrak{g}(\bm{N})\in\mathcal{L}(p,K)$ is defined by $\bm{N}^{\prime}_{ii}=\log(\bm{N}_{ii}),\forall\,i$ and $\bm{N}^{\prime}_{ij}=\bm{N}_{ij},\forall\,i>j$. We refer to $\bm{N}^{\prime}=\mathfrak{g}\circ\mathfrak{h}(\bm{A})$ as the reduced log-Cholesky factor of $\bm{A}$. The Karcher mean $\widetilde{\bm{A}}$ of $M$ restricted PSD matrices $\{\bm{A}^{m}\}_{m=1}^{M}\subset S^{*}(p,K)$ has a closed-form solution, which is given by

[TABLE]

The algorithm computing $\widetilde{\bm{A}}$ is referred to as the Low Rank Cholesky (LRC) algorithm (Neuman et al.,, 2021).

3 Statistical Analysis of the Karcher Mean

In this section, we provide the first statistical analysis of the Karcher mean under an intrinsic model on the restricted PSD manifold. Then we consider a general signal-plus-noise model under which a deterministic error bound of the Karcher mean is given.

3.1 An intrinsic model

Inspired by the isometry stated in equations (2.1) and (2.2) between the manifold $S^{*}(p,K)$ of restricted PSD matrices and the Euclidean space $\mathcal{L}(p,K)$, we propose the following intrinsic model. Suppose $\bm{A}\in S^{*}(p,K)$ is the signal matrix and denote by $\bm{N}^{\prime}=\mathfrak{g}\circ\mathfrak{h}(\bm{A})$ its reduced log-Cholesky factor. The observations $\{\bm{A}^{m}\}_{m=1}^{M}$ are generated as follows:

[TABLE]

where $\{\bm{E}^{m}\}_{m=1}^{M}\subset\mathcal{L}(p,K)$ are independent and the lower triangular entries of $\bm{E}^{m}$ are independent normal variables with mean zero and variance $\sigma^{2}$. Under this intrinsic model, the Karcher mean of $\{\bm{A}^{m}\}_{m=1}^{M}$ can be rewritten as

[TABLE]

Using measure concentration, we can obtain a non-asymptotic error bound for the Karcher mean $\widetilde{\bm{A}}$.

Theorem 3.1** (Intrinsic Model).**

Suppose $\bm{A}\in S^{*}(p,K)$ is the signal matrix and assume $\|\bm{A}\|_{2}\leq C$ for some constant $C>0$. Assume samples $\{\bm{A}^{m}\}_{m=1}^{M}$ are generated from the intrinsic model (3.1) and denote by $\widetilde{\bm{A}}$ the Karcher mean of $\{\bm{A}^{m}\}_{m=1}^{M}$. Then there exist some constants $c_{1},c_{2}>0$ such that the following inequality

[TABLE]

holds with probability at least $1-e^{-c_{1}pK}$.

Remark 3.2**.**

It is worth noting that (3.3) achieves the optimal rate $M^{-1/2}$. In addition, it only depends on the intrinsic dimension $\mathcal{O}(pK)$ of the manifold, which can be much smaller than the ambient dimension $p^{2}$.

3.2 A general signal-plus-noise model

The intrinsic model may be too restricted, so this subsection introduces a general signal-plus-noise model and then provides a deterministic characterization of the Karcher mean. An application of this deterministic error bound to the distributed PCA problem will be given in Section 4. Similar to the intrinsic model, we denote by $\bm{A}\in S^{*}(p,K)$ the signal matrix and $\bm{N}=\mathfrak{h}(\bm{A})$ its reduced Cholesky factor. The observations $\{\bm{A}^{m}\}_{m=1}^{M}\subset S^{*}(p,K)$ are given by

[TABLE]

where $\bm{E}^{m}\in\mathbb{R}^{p\times K}$ represents the $m$-th noise matrix. Here $\bm{E}^{m}$ is not necessarily a mock lower triangular matrix, so the model is quite general. Also, the reduced Cholesky factor of $\bm{A}^{m}$ is not necessarily $\bm{N}+\bm{E}^{m}$, but rather $(\bm{N}+\bm{E}^{m})\bm{Q}^{m}$ for some orthogonal matrix $\bm{Q}^{m}\in\mathcal{O}_{K\times K}$. Denote by $\bm{N}^{m}$ the reduced Cholesky factor of $\bm{A}^{m}$.

To characterize the Karcher mean (2.3) of $\{\bm{A}^{m}\}_{m=1}^{M}$, we first establish a linear perturbation expansion of QR decomposition below.

Lemma 3.3** (Linear Perturbation Expansion).**

Suppose $\bm{Q}\in\mathcal{O}_{K\times K}$ and $\bm{R}\in\mathbb{R}^{K\times K}$ is a lower triangular matrix with positive diagonal elements. Given a noise matrix $\bm{E}\in\mathbb{R}^{K\times K}$, there exist a unique orthogonal matrix $\widecheck{\bm{Q}}\in\mathcal{O}_{K\times K}$ and a lower triangular matrix $\widecheck{\bm{R}}\in\mathbb{R}^{K\times K}$ with non-negative diagonal elements such that $\widecheck{\bm{R}}\widecheck{\bm{Q}}=\bm{R}\bm{Q}+\bm{E}$. When $\epsilon_{0}=\|\bm{E}\|_{\max}$ is sufficiently small, we have

[TABLE]

where $f_{\bm{R}}:\mathbb{R}^{K\times K}\mapsto\mathbb{R}^{K\times K}$ is given by

[TABLE]

It is worth emphasizing the following properties of $f_{\bm{R}}$. First, $f_{\bm{R}}$ is linear in its argument, i.e., $f_{\bm{R}}(a\bm{E}+b\bm{F})=af_{\bm{R}}(\bm{E})+bf_{\bm{R}}(\bm{F})$ for any $a,b\in\mathbb{R}$ and $\bm{E},\bm{F}\in\mathbb{R}^{K\times K}$. Second, $f_{\bm{R}}(\bm{E})$ is a skew-symmetric matrix, i.e., $(f_{\bm{R}}(\bm{E}))^{\top}=-f_{\bm{R}}(\bm{E})$. Last, $f_{\bm{R}}$ is bounded in the sense that $\|f_{\bm{R}}(\cdot)\|_{\textrm{F}}\leq\sqrt{2}\|\bm{R}^{-1}\|_{2}\|\cdot\|_{\textrm{F}}$. Similar first-order perturbation theories exist in the literature for QR, Cholesky, and LU factorization (Chang et al.,, 1996; Stewart,, 1997, 1977; Chang et al.,, 1997), but none of them provides a linear perturbation expansion with a max-norm control on the remainder term, which is necessary for our development of the error bound on the Karcher mean and the subsequent applications in distributed PCA.

Now we are ready to present a deterministic characterization of the Karcher mean. For convenience, we write $\bm{N}=(\bm{R}^{\top}\ \bm{B}^{\top})^{\top}$ and $\bm{E}^{m}=(\bm{E}^{1,m^{\top}}\ \bm{E}^{2,m^{\top}})^{\top}$ such that $\bm{R},\bm{E}^{1,m}\in\mathbb{R}^{K\times K}$ and $\bm{B},\bm{E}^{2,m}\in\mathbb{R}^{(p-K)\times K}$. Note that $\bm{R}$ is a lower triangular matrix with positive diagonal elements since $\bm{N}\in\mathcal{L}^{*}(p,K)$. In the following theorem, we will show that when $\epsilon_{0}=\max_{m}\|\bm{E}^{m}\|_{\max}$ is sufficiently small, the reduced Cholesky factor $\widetilde{\bm{N}}$ of $\widetilde{\bm{A}}$ differs from $\bm{N}$ by a term linear in $\frac{1}{M}\sum_{m=1}^{M}\bm{E}^{m}$ and an extra term of order $\mathcal{O}_{\max}(\epsilon_{0}^{2})$. Recall that $S^{*}(p,K)$ is the manifold of restricted PSD matrices.

Theorem 3.4** (Karcher Mean on $S^{*}(p,K)$).**

When $\epsilon_{0}=\max_{m}\|\bm{E}^{m}\|_{\max}$ is sufficiently small, the reduced Cholesky factor $\widetilde{\bm{N}}$ of the Karcher mean $\widetilde{\bm{A}}$ of $\{\bm{A}^{m}=(\bm{N}+\bm{E}^{m})(\bm{N}+\bm{E}^{m})^{\top}\}_{m=1}^{M}$ on $S^{*}(p,K)$ is

[TABLE]

where $f_{\bm{R}}(\cdot)$ is given in Lemma 3.3.

From Theorem 3.4, one may easily derive a deterministic upper bound on $\|\widetilde{N}-N\|_{\textrm{F}}$ using the triangular inequality, which depends on $\|\frac{1}{M}\sum_{m=1}^{M}\bm{E}^{m}\|_{\textrm{F}}$ and $pK\epsilon_{0}^{2}$.

Corollary 3.5**.**

Under the same conditions of Theorem 3.4, if $\|\bm{N}\|_{2}\leq C$ and $\|\bm{R}^{-1}\|\leq C$ for some constant $C>0$, then we have

[TABLE]

In applications such as distributed PCA, the Frobenius norm of the average $\frac{1}{M}\sum_{m=1}^{M}\bm{E}^{m}$ is much smaller than that of $\bm{E}^{m}$. Thus, by Corollary 3.5, the Karcher mean $\widetilde{\bm{N}}$ is a better approximation of $\bm{N}$ than any $\bm{N}^{m}$ (the reduced Cholesky factor of $\bm{A}^{m}$).

4 Applications to Distributed PCA

This section applies Theorem 3.4 to show that the distributed PCA algorithm, LRC-dPCA proposed by Neuman et al., (2021), achieves the same performance as the full sample PCA when the local sample size is sufficiently large.

4.1 Distributed PCA and LRC-dPCA

We start with the distributed PCA setting as well as the LRC-dPCA algorithm. For simplicity, we consider a balanced setting, in which we have $M$ machines and the $m$-th machine has $n$ samples $\{\bm{x}_{i}^{m}\}_{i=1}^{n}\subset\mathbb{R}^{p}$. Denote by $N=Mn$ the total number of samples. Assume all samples are i.i.d. sub-Gaussian with mean ${\bm{0}}$ and covariance $\bm{\Sigma}$.

Definition 4.1** (sub-Gaussian).**

We say a random vector $\bm{x}\in\mathbb{R}^{p}$ is sub-Gaussian with mean ${\bm{0}}$ and covariance $\bm{\Sigma}$ if $\bm{z}=\bm{\Sigma}^{-1/2}\bm{x}$ is sub-Gaussian with mean ${\bm{0}}$ and covariance $\bm{I}_{p}$, i.e., there exists a constant $\sigma>0$ such that the following inequality holds,

[TABLE]

Remark 4.2**.**

Fan et al., (2019)** and Bhaskara and Wijewardena, (2019) use the following equivalent definition of a sub-Gaussian vector: $\bm{x}\in\mathbb{R}^{d}$ is sub-Gaussian with mean ${\bm{0}}$ and covariance $\bm{\Sigma}$ if there exists a constant $C>0$ such that $\|\bm{u}^{\top}\bm{x}\|_{\psi_{2}}\leq C\sqrt{\mathbb{E}(\bm{u}^{\top}\bm{x})^{2}},\forall\bm{u}\in\mathbb{R}^{d}$. For more information on the equivalent definitions of sub-Gaussian vectors, one may refer to Vershynin, (2012).

Given a positive integer $K$, the goal is to compute the top $K$ eigenspace of $\bm{\Sigma}$ using all data on $M$ machines with small communication cost. We consider the LRC-dPCA algorithm proposed by Neuman et al., (2021), collected in Algorithm 1. Following this algorithm, we first compute $\widehat{\bm{\Sigma}}^{m}=\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}^{m}\bm{x}_{i}^{m\top}$ on each local machine and then compute the top $K$ eigenvectors $\widehat{\bm{V}}^{m}=(\widehat{\bm{v}}_{1}^{m},\ldots,\widehat{\bm{v}}_{K}^{m})\in\mathcal{O}_{p\times K}$ and eigenvalues $\widehat{\bm{\Lambda}}^{m}={\rm diag}(\lambda_{1}^{m},\ldots,\lambda^{m})$ of $\widehat{\bm{\Sigma}}^{m}$. After communicating these local estimators $\widehat{\bm{V}}^{m},\widehat{\bm{\Lambda}}^{m}$ to a central server, we compute the Karcher mean $\widetilde{\bm{A}}$ of $\{\widehat{\bm{V}}^{m}(\bm{\Lambda}^{m})^{2}\widehat{\bm{V}}^{m\top}\}_{m=1}^{M}$ on the manifold of restricted PSD matrices111Here we choose $(\widehat{\bm{\Lambda}}^{m})^{2}$ rather than $\widehat{\bm{\Lambda}}^{m}$ only for technical reasons in the theoretical proofs.. Finally, the top $K$ eigenvectors $\widetilde{\bm{V}}\in\mathcal{O}_{p\times K}$ of $\widetilde{\bm{A}}$ is returned.

4.2 Theoretical analysis

A statistical analysis of the LRC-dPCA algorithm is missing in its original paper (Neuman et al.,, 2021). In this subsection, we will utilize our deterministic characterization of the Karcher mean on $S^{*}(p,K)$, i.e., Theorem 3.4, to show that given sufficiently large sub-sample size, LRC-dPCA matches the performance of the full sample PCA. Denote by $\bm{V}=(\bm{v}_{1},\ldots,\bm{v}_{K})\in\mathcal{O}_{p\times K}$ and $\bm{\Lambda}={\rm diag}(\lambda_{1},\ldots,\lambda_{K})$ the top $K$ eigenvectors and eigenvalues of $\bm{\Sigma}$, respectively. To ensure the uniqueness of $\textnormal{span}(\bm{V})$, we assume $\Delta_{K}=\lambda_{K}(\bm{\Sigma})-\lambda_{K+1}(\bm{\Sigma})>0$. Write $\bm{A}=\bm{V}\bm{\Lambda}^{2}\bm{A}^{\top}$ and assume the first $K$ columns of $\bm{A}$ are linearly independent, i.e., $\bm{A}\in S^{*}(p,K)$. When the subsample size $n$ is sufficiently large, we will show that $\widehat{\bm{A}}^{m}=\widehat{\bm{V}}^{m}(\widehat{\bm{\Lambda}}^{m})^{2}\widehat{\bm{V}}^{m\top}$ also belongs to $S^{*}(p,K)$ with high probability. Here $\widehat{\bm{V}}^{m}$ and $\widehat{\bm{\Lambda}}^{m}$ denote the top $K$ eigenvectors and eigenvalues of $\widehat{\bm{\Sigma}}^{m}$ respectively. Denote by $\widetilde{\bm{A}}$ the Karcher mean of $\{\widehat{\bm{A}}^{m}\}_{m=1}^{M}$ on $S^{*}(p,K)$. We further denote by $\bm{N},\widehat{\bm{N}}^{m},\widetilde{\bm{N}}$ the reduced Cholesky factors of $\bm{A},\widehat{\bm{A}}^{m},\widetilde{\bm{A}}$. In addition, we define $\bm{Q}^{*}\in\mathcal{O}_{K\times K}$ by the equality $\bm{N}=\bm{V}\bm{\Lambda}\bm{Q}^{*}$.

In the rest of this subsection, we will apply Theorem 3.4 to study the properties of $\widetilde{\bm{A}}$. First, we show that $\{\widehat{\bm{A}}^{m}\}_{m=1}^{M}$ follow the general signal-plus-noise model (3.4).

Lemma 4.3**.**

Let $\widehat{\bm{E}}^{m}=\widehat{\bm{\Sigma}}^{m}\widehat{\bm{V}}^{m}\widehat{\bm{H}}^{m}\bm{Q}^{*}-\bm{\Sigma}\bm{V}\bm{Q}^{*}$, where $\bm{H}^{m}=\widehat{\bm{V}}^{m\top}\bm{V}$ and $\widehat{\bm{H}}^{m}=\textnormal{sgn}(\bm{H}^{m})\overset{\rm def}{=}\bm{U}_{1}\bm{U}_{2}^{\top}$ with $\bm{U}_{1},\bm{U}_{2}$ given by the singular value decomposition $\bm{H}^{m}=\bm{U}_{1}\bm{\Gamma}\bm{U}_{2}^{\top}$ of $\bm{H}^{m}$. Then $\widehat{\bm{A}}^{m}=(\bm{N}+\widehat{\bm{E}}^{m})(\bm{N}+\widehat{\bm{E}}^{m})^{\top}$.

Let us make several remarks on $\widehat{\bm{E}}^{m}$. It is well-known that

[TABLE]

and thus $\widehat{\bm{V}}^{m}\widehat{\bm{H}}^{m}$ is a good estimator of $\bm{V}$ (Chen et al.,, 2020). Furthermore, by Lemma H.2, when $\epsilon=\max_{m}\|\mathcal{E}^{m}\|_{2}/\Delta_{K}\leq 1/10$ with $\mathcal{E}^{m}=\widehat{\bm{\Sigma}}^{m}-\bm{\Sigma}$, $\widehat{\bm{V}}^{m}\widehat{\bm{H}}^{m}$ has the following first-order expansion around $\bm{V}$,

[TABLE]

where $g$ is a linear function defined in Lemma H.2. Substituting (4.1) into the definition of $\widehat{\bm{E}}^{m}$, we obtain the following linear expansion of $\widehat{\bm{E}}^{m}$ in terms of $\mathcal{E}^{m}$,

[TABLE]

Since $g$ is linear in its argument, the leading term of $\frac{1}{M}\sum_{m=1}^{M}\widehat{\bm{E}}^{m}$ is linear in $\frac{1}{M}\sum_{m=1}^{M}\mathcal{E}^{m}$. This enables an upper bound for $\|\frac{1}{M}\sum_{m=1}^{M}\widehat{\bm{E}}^{m}\|_{\textrm{F}}$, provided by the following lemma.

Lemma 4.4** (Bounding $\|M^{-1}\sum_{m=1}^{M}\widehat{\bm{E}}^{m}\|_{\textrm{F}}$).**

Suppose $\Delta_{K}>0$ and $\|\bm{\Sigma}\|_{2}$ is bounded. Let $\mathcal{E}^{m}=\widehat{\bm{\Sigma}}^{m}-\bm{\Sigma}$ and $\epsilon=\max_{m}\|\mathcal{E}^{m}\|_{2}/\Delta_{K}$. When $\epsilon\leq 1/10$, the following bound

[TABLE]

holds for some constant $C>0$.

To apply Theorem 3.4, we also need to upper bound the max norm $\epsilon_{0}=\max_{m}\|\widehat{\bm{E}}^{m}\|_{\max}$. Again, this is based on the first-order expansion (4.2) of $\bm{E}^{m}$.

Lemma 4.5** (Bounding $\max_{m}\|\widehat{\bm{E}}^{m}\|_{\max}$).**

Assume $\Delta_{K}>0$ and $\|\bm{\Sigma}\|_{2}$ is bounded. When $\epsilon=\max_{m}\|\mathcal{E}^{m}\|_{2}/\Delta_{K}\leq 1/10$, we have with probability at least $1-2Me^{-C_{1}n\delta_{1}^{2}}-Me^{-C_{2}\sqrt{\delta_{2}n/r}}$ that

[TABLE]

for some constants $C_{1},C_{2},C_{3}>0$ and $r=\textnormal{Tr}(\bm{\Sigma})/\lambda_{1}(\bm{\Sigma})$. In addition, when $n\gtrsim\log^{3}(pM)r^{2}$, we have with probability at least $1-2p^{-1}$ that

[TABLE]

for some constant $C>0$.

In Lemma 4.5, we show that $\|\widehat{\bm{E}}^{m}\|_{\max}\lesssim\sqrt{\log(p)/n}$ with high probability when $n\gtrsim\log^{3}(p)r^{2}$. By (4.2) and Lemma H.1, we can show that $\|\widehat{\bm{E}}^{m}\|_{\textrm{F}}\lesssim\sqrt{p/n}$ with high probability. The upper bound on $\|\widehat{\bm{E}}^{m}\|_{\max}$ is thus smaller by a factor of $\sqrt{p/\log(p)}$ than the upper bound on $\|\widehat{\bm{E}}^{m}\|_{F}$. This implies that $\widehat{\bm{E}}^{m}$ is delocalized across the entries. Moreover, Lemma 4.5 implies that when we apply Theorem 3.4 to the LRC-dPCA algorithm, the remainder term $\mathcal{O}_{\max}(\epsilon_{0}^{2})$ is negligible compared to the leading term $\frac{1}{M}\sum_{m=1}^{M}\widehat{\bm{E}}^{m}-\bm{N}f_{\bm{R}}(\frac{1}{M}\sum_{m=1}^{M}\widehat{\bm{E}}^{1,m})$. This provides the last key ingredient to the following theorem, which gives an upper bound for $\|\widetilde{\bm{N}}-\bm{N}\|_{F}$. Here $\widetilde{\bm{N}}$ is the reduced Cholesky factor of the Karcher mean $\widetilde{\bm{A}}$.

Theorem 4.6** (Bounding $\|\widetilde{\bm{N}}-\bm{N}\|_{\textrm{F}}$).**

Assume $\Delta_{K}>0$ and $\|\bm{\Sigma}\|_{2}$ is bounded. Partition $\bm{N}=(\bm{R}^{\top}\ \bm{B}^{\top})^{\top}$ such that $\bm{R}\in\mathbb{R}^{K\times K}$ and $\bm{B}\in\mathbb{R}^{(p-K)\times K}$ and assume $\|\bm{R}^{-1}\|_{2}\leq C$ for some constant $C>0$. When $\epsilon=\max_{m}\|\mathcal{E}^{m}\|_{2}/\Delta_{K}\leq 1/10$ and $\epsilon_{0}=\max_{m}\|\widehat{\bm{E}}^{m}\|_{\max}$ is sufficiently small, the following bound

[TABLE]

holds. Define $r=\textnormal{Tr}(\bm{\Sigma})/\lambda_{1}(\bm{\Sigma})$, $\widetilde{r}_{1}=(\log^{2}(pM)r)\vee(\log(pM)\sqrt{p})$ and $\widetilde{r}_{2}=\sqrt{p}\log^{4}(pM)r^{2}$. Then we have with probability at least $1-4p^{-1}$ that

[TABLE]

When $n\gtrsim\widetilde{r}_{2}/\widetilde{r}_{1}$, the third term is negligible. When we further assume $n\gtrsim M\widetilde{r}_{1}^{2}/(\log^{2}(p)r)$, the upper bound reduces to

[TABLE]

Theorem 4.6 shows that given sufficiently large local sample size, i.e., $n\gtrsim M\widetilde{r}_{1}^{2}/(\log^{2}(p)r)$, $\widetilde{\bm{N}}$ is as good as the full sample estimator of $\bm{N}$ in terms of the Frobenius norm. Moreover, $\|\widetilde{\bm{N}}-\bm{N}\|_{\textrm{F}}$ is of the same order as $\|M^{-1}\sum_{m=1}^{M}\mathcal{E}^{m}\|_{2}$ (see Lemma H.1). Note that the singular vectors of $\bm{N}$ are equal to $\bm{V}$, the singular values of $\bm{N}$ are equal to $\bm{\Lambda}$, and LRC-dPCA uses the singular vectors of $\widetilde{\bm{N}}$ as an estimator of $\bm{V}$. Then it follows from Wedin’s sin($\bm{\Theta}$) theorem (Chen et al.,, 2020) that LRC-dPCA and full sample PCA share the same performance in eigenvector estimation.

Remark 4.7**.**

*Similar to Lemma 4.5, we can show that $\|\widehat{\bm{V}}^{m}\widehat{\bm{H}}^{m}-\bm{V}\|_{\max}\lesssim\sqrt{{\log(p)}/{n}}$ with high probability when $n\gtrsim\log^{3}(p)r^{2}$. Compared to the upper bound $\|\widehat{\bm{V}}^{m}\widehat{\bm{H}}^{m}-\bm{V}\|_{F}\lesssim\sqrt{{p}/{n}}$, the max norm bound again implies that the residual matrix $\widehat{\bm{V}}^{m}\widehat{\bm{H}}^{m}-\bm{V}$ does not concentrate on a few coordinates. This has connections to the infinity norm eigenvector perturbation theory (Fan et al.,, 2018; Chen et al.,, 2020; Abbe et al.,, 2020; Damle and Sun,, 2020; Cape et al.,, 2019). However, most applications in their works require incoherence conditions on the eigenvectors. In contrast, we do not require such conditions. *

4.3 Manifold selection

As one may notice, $\bm{A}=\bm{V}\bm{\Lambda}^{2}\bm{V}^{\top}$ may not belong to $S^{*}(p,K)$, i.e., the first $K$ columns of $\bm{A}$ may be linearly dependent. If we decompose $\bm{A}=\bm{F}\bm{F}^{\top}$ for some $\bm{F}\in\mathbb{R}^{p\times K}$ and write $\bm{F}=(\bm{F}_{1}^{\top}\ \bm{F}_{2}^{\top})^{\top}$ with $\bm{F}_{1}\in\mathbb{R}^{K\times K}$ and $\bm{F}_{2}\in\mathbb{R}^{(p-K)\times K}$. Then the smallest singular value $\sigma_{\min}(\bm{F}_{1})$ of $\bm{F}_{1}$ may be zero or very small depending on $p$. In these cases, the condition $\|\bm{R}^{-1}\|_{2}\leq C$ for some constant $C>0$ in Theorem 4.6 may not hold, and it is not suitable to directly use the manifold $S^{*}(p,K)$ in the LRC-dPCA algorithm.

To fix this issue, we will utilize $\frac{p!}{(p-K)!}$ cousins of the manifold $S^{*}(p,K)$, or equivalently $\mathcal{L}^{*}(p,K)$. Let us introduce these cousin manifolds first. Recall that $\mathcal{L}^{*}(p,K)$ consists of $\bm{N}\in\mathbb{R}^{p\times K}$ such that $\bm{N}_{1:K,1:K}$ is a lower triangular matrix with positive diagonal elements. Here $\bm{N}_{1:K,1:K}\in\mathbb{R}^{K\times K}$ represents the sub-matrix of $\bm{N}$ with row index $[1,\ldots,K]$ and column index $[1,\ldots,K]$. Let $\mathcal{I}=[i_{1},\ldots,i_{K}]$ be an ordered index set of size $K$. A cousin $\mathcal{L}^{*}_{\mathcal{I}}(p,K)$ of $\mathcal{L}^{*}(p,K)$ consists of $\bm{N}\in\mathbb{R}^{p\times K}$ such that $\bm{N}_{\mathcal{I},1:K}$ is a lower triangular matrix with positive diagonal elements. Similarly, we define $S^{*}_{\mathcal{I}}(p,K)$ as the set of all matrices in $S(p,K)$ with the $\mathcal{I}$-th rows linearly independent. Similar to the relationship between $S^{*}(p,K)$ and $\mathcal{L}^{*}(p,K)$, for any $\bm{A}\in S^{*}_{\mathcal{I}}(p,K)$, there exists a unique element $\bm{N}\in\mathcal{L}^{*}_{\mathcal{I}}(p,K)$ such that $\bm{A}=\bm{N}\bm{N}^{\top}$. Also, we define the Riemannian structure on $S^{*}_{\mathcal{I}}(p,K)$ and $\mathcal{L}^{*}_{\mathcal{I}}(p,K)$ in a way similar to (2.1) and (2.2). In addition, all theory established in Section 3 and 4 can be rephrased in the language of $S^{*}_{\mathcal{I}}(p,K)$. The only difference is that the row index set $[1,\dots,K]$ is replaced by $\mathcal{I}$.

Now we are in a position to solve the challenge raised at the beginning of this subsection. If $\bm{A}=\bm{V}\bm{\Lambda}^{2}\bm{V}^{\top}$ does not belong to $S^{*}(p,K)$, then we should choose a suitable ordered index set $\mathcal{I}$ rather than $[1,\ldots,K]$, and then apply the LRC-dPCA algorithm on the manifold $S^{*}_{\mathcal{I}}(p,K)$. Motivated by the condition $\|\bm{R}^{-1}\|\leq C$ in Theorem 4.6, we propose the find_index method in Algorithm 2. Given $\bm{V}$, $\bm{\Lambda}$, and $K$, the algorithm outputs an ordered index set $\mathcal{I}$ of size $K$. To avoid exhaustive search, the algorithm determines $\mathcal{I}$ in a sequential manner. In the $k$th step, we choose an index $i\in[p]$ such that the $k$-by-$k$ matrix $\bm{T}_{k}=\bm{T}[c(\mathcal{I}[1:(k-1)],i),c(1:k)]$ has the largest $\sigma_{k}(\bm{T}_{k})$ among all $p$ candidates, where $c(\cdot)$ indicates the index set. In practice when $\bm{V}$ and $\bm{\Lambda}$ is unknown, we can use $\widehat{\bm{V}}^{1}$ and $\widehat{\bm{\Lambda}}^{1}$ to find a suitable index set and this index set is then shared by all machines.

5 Numerical Experiments

In this section, we present numerical experiments on three synthetic examples: averaging PSD matrices under the intrinsic model, the distributed PCA problems, and averaging PSD matrices under an extrinsic model.

5.1 Averaging PSD matrices

Our first experiment is to illustrate the concentration of the Karcher mean (2.3) under the intrinsic model (3.1), i.e., Theorem 3.1. We set $K=5,\sigma^{2}=1$ and let $p$ vary across $[100,200,300,400]$ and let $M$ range from 30 to 270 with an increment of 30. For each $p$, we generate a $p\times p$ matrix $\bm{\Sigma}$ with elements i.i.d. $\mathcal{N}(0,1)$, and then take $\bm{A}=\bm{V}\bm{\Lambda}\bm{V}^{\top}$, where $\bm{V}=(\bm{v}_{1},\ldots,\bm{v}_{K})$ and $\bm{\Lambda}=(\lambda_{1},\ldots,\lambda_{K})$ are the top $K$ left singular vectors and singular values of $\bm{\Sigma}$, respectively. Given $M$, we generate $\{\bm{A}^{m}\}_{m=1}^{M}$ from the intrinsic model (3.1). Then the Karcher mean $\widetilde{\bm{A}}$ of $\{\bm{A}^{m}\}_{m=1}^{M}$ is computed and the error $\|\widetilde{\bm{A}}-\bm{A}\|_{F}$ is reported in the top figure in Figure 1. As our theory shows, the estimation error turns smaller as $M$ increases or $p$ decreases.

In addition, we compare the Karcher mean $\widetilde{\bm{A}}$, referred to as LRC, with the usual Euclidean method $\widetilde{A}^{\rm eu}$, which is defined as the best rank-$K$ approximation of $M^{-1}\sum_{m=1}^{M}\bm{A}^{m}$. We take $p=100$ and repeat the above data generation processing. The errors $\|\widetilde{\bm{A}}-\bm{A}\|_{\textrm{F}}$ and $\|\widetilde{\bm{A}}^{\rm eu}-\bm{A}\|_{\textrm{F}}$ are reported in the bottom figure in Figure 1. As displayed in the figure, under the intrinsic model, the geometry-aware method, LRC, outperforms the Euclidean method. This justifies the intuition that for models with specific geometric structures, it is better to take that geometric information into account.

5.2 Distributed PCA

Our second experiment studies the Karcher mean under a general signal-plus-noise model. Specifically, we consider the distributed PCA problems and numerically verify Theorem 4.6, which shows that LRC-dPCA achieves the same performance as full sample PCA (fPCA). In our setting, $p=100$, $K=5$, the population covariance $\bm{\Sigma}$ is generated by $\bm{\Sigma}=\bm{V}\bm{V}^{\top}+0.3\bm{I}_{p}$, where $\bm{V}\in\mathbb{R}^{p\times K}$ with elements $\rm{i.i.d.}$ $\mathcal{N}(0,1)$. We first fix the number of machines $M=50$ and let the sub-sample size $n$ vary across $[500,1000,\ldots,2500]$. On the $m$-th machine, we generate $n$ $\rm{i.i.d.}$ samples $\{\bm{x}_{i}^{m}\}_{i=1}^{n}$ from $\mathcal{N}({\bm{0}},\bm{\Sigma})$ and compute the local sample covariance matrix $\widehat{\bm{\Sigma}}^{m}=\sum_{i=1}^{n}\bm{x}_{i}^{m}\bm{x}_{i}^{m\top}/n$. Then we apply four methods , namely fPCA, LRC-dPCA, dPCA-Fan (Fan et al.,, 2019), and dPCA-BW (Bhaskara and Wijewardena,, 2019), to compute the top $K$ eigenvectors of $\bm{\Sigma}$. Let $\widehat{\bm{V}}\in\mathcal{O}_{p\times K}$ be the estimated top $K$ eigenvectors. The error is defined as $\|\widehat{\bm{V}}\widehat{\bm{V}}^{\top}-\bm{V}(\bm{V}^{\top}\bm{V})^{-1}\bm{V}^{\top}\|_{F}$, which is the distance between the population projection matrix $\bm{V}(\bm{V}^{\top}\bm{V})^{-1}\bm{V}^{\top}$ and the estimated projection matrix $\widehat{\bm{V}}\widehat{\bm{V}}^{\top}$. For each $n$ and each method, the experiment is repeated 100 times and the average of error is recorded. The top figure in Figure 2 displays the relationship between $\log(\rm{error})$ and $\log(n)$ for all methods. It turns out that all methods share similar performance and there is a linear relationship between $\log(\rm{error})$ and $\log(n)$ with slope $-1/2$, which verifies the relationship $\rm{error}\sim n^{-1/2}$. Next, we fix the sub-sample size $n=1000$ and let $M$ vary across $[50,100,\ldots,200]$ and repeat the above procedures. The relationship between $\log(\rm{error})$ and $\log(M)$ is reported in the bottom figure of Figure 2. As it displayed, all four methods are almost the same and there is also a linear relationship between $\log(\rm{error})$ and $\log(M)$ with slope $-1/2$, which indicates $\rm{error}\sim M^{-1/2}$. Since there is no specific geometric information in the setting, it is expected that LRC-dPCA only matches (rather than surpasses) the performance of the state-of-the-art methods, dPCA-Fan, dPCA-BW, and the optimal method, fPCA.

5.3 Averaging PSD matrices (extrinsic)

Our third experiment considers another signal-plus-noise model, which adds extrinsic noises to the intrinsic model. Specifically, we set $p=100$, $K=5$, and we generate a $p\times p$ matrix $\bm{\Sigma}$ with elements i.i.d. $\mathcal{N}(0,1)$, and then take $\bm{A}=\bm{V}\bm{\Lambda}\bm{V}^{\top}$, where $\bm{V}=(\bm{v}_{1},\ldots,\bm{v}_{K})$ and $\bm{\Lambda}=(\lambda_{1},\ldots,\lambda_{K})$ are the top $K$ left singular vectors and singular values of $\bm{\Sigma}$, respectively. Given $M$ and $\sigma^{2}$, we generate $\{\bm{A}^{m}\}_{m=1}^{M}$ from the intrinsic model (3.1). Then we add extrinsic noises to $\bm{A}^{m}$ as follows. For each $m$, we generate $\{\bm{x}_{i}^{m}\}_{i=1}^{2000}$ $\rm{i.i.d.}$ from $\mathcal{N}({\bm{0}},\bm{A}^{m}+0.01\bm{I}_{p})$, compute $\widehat{\bm{\Sigma}}^{m}=\sum_{i=1}^{2000}\bm{x}_{i}^{m}\bm{x}_{i}^{m\top}/2000$, and set $\bm{A}^{\prime m}$ as the best rank-$K$ approximation of $\widehat{\bm{\Sigma}}^{m}$. We compute the Karcher mean $\widetilde{\bm{A}}$ of $\{\bm{A}^{\prime m}\}_{m=1}^{M}$, which is referred to as LRC, and report the error $\|\widetilde{\bm{A}}-\bm{A}\|_{\textrm{F}}$. In contrast, we also apply the Euclidean method, which computes the best rank-$K$ approximation $\widetilde{\bm{A}}^{\rm eu}$ of $\sum_{m=1}^{M}\bm{A}^{\prime m}/M$, and report the error $\|\widetilde{\bm{A}}^{\rm eu}-\bm{A}\|_{\textrm{F}}$. First, we set $\sigma^{2}=0.5$ and let $M$ vary across $[100,200,\ldots,1000]$. The errors of both methods are displayed in the top figure of Figure 3. As shown in the figure, the geometry-aware method, LRC, still outperforms the Euclidean method even if extrinsic noises are added to the intrinsic model. Next, we fix $M=400$ and let $\sigma^{2}$ range from $[0,0.1\ldots,0.7]$. The errors of both methods are shown in the bottom figure of Figure 3. Recall that $\sigma^{2}$ denotes the strength of intrinsic noises. The bottom figure indicates that when the intrinsic noises are small, then the geometry-aware method and the Euclidean method are comparable, but when the intrinsic noises becomes large, the geometry-aware method tends to outperform the Euclidean method. Overall, this experiment shows that, in a general signal-plus-noise model, if there exist large intrinsic noises, then it is better to utilize the geometry-aware method.

6 Concluding Remarks

This paper considers the geometry of restricted PSD matrices proposed by Neuman et al., (2021). In particular, we provide a non-asymptotic statistical analysis of the Karcher mean of restricted PSD matrices under an intrinsic model. Moreover, for general signal-plus-noise models, we establish a deterministic error bound concerning the Karcher mean. This is based on a linear perturbation expansion of the QR decomposition, which may be of independent interest. As an application, we use the deterministic error analysis of the Karcher mean to prove that the distributed PCA algorithm, LRC-dPCA, achieves the same performance as the full sample PCA. Motivated by the established theory, we propose a manifold selection procedure for the LRC-dPCA algorithm. Finally, we carry out three synthetic numerical experiments to verify our theories. One observation in the experiment is that if data model has certain geometric structure, then it is better to utilize the geometry-aware method.

Several interesting topics are worth of future studies. In manifold-valued data analysis (Patrangenaru and Ellingson,, 2016), it remains to determine which statistical model is more suitable for the given data. For example, the highly anisotropic diffusion tensor images are modelled as PSD matrices (Bonnabel et al.,, 2013), so it is interesting to investigate the performances of the proposed intrinsic model for such data. Second, it is interesting to extend our study to regression, classification, and clustering problems.

APPENDIX

Appendix A Proof of Theorem 3.1

Proof of Theorem 3.1.

Recall that the Karcher mean $\widetilde{\bm{A}}$ of $\{\bm{A}^{m}\}_{m=1}^{M}$ under the intrinsic model is given by (3.1). First, we give an upper bound on the Frobenius norm of $\frac{1}{M}\sum_{m=1}^{M}\bm{E}^{m}$. By the intrinsic model, we know $\frac{1}{M}\sum_{m=1}^{M}\bm{E}^{m}$ is a mock lower triangular matrix with lower triangular elements i.i.d. $\mathcal{N}(0,\sigma^{2}/M)$. Therefore, by the concentration of $\chi^{2}$ ((2.19) in Wainwright, (2019)), for all $t\in(0,1)$, we have

[TABLE]

with probability at least $1-e^{-pKt^{2}/8}$. In a similar spirit, using union bound, we have for $t\in(0,1)$,

[TABLE]

with probability at least $1-Ke^{-t^{2}/8}$. Here $\bm{E}_{ii}^{m}$ is the $(i,i)$-th element of $\bm{E}^{m}$. Thus with high probability, $\max_{i=1,\ldots,K}|\frac{1}{M}\sum_{m=1}^{M}\bm{E}^{m}_{ii}|^{2}\leq 1/2$ and

[TABLE]

where we use the inequality $|\exp(x)-1|\leq 2|x|$ for $x\leq 1/2$. Denote by $\bm{N}$ the reduced Cholesky factor of $\bm{A}$. Then it holds that $\|\bm{N}\|_{2}=\|\bm{A}\|_{2}^{1/2}\leq C^{1/2}$ for some constant $C>0$. Furthermore, by (3.1), we have for some constants $c_{1},c_{2}>0$ that

[TABLE]

with probability at least $1-e^{-c_{1}pK}$. ∎

Appendix B Proof of Lemma 3.3

Proof of Lemma 3.3.

The proof of this lemma is split up into three steps. First, we assume $\bm{Q}=\bm{I}_{K}$ and show that $\widecheck{\bm{Q}}$ has the form of $\bm{I}_{K}+\widecheck{\bm{P}}+\mathcal{O}_{\max}(\epsilon_{0}^{2})$, where $\widecheck{\bm{P}}\in\mathbb{R}^{K\times K}$ is a skew-symmetric matrix of order $\mathcal{O}_{\max}(\epsilon_{0})$. Second, by taking upper triangular off-diagonal elements of $(\bm{R}+\bm{E})\widecheck{\bm{Q}}^{\top}$ as zero, we derive a closed-form expression of $\widecheck{\bm{P}}$ (up to a higher-order term). Motivated by this closed-form expression, we define a function $f_{\bm{R}}:\mathbb{R}^{K\times K}\mapsto\mathbb{R}^{K\times K}$ satisfying several desired conditions. For example, we have $\widecheck{\bm{Q}}=\bm{I}_{K}+f_{\bm{R}}(\bm{E})+\mathcal{O}_{\max}(\epsilon_{0}^{2})$ and $f_{\bm{R}}$ is linear in its argument. Third, we extend the results to the general case when $\bm{Q}\in\mathcal{O}_{K\times K}$ may differ from $\bm{I}_{K}$.

Step 1. When $\epsilon_{0}=\|\bm{E}\|_{\max}$ is sufficiently small, the matrix $\bm{R}+\bm{E}$ is still non-singular and by QR decomposition there exists a unique orthogonal matrix $\widecheck{\bm{Q}}\in\mathcal{O}_{K\times K}$ such that $\widecheck{\bm{R}}=(\bm{R}+\bm{E})\widecheck{\bm{Q}}^{\top}$ is a lower triangular matrix with positive diagonal elements. In this step, we will show that $\widecheck{\bm{Q}}$ has a form of $\bm{I}_{K}+\widecheck{\bm{P}}+\mathcal{O}_{\max}(\epsilon_{0}^{2})$, where $\widecheck{\bm{P}}\in\mathbb{R}^{K\times K}$ is a skew-symmetric matrix of order $\mathcal{O}_{\max}(\epsilon_{0})$. To that end, we construct $\widecheck{\bm{Q}}$ as a product of $K(K-1)/2$ rotation matrices $\{\widecheck{\bm{Q}}^{ij},1\leq i<j\leq K\}$, which set the upper triangular off-diagonal elements as zero in a sequential fashion. In specific, we arrange these $K(K-1)/2$ rotation matrices in a prescribed order, i.e., $\{(1,2),\ldots,(1,K),(2,3),\ldots,(K-1,K)\}$. In this way, we may relabel $\{\widecheck{\bm{Q}}^{ij},1\leq i<j\leq K\}$ as $\{\widecheck{\bm{Q}}^{(s)},1\leq s\leq K(K-1)/2\}$ and write $\widecheck{\bm{Q}}=\widecheck{\bm{Q}}^{(K(K-1)/2)}\cdots\widecheck{\bm{Q}}^{(1)}$. In the remainder of this proof, we will use $s(i,j)$ to represent the s-index of the $(i,j)$th rotation matrix $\widecheck{\bm{Q}}^{ij}$.

For each $(i,j)$, we set the rotation matrix $\widecheck{\bm{Q}}^{ij}$ as

[TABLE]

where $\theta^{ij}$ is chosen in a sequential fashion such that the $(i,j)$th element of $\widecheck{\bm{R}}^{(s(i,j))}\coloneqq(\bm{R}+\bm{E})\widecheck{\bm{Q}}^{(1)^{\top}}\cdots\widecheck{\bm{Q}}^{(s(i,j))^{\top}}$ is zero and the diagonal elements of $\widecheck{\bm{R}}^{(s(i,j))}$ keep positive. Note that $\widecheck{\bm{Q}}^{ij}$ is by definition an orthogonal matrix. A simple calculation gives that $\theta^{ij}=\theta^{(s(i,j))}=\arctan(\widecheck{\bm{R}}^{(s(i,j)-1)}_{ij}/\widecheck{\bm{R}}^{(s(i,j)-1)}_{ii})$.

Next, we show that $\theta^{ij}$ is a small quantity of order $\mathcal{O}(\epsilon_{0})$ via an deductive argument. First, when $\epsilon_{0}=\|\bm{E}\|_{\max}$ is sufficiently small, $\theta^{(1)}=\theta^{12}=\arctan(\bm{E}_{12}/(\bm{R}_{11}+\bm{E}_{11}))$ is a small quantity of order $\mathcal{O}(\epsilon_{0})$. Thus, by definition of $\widecheck{\bm{Q}}^{12}$, we have $\widecheck{\bm{Q}}^{(1)}=\widecheck{\bm{Q}}^{12}=\bm{I}_{K}+\mathcal{O}_{\max}(\epsilon_{0})$ and

[TABLE]

Note that the error matrix $\bm{E}^{(1)}\coloneqq\widecheck{\bm{R}}^{(1)}-\bm{R}$ is again of order $\mathcal{O}_{\max}(\epsilon_{0})$. This implies that $\theta^{(2)}=\theta^{13}=\arctan(\bm{E}^{(1)}_{13}/(\bm{R}_{11}+\bm{E}^{(1)}_{11}))$ is also a small quantity of order $\mathcal{O}(\epsilon_{0})$. Applying this deductive argument $K(K-1)/2$ times, we conclude that all $\theta^{(s)}$, $1\leq s\leq K(K-1)/2$, are small quantities of order $\mathcal{O}(\epsilon_{0})$.

Now we are able to show that $\widecheck{\bm{Q}}$ has a form of $\bm{I}_{K}+\widecheck{\bm{P}}+\mathcal{O}_{\max}(\epsilon_{0}^{2})$, where $\widecheck{\bm{P}}\in\mathbb{R}^{K\times K}$ is a skew-symmetric matrix of order $\mathcal{O}_{\max}(\epsilon_{0})$. Since $\theta^{ij}$ is of order $\mathcal{O}(\epsilon_{0})$, by Taylor expansion, we have $\sin(\theta^{ij})=\theta^{ij}+\mathcal{O}(\epsilon_{0}^{3})$ and $1-\cos(\theta^{ij})=\mathcal{O}(\epsilon_{0}^{2})$. Thus, we can rewrite $\widecheck{\bm{Q}}^{ij}$ as

[TABLE]

where $\widecheck{\bm{P}}^{ij}_{ij}=-\widecheck{\bm{P}}^{ij}_{ij}=\theta^{ij}$ and $\widecheck{\bm{P}}^{ij}_{kl}=0$ otherwise. Since $\widecheck{\bm{P}}^{ij}$ is of order $\mathcal{O}_{\max}(\epsilon_{0})$, we have

[TABLE]

where $\widecheck{\bm{P}}^{(s(i,j))}=\widecheck{\bm{P}}^{ij}$ and $\widecheck{\bm{P}}=\sum_{s=1}^{K(K-1)/2}\widecheck{\bm{P}}^{(s)}$. Since $\widecheck{\bm{P}}^{(s)}$ is skew-symmetric and of order $\mathcal{O}_{\max}(\epsilon_{0})$ for all $s$, $\widecheck{\bm{P}}$ is also a skew-symmetric matrix of order $\mathcal{O}_{\max}(\epsilon_{0})$, which concludes the proof of step 1.

Step 2. Now we are ready to derive a closed-form expression of $\widecheck{\bm{P}}$ (maybe up to a higher-order term) by taking upper triangular off-diagonal elements of $(\bm{R}+\bm{E})\widecheck{\bm{Q}}^{\top}$ as zero. Substituting (B) into $(\bm{R}+\bm{E})\widecheck{\bm{Q}}^{\top}$, we obtain

[TABLE]

where the second equality follows from the skew-symmetry of $\widecheck{\bm{P}}$. Since $\bm{R}$ is a lower triangular matrix with positive diagonal elements, $\bm{R}$ is invertible and $\bm{R}^{-1}$ is also a lower triangular matrix. As a result, the matrix $\bm{R}^{-1}(\bm{R}+\bm{E})\widecheck{\bm{Q}}^{\top}$ is also a lower triangular matrix. Multiplying LHS and RHS of (B) by $\bm{R}^{-1}$ simultaneously, we obtain

[TABLE]

For convenience, we define a function $\mathcal{U}(\cdot):\mathbb{R}^{K\times K}\mapsto\mathbb{R}^{K\times K},\bm{P}\mapsto\mathcal{U}(\bm{P})$, where $\mathcal{U}(\bm{P})$ takes the upper triangular off-diagonal elements of $\bm{P}$, i.e.,

[TABLE]

Since $\bm{R}^{-1}(\bm{R}+\bm{E})\widecheck{\bm{Q}}^{\top}$ is a lower triangular matrix, we have by (B.3) that

[TABLE]

Since $\widecheck{\bm{P}}$ is skew-symmetric, we get the following closed-form solution of $\widecheck{\bm{P}}$ (up to a higher-order term),

[TABLE]

Motivated by the linear expansion of $\widecheck{\bm{P}}$, we define the following function,

[TABLE]

Note that $f_{\bm{R}}$ is linear in the sense that $f_{\bm{R}}(a\bm{E}+b\bm{F})=af_{\bm{R}}(\bm{E})+bf_{\bm{R}}(\bm{F})$ for all $a,b\in\mathbb{R}$ and $\bm{E},\bm{F}\in\mathbb{R}^{K\times K}$. Also, the image $f_{\bm{R}}(\bm{E})$ is a skew-symmetric matrix, i.e., $(f_{\bm{R}}(\bm{E}))^{\top}=-f_{\bm{R}}(\bm{E})$. Moreover, we have $\|f_{\bm{R}}(\bm{E})\|_{\textrm{F}}\leq\sqrt{2}\|\bm{R}^{-1}\|_{2}\|\bm{E}\|_{\textrm{F}}$. By (B), we can rewrite $\widecheck{\bm{Q}}$ as follows,

[TABLE]

Moreover, by definition of $\widecheck{\bm{R}}$, we have

[TABLE]

Step 3. In general, when $\bm{Q}\in\mathcal{O}_{K\times K}$ may differ from $\bm{I}_{K}$, we can transform the QR decomposition $\widecheck{\bm{R}}\widecheck{\bm{Q}}=\bm{R}\bm{Q}+\bm{E}$ suitably and apply the results in the previous two steps to prove the lemma. In specific, we have

[TABLE]

When $\epsilon_{0}=\|\bm{E}\|_{\max}$ is sufficiently small, $\|\bm{E}\bm{Q}^{\top}\|_{\max}\leq\sqrt{K}\epsilon_{0}$ can also be sufficiently small. In addition, $\widecheck{\bm{Q}}\bm{Q}^{\top}$ is still an orthogonal matrix that appears in the QR decomposition of $\bm{R}+\bm{E}\bm{Q}^{\top}$. Therefore, by (B.4), we have

[TABLE]

By multiplying both LHS and RHS of this equation by $\bm{Q}$, we obtain that

[TABLE]

In addition, by definition of $\widecheck{\bm{R}}$, we have

[TABLE]

which concludes the proof. ∎

Appendix C Proof of Theorem 3.4

Proof of Theorem 3.4.

First, we use Lemma 3.3 to give a first-order perturbation expansion for the reduced Cholesky factor $\bm{N}^{m}$ of $\bm{A}^{m}=(\bm{N}+\bm{E}^{m})(\bm{N}+\bm{E}^{m})^{\top}$. Define $\bm{Q}^{m}\in\mathcal{O}_{K\times K}$ as an orthogonal matrix such that $\bm{N}^{m}=(\bm{N}+\bm{E}^{m})\bm{Q}^{m\top}$, or equivalently, $(\bm{R}+\bm{E}^{1,m})\bm{Q}^{m\top}$ is a lower triangular matrix with positive diagonal elements. By Lemma 3.3, when $\|\bm{E}^{1,m}\|_{\max}\leq\epsilon_{0}$ is sufficiently small, we have

[TABLE]

where $f_{\bm{R}}$ is defined in Lemma 3.3. By definition of $\bm{Q}^{m}$, we have

[TABLE]

where we use the property $f_{\bm{R}}(\bm{E}^{1,m})^{\top}=-f_{\bm{R}}(\bm{E}^{1,m})$. Using this linear perturbation expansion, we are now able to characterize the Karcher mean $\widetilde{\bm{A}}=\widetilde{\bm{N}}\widetilde{\bm{N}}^{\top}$ (or $\widetilde{\bm{N}}$) of $\{\bm{A}^{m}\}_{m=1}^{M}$ (or $\{\bm{N}^{m}\}_{m=1}^{M}$) on the manifold $S^{*}(p,K)$ (or $\mathcal{L}^{*}(p,K)$). By the LRC algorithm, i.e., (2.3), we have $\widetilde{\bm{N}}$ is equal to $\frac{1}{M}\sum_{m=1}^{M}\bm{N}^{m}$ except that the diagonal elements of $\widetilde{\bm{N}}$ are given by

[TABLE]

However, when $\epsilon_{0}$ is sufficiently small, $|\bm{N}^{m}_{ii}-\bm{N}_{ii}|$ is of order $\mathcal{O}(\epsilon_{0})$ and thus

[TABLE]

Therefore, we have

[TABLE]

where the last equality follows from the linear property of $f_{\bm{R}}(\cdot)$. ∎

Appendix D Proof of Lemma 4.3

Proof of Lemma 4.3.

Since $\bm{V},\bm{\Lambda}$ denote the top $K$ eigenvectors and eigenvalues of $\bm{\Sigma}$, respectively, we have $\bm{\Sigma}\bm{V}=\bm{V}\bm{\Lambda}$ and thus $\bm{N}=\bm{V}\bm{\Lambda}\bm{Q}^{*}=\bm{\Sigma}\bm{V}\bm{Q}^{*}$. Similarly, we have $\widehat{\bm{\Sigma}}^{m}\widehat{\bm{V}}^{m}=\widehat{\bm{V}}^{m}\widehat{\bm{\Lambda}}^{m}$. Since $\widehat{\bm{H}}^{m}$ and $\bm{Q}^{*}$ are both orthogonal matrices, we have

[TABLE]

which concludes our proof. ∎

Appendix E Proof of Lemma 4.4

Proof of Lemma 4.4.

The proof of this lemma is based on a first-order expansion of $\widehat{\bm{E}}^{m}$. Define $\mathcal{E}^{m}=\widehat{\bm{\Sigma}}^{m}-\bm{\Sigma}$ and $\epsilon=\max_{m}\|\mathcal{E}^{m}\|_{2}/\Delta_{K}$. When $\epsilon\leq 1/10$, by Lemma H.2, we have

[TABLE]

where

[TABLE]

with $\bm{G}_{j}=\sum_{i>K}(\lambda_{i}-\lambda_{j})^{-1}\bm{v}_{i}\bm{v}_{i}^{\top}$ for $j\in[K]$ and $\lambda_{i}$/$\bm{v}_{i}$ being the $i$th eigenvalue/eigenvector of $\bm{\Sigma}$. By definition of $\widehat{\bm{E}}^{m}$, we have

[TABLE]

Since $\textnormal{vec}\circ g\circ\textnormal{vec}^{-1}$ is a linear mapping from $\mathbb{R}^{pK}$ to $\mathbb{R}^{pK}$, where $\textnormal{vec}:\mathbb{R}^{p\times K}\mapsto\mathbb{R}^{pK}$ is the vectorization mapping, the average of $\widehat{\bm{E}}^{m}$ can be expressed as

[TABLE]

Thus, by the triangular inequality, we have

[TABLE]

Since $\|\bm{Q}^{*}\|_{2}=1$ and $\|\bm{V}\|_{\textrm{F}}=\sqrt{K}$, we have

[TABLE]

In addition, since $\|\bm{\Sigma}\bm{G}_{j}\|_{2}\leq\Delta_{K}^{-1}\lambda_{K}$ for all $j\in[K]$, we have

[TABLE]

Thus, we have

[TABLE]

for some constant $C>0$. ∎

Appendix F Proof of Lemma 4.5

Proof of Lemma 4.5.

Similar to Lemma 4.4, the proof of this lemma is also based on the first-order expansion of $\widehat{\bm{E}}^{m}$. Define $\mathcal{E}^{m}=\widehat{\bm{\Sigma}}^{m}-\bm{\Sigma}$ and $\epsilon=\max_{m}\|\mathcal{E}^{m}\|_{2}/\Delta_{K}$. When $\epsilon\leq 1/10$, by Lemma H.2, we have

[TABLE]

where

[TABLE]

with $\bm{G}_{j}=\sum_{i>K}(\lambda_{i}-\lambda_{j})^{-1}\bm{v}_{i}\bm{v}_{i}^{\top}$ for $j\in[K]$ and $\lambda_{i}$/$\bm{v}_{i}$ being the $i$th eigenvalue/eigenvector of $\bm{\Sigma}$. By definition of $\widehat{\bm{E}}^{m}$, we have

[TABLE]

By the triangular inequality and the fact that $\|\cdot\|_{\max}\leq\|\cdot\|_{\textrm{F}}$, we have

[TABLE]

where the second inequality follows from the inequality $\|\cdot\bm{Q}^{*}\|_{\max}\leq\sqrt{K}\|\cdot\|_{\max}$. Thus, to bound $\|\widehat{\bm{E}}^{m}\|_{\max}$, it suffices to bound $\|\mathcal{E}^{m}\bm{V}\|_{\max}$, $\|\bm{\Sigma}g(\mathcal{E}^{m}\bm{V})\|_{\max}$, and $\epsilon^{2}=(\max_{m}\|\mathcal{E}^{m}\|_{2}/\Delta_{K})^{2}$ separately. To give an upper bound on the first two terms, we will need Proposition 1 in Pilanci and Wainwright, (2015), which is presented below for reader’s convenience.

Proposition F.1** (Proposition 1 in Pilanci and Wainwright, (2015)).**

Let $\left\{\bm{z}_{i}\right\}_{i=1}^{n}\subset\mathbb{R}^{p}$ be i.i.d. samples generated from a zero-mean sub-Gaussian distribution with $\operatorname*{\rm Cov}(\bm{z}_{i})=\bm{I}_{p}$. Then there exist some universal constants $C_{1},C_{2}>0$ such that for any subset $\mathcal{Y}\subset\mathbb{S}^{p-1}$, we have with probability at least $1-e^{-C_{2}n\delta^{2}}$,

[TABLE]

where $\bm{Z}^{\top}=\left(\bm{z}_{1},\ldots,\bm{z}_{n}\right)\in\mathbb{R}^{p\times n}$ and $\mathbb{W}(\mathcal{Y})$ is the Gaussian width of the subset $\mathcal{Y}$. Specifically, $\mathbb{W}(\mathcal{Y})$ is defined by

[TABLE]

where the expectation is taken on $\bm{h}\in\mathbb{R}^{p}$, which is a standard normal random vector.

1: bound $\|\mathcal{E}^{m}\bm{V}\|_{\max}$. Denote by $\bm{e}_{l}\in\mathbb{R}^{p}$ the basis vector with value 1 at the $l$th entry and 0 at other entries. Then the max norm has the following expression,

[TABLE]

Let $\bm{z}_{i}^{m}=\bm{\Sigma}^{-1/2}\bm{x}_{i}^{m}$ and $\bm{Z}^{m}=(\bm{z}_{1}^{m},\ldots,\bm{z}_{n}^{m})^{\top}$. Since $\{\bm{x}_{i}^{m}\}_{i=1}^{n}$ are i.i.d. sub-Gaussian with mean ${\bm{0}}$ and covariance $\bm{\Sigma}$, $\{\bm{z}_{i}^{m}\}_{i=1}^{n}$ are i.i.d. sub-Gaussian with mean ${\bm{0}}$ and covariance $\bm{I}_{p}$. By definition of $\mathcal{E}^{m}$, we have

[TABLE]

By the polarization equality, we have

[TABLE]

Then by the triangular inequality, we have

[TABLE]

Since $\|\bm{e}_{l}\|_{2}=\|\bm{v}_{k}\|_{2}=1$, we have

[TABLE]

for all $l\in[p]$ and $k\in[K]$. Define $\mathcal{Y}_{1}\subset\mathbb{S}^{p-1}$ as follows,

[TABLE]

then we have

[TABLE]

where $C_{1}>0$ is a constant dependent on $\|\bm{\Sigma}\|_{2}$. We remark here that in the proof notations $C,C_{1},C_{2}$ represent some universal constants, which may vary according to the context. By Proposition F.1, there exist some universal constants $C_{1},C_{2}>0$ such that the following inequality holds

[TABLE]

with probability at least $1-e^{-C_{2}n\delta^{2}}$. Since $\mathcal{Y}_{1}$ is a finite set with cardinality $|\mathcal{Y}_{1}|\leq Cp$ for some constant $C>0$, by the maximal inequality (Mohri et al.,, 2018), the following inquality

[TABLE]

holds for some constant $C_{1}>0$. Thus with probability at least $1-e^{-C_{2}n\delta^{2}}$, we have

[TABLE]

2: bound $\|\bm{\Sigma}g(\mathcal{E}^{m}\bm{V})\|_{\max}$. The proof of this step is similar to that of step one. By definition of $g$, we have

[TABLE]

Then we may write the max norm as

[TABLE]

Similar to step one, by the polarization equality, the triangular inequality, and the fact $\bm{G}_{k}=\bm{G}_{k}^{\top}$, we have

[TABLE]

By definition of $\bm{G}_{k}$, we have

[TABLE]

and thus $\|\bm{\Sigma}^{1/2}\bm{G}_{k}\bm{\Sigma}\|_{2}\leq\lambda_{K}^{3/2}\Delta_{K}^{-1}$. Since $\|\bm{e}_{l}\|_{2}=\|\bm{v}_{k}\|_{2}=1$, we have

[TABLE]

Define the following set $\mathcal{Y}_{2}\subset\mathbb{S}^{p-1}$,

[TABLE]

Then we have

[TABLE]

for some constant $C_{1}>0$ dependent on $\|\bm{\Sigma}\|_{2}$ and $\Delta_{K}$. Again by Proposition F.1, we have with probability at least $1-e^{-C_{2}n\delta^{2}}$ that

[TABLE]

for some universal constants $C_{1},C_{2}>0$. Since $\mathcal{Y}_{2}$ is a finite set with cardinality $|\mathcal{Y}_{2}|\leq Cp$ for some constant $C>0$, by the maximal inequality (Mohri et al.,, 2018), we have

[TABLE]

for some constant $C_{1}>0$. Thus with probability at least $1-e^{-C_{2}n\delta^{2}}$, we have

[TABLE]

3: bound $\epsilon^{2}$. We will use the tail bound of $\|\mathcal{E}^{m}\|_{2}$ in Lemma H.1 to bound $\epsilon^{2}$. By Lemma H.1, we have with probability at least $1-e^{-C\sqrt{\frac{\delta n}{r}}}$ that

[TABLE]

for some constant $C>0$ and $r=\textnormal{Tr}(\bm{\Sigma})/\lambda_{1}\leq p$. Then by union bound, we have with probability at least $1-Me^{-C\sqrt{\frac{\delta n}{r}}}$ that

[TABLE]

for some constant $C>0$.

Last: bound $\max_{m}\|\widehat{\bm{E}}^{m}\|_{\max}$. We will combine the results from 1 to 3 and apply a union bound to obtain the upper bound on $\max_{m}\|\widehat{\bm{E}}^{m}\|_{\max}$. In specific, by union bound, we have with probability at least $1-2Me^{-C_{1}n\delta_{1}^{2}}-Me^{-C_{2}\sqrt{\delta_{2}n/r}}$ that

[TABLE]

for some constants $C_{1},C_{2},C_{3}>0$. Take $\delta_{1}=\sqrt{\frac{\log(2Mp)}{C_{1}n}}$ and $\delta_{2}=\frac{\log(Mp)^{2}r}{C_{2}^{2}n}$, then we have with probability at least $1-2p^{-1}$ that,

[TABLE]

for some constants $C_{1},C_{2}>0$. When $n\gtrsim\log^{3}(pM)r^{2}$, with probability at least $1-2p^{-1}$ the following bound

[TABLE]

holds for some constant $C>0$. ∎

Appendix G Proof of Theorem 4.6

Proof of Theorem 4.6.

We will combine the results in Theorem 3.4, Lemma 4.3, Lemma 4.4, and Lemma 4.5 to prove this theorem. Recall that $\widehat{\bm{E}}^{m}=\widehat{\bm{\Sigma}}^{m}\widehat{\bm{V}}^{m}\widehat{\bm{H}}^{m}\bm{Q}^{*}-\bm{\Sigma}\bm{V}\bm{Q}^{*}$, $\mathcal{E}^{m}=\widehat{\bm{\Sigma}}^{m}-\bm{\Sigma}$, $\epsilon=\max_{m}\|\mathcal{E}^{m}\|_{2}/\Delta_{K}$, and $\epsilon_{0}=\max_{m}\|\widehat{\bm{E}}^{m}\|_{\max}$. In addition, we partition $\bm{N}=(\bm{R}^{\top}\ \bm{B}^{\top})^{\top}$ and $\widehat{\bm{E}}^{m}=(\widehat{\bm{E}}^{1,m^{\top}}\ \widehat{\bm{E}}^{2,m^{\top}})^{\top}$ such that $\bm{R},\widehat{\bm{E}}^{1,m}\in\mathbb{R}^{K\times K}$ and $\bm{B},\widehat{\bm{E}}^{2,m}\in\mathbb{R}^{(p-K)\times K}$. By Lemma 4.3, we have $\widehat{\bm{A}}^{m}=(\bm{N}+\widehat{\bm{E}}^{m})(\bm{N}+\widehat{\bm{E}}^{m})^{\top}$ for all $m\in[M]$, where $\bm{N}$ is the reduced Cholesky factor of $\bm{A}=\bm{V}\bm{\Lambda}^{2}\bm{V}^{\top}$. By Lemma 4.5, $\epsilon_{0}=\max_{m}\|\widehat{\bm{E}}^{m}\|_{\max}$ is sufficiently small with high probability when $n$ is sufficiently large. Thus, we can apply Theorem 3.4 to the LRC-dPCA algorithm to obtain the desired result. In specific, by Theorem 3.4, we have

[TABLE]

where $f_{\bm{R}}$ is defined in Lemma 3.3. By the triangular inequality, we have

[TABLE]

where the second inequality is due to the property $\|f_{\bm{R}}(\cdot)\|_{\textrm{F}}\leq\sqrt{2}\|\bm{R}^{-1}\|_{2}\|\cdot\|_{\textrm{F}}$. By Lemma 4.4, we have

[TABLE]

for some constant $C>0$.

Next, we give a high probability bound on $\|\widetilde{\bm{N}}-\bm{N}\|_{\textrm{F}}$ in three steps. First, by Lemma H.1, we have with probability at least $1-e^{-\frac{\delta_{1}}{C_{1}\lambda_{1}\sqrt{r/(Mn)}}}$ that

[TABLE]

for some constant $C_{1}>0$ and $r=\textnormal{Tr}(\bm{\Sigma})/\lambda_{1}(\bm{\Sigma})$. Second, by Lemma H.1 and the union bound, we have with probability at least $1-Me^{-\frac{\delta_{2}}{C_{2}\lambda_{1}\sqrt{r/n}}}$ that

[TABLE]

for some constant $C_{2}>0$. Third, by Lemma 4.5, we have with probability at least $1-2Me^{-C_{3}n\delta_{3}^{2}}-Me^{-C_{4}\sqrt{\delta_{4}n/r}}$ that

[TABLE]

for some constants $C_{3},C_{4},C_{5}>0$. By the union bound, we combine these three high probability bounds with (G.1) to obtain the desired result. In specific, with probability at least $1-e^{-\frac{\delta_{1}}{C_{1}\lambda_{1}\sqrt{r/(Mn)}}}-Me^{-\frac{\delta_{2}}{C_{2}\lambda_{1}\sqrt{r/n}}}-2Me^{-C_{3}n\delta_{3}^{2}}-Me^{-C_{4}\sqrt{\delta_{4}n/r}}$, the following inequality

[TABLE]

holds for some constants $C_{1},C_{2},C_{3},C_{4}>0$. Take $\delta_{1}=C_{1}\lambda_{1}\log(p)\sqrt{r/(Mn)}$, $\delta_{2}=C_{2}\lambda_{1}\log(pM)\sqrt{r/n}$, $\delta_{3}=\sqrt{\frac{\log(2pM)}{C_{3}n}}$, and $\delta_{4}=\frac{\log^{2}(pM)r}{C_{4}^{2}n}$, then we have with probability at least $1-4p^{-1}$ that

[TABLE]

When $n\gtrsim(\log^{2}(pM)\sqrt{p}r)\vee(\log^{3}(pM)r^{2})$, we have with probability at least $1-4p^{-1}$ that

[TABLE]

In addition, when $n\gtrsim\frac{M((\log^{4}(pM)r^{2})\vee(\log^{2}(pM)p))}{\log^{2}(p)r}$, we have with probability at least $1-4p^{-1}$ that

[TABLE]

which concludes the proof. ∎

Appendix H Auxiliary Lemmas

The following lemma gives a tail bound of $\|\widehat{\bm{\Sigma}}-\bm{\Sigma}\|_{2}$ in the sub-Gaussian case.

Lemma H.1** (Lemma 3 in Fan et al., (2019)).**

Suppose $\{\bm{x}_{i}\}_{i=1}^{n}\subset\mathbb{R}^{p}$ are i.i.d. sub-Gaussian with mean ${\bm{0}}$ and covariance $\bm{\Sigma}$. Let $\widehat{\bm{\Sigma}}=\frac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{\top}$ be the sample covariance matrix, $\{\lambda_{j}\}_{j=1}^{p}$ be the eigenvalues of $\bm{\Sigma}$ sorted in descending order, and $r=\textnormal{Tr}(\bm{\Sigma})/\lambda_{1}$. There exist constants $C_{1}\geq 1$ and $C_{2}\geq 0$ such that when $n\geq r$, we have

[TABLE]

and $\|\|\widehat{\bm{\Sigma}}-\bm{\Sigma}\|_{2}\|_{\psi_{1}}\leq C_{2}\lambda_{1}\sqrt{r/n}$.

The following lemma provides a first-order expansion of $\widehat{\bm{V}}\widehat{\bm{H}}$ around $\bm{V}$, where $\widehat{\bm{H}}=\mathop{\mathrm{argmin}}_{\bm{O}\in\mathcal{O}_{K\times K}}\|\widehat{\bm{V}}\bm{O}-\bm{V}\|_{\textrm{F}}$.

Lemma H.2** (Lemma 8 in Fan et al., (2019)).**

Let $\bm{\Sigma},\widehat{\bm{\Sigma}}\in\mathbb{R}^{p\times p}$ be symmetric matrices with eigenvalues $\{\lambda_{i}\}_{i=1}^{p}$ and $\{\widehat{\lambda}_{i}\}_{i=1}^{p}$ (in descending order) and eigenvectors $\{\bm{v}_{j}\}_{j=1}^{p},\{\widehat{\bm{v}}_{j}\}_{j=1}^{p}$ such that $\bm{\Sigma}\bm{v}_{j}=\lambda_{j}\bm{v}_{j}$ and $\widehat{\bm{\Sigma}}\widehat{\bm{v}}_{j}=\widehat{\lambda}_{j}\widehat{\bm{v}}_{j}$ for $j\in[p]$. Define $\mathcal{E}=\widehat{\bm{\Sigma}}-\bm{\Sigma}$, $S=\{s+1,\ldots,s+K\}$ for some fixed $s\in\{0,1,\ldots,p-K\}$, $\bm{G}_{j}=\sum_{i\notin S}(\lambda_{i}-\lambda_{j})^{-1}\bm{v}_{i}\bm{v}_{i}^{\top}$ for $j\in[K]$, and

[TABLE]

Let $\bm{V}=(\bm{v}_{s+1},\ldots,\bm{v}_{s+K})$, $\widehat{\bm{V}}=(\widehat{\bm{v}}_{s+1},\ldots,\widehat{\bm{v}}_{s+K})$, $\bm{H}=\widehat{\bm{V}}^{\top}\bm{V}$, and $\widehat{\bm{H}}=\textnormal{sgn}(\bm{H})\coloneqq\bm{U}_{1}\bm{U}_{2}^{\top}$, where $\bm{H}=\bm{U}_{1}\bm{\Gamma}\bm{U}_{2}^{\top}$ is the unique singular value decomposition of $\bm{H}$. If $\Delta=\min\{\lambda_{s}-\lambda_{s+1},\lambda_{s+K}-\lambda_{s+K+1}\}>0$ and $\epsilon=\|\mathcal{E}\|_{2}/\Delta\leq 1/10$, where $\lambda_{0}=\infty$ and $\lambda_{p+1}=-\infty$, we have

[TABLE]

It is worth noting that since $\|g(\cdot)\|_{\textrm{F}}\leq\Delta^{-1}\|\cdot\|_{\textrm{F}}$, the Frobenius norm of the remainder term $\widehat{\bm{V}}\widehat{\bm{H}}-\bm{V}-g(\mathcal{E}\bm{V})$ is of order

[TABLE]