A Communication-Efficient And Privacy-Aware Distributed Algorithm For Sparse PCA

We consider the following distributed setting. The data matrix $A$ is divided into $d$ blocks, each containing a number of samples; namely, $A=[A_{1},A_{2},\dotsc,A_{d}]$ where $A_{i}\in\mathbb{R}^{n\times m_{i}}$ so that $m_{1}+\dotsb+m_{d}=m$. This is a natural setting since each sample is a column of $A$. These submatrices $A_{i}$, $i=1,\dotsc,d$, are stored locally in $d$ locations, possibly having been collected at different locations by different agents, and all the agents are connected through a communication network. According to the splitting scheme of $A$, the function $f(Z)$ can also be distributed into $d$ agents, namely,

In terms of network topology, we only need to assume that the network allows global summation operations (say, through all-reduce type of communications [43]) which are required by our algorithm. In particular, our algorithm will operate well in the federated-learning setting [41] where all the agents are connected to a center server so that global summations can be readily achieved in the network. To this extent, we say our algorithm is a centralized algorithm.

For most distributed algorithms, since the amount of communications per iteration remains essentially the same, the total communication overhead is proportional to the required number of iterations regardless of the underlying network topology. Therefore, our goal is to devise an algorithm that converges fast in terms of iteration counts, while considerations on other network topologies and communication patterns are beyond the scope of the current paper.

In certain applications, such as those in healthcare and financial industry [38, 60], preserving privacy of local data is of primary importance. In this paper, we consider the following scenario: each agent $i$ wants to keep its local dataset (i.e., $A_{i}A_{i}^{\top}$) from being discovered by any other agents including the center. In this situation, it is not an option to implement a pre-agreed encryption or a coordinated masking operation. For convenience, we will call such a privacy situation of intrinsic privacy. For an algorithm to preserve intrinsic privacy, publicly exchanged information must be safe so that none of the local-data matrices can be figured out by any means. We will show later that, in general, algorithms based on traditional methods with distributed matrix multiplications are not intrinsically private.

Optimization problems with orthogonality constraints have been actively investigated in recent decades, for which many algorithms and solvers have been developed, such as, gradient approaches [40, 42, 1, 3], conjugate gradient approaches [16, 45, 61], constraint preserving updating schemes [52, 32], Newton methods [30, 29], trust-region methods [2], multipliers correction frameworks [21, 49], and orthonormalization-free approaches [22, 55]. These aforementioned algorithms are designed for smooth objective functions, and are generally not suitable for problem (2).

There exist some algorithms specifically designed for solving non-smooth optimization problems with orthogonality constraints; most of which adopt certain non-smooth optimization techniques to the Stiefel manifold; for instances, Riemannian subgradient methods [17, 18], proximal point algorithms [6], non-smooth trust-region methods [25], gradient sampling methods [28], and proximal gradient methods [11, 31]. Out of these algorithms, proximal gradient methods are among the most efficient. Specifically, Chen et al. [11] propose a manifold proximal gradient method (ManPG) and its accelerated version ManPG-Ada for non-smooth optimization problems over the Stiefel manifold. Starting from a current iterate, say $Z^{(k)}\in\mathcal{S}_{n,p}$, ManPG and ManPG-Ada generate the next iterate by solving the following subproblem restricted to the tangent space $\mathcal{T}_{Z^{(k)}}\mathcal{S}_{n,p}=\{D\in\mathbb{R}^{n\times p}\mid D^{\top}Z^{(k)}+(Z^{(k)})^{\top}D=0\}$:

which consumes the most computational time in the algorithm. In ManPG [11], the semi-smooth Newton (SSN) method [57] is deployed to solve the above subproblem, and global convergence to stationary points is established with the convergence rate $O(1/\sqrt{k})$. Huang and Wei [31] later extend the framework of ManPG to general Riemannian manifolds beyond the Stiefel manifold. Another class of algorithms first introduces auxiliary variables to split the objective function and the orthogonality constraint and then utilizes alternating minimization techniques to solve the resulting model, including SOC [37], MADMM [36], and PAMAL [12]. It is worth mentioning that SOC and MADMM lack a convergence guarantee. As for PAMAL, although convergence is guaranteed, its numerical performance is heavily dependent on its parameters as discussed in [11].

In principle, all the aforementioned algorithms can be adapted to the distributed setting considered in this paper. We take the algorithm ManPG as an example. Under the distributed setting, each agent computes the local product $A_{i}A_{i}^{\top}Z^{(k)}$ individually, then one round of communications will gather the global sum $\sum_{i=1}^{d}A_{i}A_{i}^{\top}Z^{(k)}$ at a center. The center then solves subproblem (4) and scatters the updated $Z^{(k+1)}$ back to all agents. At each iteration, distributed computation is basically limited to the matrix-multiplication level.

We point out that, without prior data-masking, distributed algorithms at the matrix-multiplication level cannot preserve local-data privacy intrinsically. Specifically, local data $A_{i}A_{i}^{\top}$, privately owned by agent $i$, can be uncovered by anyone who has access to publicly exchanged products $A_{i}A_{i}^{\top}Z^{(k)}$, after collecting enough such products and then solving a system of linear equations for the “unknown” $A_{i}A_{i}^{\top}$. This idea works even for more complicated iterative procedures as long as a mapping from data $A_{i}A_{i}^{\top}$ to a publicly exchanged quantity is deterministic and known.

To illustrate this data leakage situation, we now present a simple experiment on algorithm ManPG-Ada [11]. The test environment well be described in Section 5. As mentioned earlier, at iteration $k$ the publicly shared quantity in ManPG-Ada by agent $i$ is $S_{i}^{(k)}:=A_{i}A_{i}^{\top}Z^{(k)}$, where $Z^{(k)}$ is the $k$-th iterate of ManPG-Ada accessible to all agents. For instance, to recover the unknown local data $A_{1}A_{1}^{\top}$ we construct the following linear system with unknown $Y\in\mathbb{R}^{n\times n}$:

Let $Y^{(k)}$ be the minimum norm least-squares solution to the linear equation (5) at iteration $k$. We perform an experiment to illustrate that $Y^{(k)}$ will quickly converge to $A_{1}A_{1}^{\top}$, when ManPG-Ada is deployed to solve the sparse PCA problem on a randomly generated test matrix with $n=100$ and $m=1280$ (other parameters are $d=10$, $p=10$, and $\mu=0.05$). Figure 1 depicts how the reconstruction error $\|Y^{(k)}-A_{1}A_{1}^{\top}\|_{\mathrm{F}}$ and the stationarity violation¹¹1 Suppose $D^{(k)}$ is the solution to (4). According to Lemma 5.3 in [11], $X^{(k)}$ is a first-order stationary point if $D^{(k)}=0$. Therefore, the stationarity violation is defined as $\|D^{(k)}\|_{\mathrm{F}}$. vary, on a logarithmic scale, as the number of iterations increases. We observe that local data $A_{1}A_{1}^{\top}$ is obtained with high accuracy much faster than solving the sparse PCA problem.

In order to handle distributed datasets, some distributed ADMM-type algorithms have been introduced to PCA and related problems. These methods achieve algorithm-level parallelization, which are generally more secure than those based on matrix-multiplication-level parallelization. For instance, [26] proposes a distributed algorithm for sparse PCA with convergence to stationary points, but only studies the special case of $p=1$. Moreover, under the distributed setting, [51] develop a subspace splitting strategy to tackle the smooth optimization problem over the Grassmann manifold without the $\ell_{1}$-norm regularizer in (2).

Recently, decentralized optimization has attracted attentions partly due to its wide applications in wireless sensor networks. Only local communications are carried out in decentralized optimization, namely, agents exchange information only with their immediate neighbors. This is quite different from the distributed setting considered in the current paper. We refer interested readers to the references [23, 19, 20, 4, 59, 10, 50] for some decentralized PCA algorithms.

Recently, a subspace splitting strategy has been proposed [51] to accelerate convergence of distributed algorithms for optimization problems over Grassmann manifolds where objective functions are orthogonal-transformation invariant²²2A function $f(X)$ is called orthogonal-transformation invariant if $f(XO)=f(X)$ for any $X\in\mathcal{S}_{n,p}$ and $O\in\mathcal{S}_{p,p}$. and smooth. In regularized problems such as sparse PCA, orthogonal-transformation invariance and smoothness no longer hold. In this paper, we present a non-trivial extension of the subspace splitting strategy to more general optimization problems over Stiefel manifolds. In particular, by incorporating the subspace splitting strategy into the framework of ManPG [11], we have developed a distributed algorithm DSSAL1 to solve the $\ell_{1}$-regularized optimization problem (2) for sparse PCA.

The main step of DSSAL1 is to solve the subproblem:

which is identical to subproblem (4) in ManPG except that each local data matrix $A_{i}A_{i}^{\top}$ is replaced, or masked, by a matrix $-Q_{i}^{(k)}$ whose expression will be derived later.

In a sense, the main contributions of this paper can be attributed to this remarkably simple replacement or masking, which brings two great benefits. Firstly, convergence will be significantly accelerated which will be shown empirically in Section 5. Secondly, publicly exchanged products $Q_{i}^{(k)}Z^{(k)}$ are no longer images of some fixed and known mapping of $A_{i}A_{i}^{\top}$, making it practically impossible to uncover $A_{i}A_{i}^{\top}$ from such products via equation-solving.

Several innovations have been made in the development of our algorithms. Firstly, in our algorithm, only the global variable can possibly converge to a solution, while the local variables will never do. The role of the local variables, which are generally dense, is to help identify the right subspace. Secondly, we devise an inexact-solution strategy to effectively handle the difficult subproblem for the global variable where orthogonality and sparsity are pursued simultaneously. Thirdly, we establish the global convergence to stationarity for our algorithm, overcoming a number of technical difficulties associated with the rather complex algorithm construction, as well as the non-convexity and non-smoothness in problem (2).

The Euclidean inner product of two matrices $Y_{1},Y_{2}$ of the same size is defined as $\left\langle Y_{1},Y_{2}\right\rangle=\mathrm{tr}(Y_{1}^{\top}Y_{2})$. The Frobenius norm and spectral norm of a given matrix $C$ are denoted by $\left\|C\right\|_{\mathrm{F}}$ and $\left\|C\right\|_{2}$, respectively. Given a differentiable function $g(X):\mathbb{R}^{n\times p}\to\mathbb{R}$, the gradient of $g$ with respect to $X$ is denoted by $\nabla g(X)$. And the subdifferential of a Lipschitz continuous function $h(X)$ is denoted by $\partial h(X)$. The tangent space to the Stiefel manifold $\mathcal{S}_{n,p}$ at $Z\in\mathcal{S}_{n,p}$, is represented by $\mathcal{T}_{Z}\mathcal{S}_{n,p}=\{Y\in\mathbb{R}^{n\times p}\mid Y^{\top}Z+Z^{\top}Y=0\}$, and the orthogonal projection of $Y$ onto $\mathcal{T}_{Z}\mathcal{S}_{n,p}$ is denoted by $\mathrm{Proj}_{\mathcal{T}_{Z}\mathcal{S}_{n,p}}\left(Y\right)=\left(I_{n}-ZZ^{\top}\right)Y+Z\left(Z^{\top}Y-Y^{\top}Z\right)/2$. For $X\in\mathcal{S}_{n,p}$, we define ${\mathbf{P}}^{\perp}_{X}:=I_{n}-XX^{\top}$ standing for the projection operator onto the null space of $X^{\top}$. Further notation will be introduced as it occurs.

The rest of this paper is organized as follows. In Section 2, we introduce a novel subspace-splitting model for sparse PCA, and investigate the structure of associated Lagrangian multipliers. Then a distributed algorithm is proposed to solve this model based on an ADMM-like framework in Section 3. Moreover, convergence properties of the proposed algorithm are studied in Section 4. Numerical experiments on a variety of test problems are presented in Section 5 to evaluate the performance of the proposed algorithm. We draw final conclusions and discuss some possible future developments in the last section.

In this subsection, we aim to present the stationarity conditions of the sparse PCA problem (2). We first introduce the definition of Clarke subgradient [13] for non-smooth functions.

As discussed in [58, 11], the first-order stationarity condition of (2) can be stated as:

We provide an equivalent description of the above first-order stationarity condition, which will be used in the theoretical analysis.

The proof of Lemma 2.2 is put into Appendix A. Then we can characterize the first-order stationary points of (7) in the following manner.

By associating dual variables $\Gamma_{i}$, $\Lambda_{i}$, and $\Theta$ to the constraints (7b), (7c), and (7d), respectively, we derive an equivalent description of first-order stationarity conditions of (7).

The proof of Proposition 2.4 is relegated to Appendix B. Actually, the equations in (9) can be viewed as the KKT conditions of (7), while $\Gamma_{i}\in\mathbb{R}^{p\times p}$, $\Lambda_{i}\in\mathbb{R}^{n\times n}$, and $\Theta\in\mathbb{R}^{p\times p}$ are the Lagrangian multipliers corresponding to the constraints (7b), (7c), and (7d), respectively.

In [51], a low-rank and closed-form multiplier is devised with respect to the subspace constraint (7c) for the PCA problems, which can be expressed as:

In Appendix (B), we further verify that the above formulation is also valid for the sparse PCA problems at any first-order stationary point $(Z,\{X_{i}\})$. In the next section, we will use (10) to update the multipliers in our framework of ADMM. This strategy simultaneously saves computational costs and storage requirements.

We now describe how to approximate subproblems (12) and (13). By rearrangments, subproblem (12) reduces to:

where $Q^{(k)}=\sum_{i=1}^{d}Q_{i}^{(k)}$ and $Q_{i}^{(k)}$ is an $n\times n$ matrix defined by

We quickly point out here that it is not necessary to construct and store these $Q$-matrices explicitly since we will only use them to multiply $n\times p$ matrices in an iterative scheme to obtain approximate solutions to subproblem (16). More details will follow later.

Apparently, subproblem (16) pursues the orthogonality and sparsity simultaneously, and is not easier to solve than the the original problem (2). However, we will demonstrate later by both theoretical analysis and numerical experiments that inexactly solving (16) by conducting one proximal gradient step is adequate for the global convergence.

Starting from the current iterate $Z^{(k)}$, we first find a decent direction $D^{(k)}$ restricted to the tangent space $\mathcal{T}_{Z^{(k)}}\mathcal{S}_{n,p}$ by solving the following subproblem

where $\eta>0$ is the step size. Since $Z^{(k)}+D^{(k)}$ does not necessarily lie on the Stiefel manifold $\mathcal{S}_{n,p}$, we then perform a projection to bring it back to $\mathcal{S}_{n,p}$, which can be represented as:

Here, the orthogonal projection of a matrix $C\in\mathbb{R}^{n\times p}$ onto $\mathcal{S}_{n,p}$ is denoted by $\mathrm{Proj}_{\mathcal{S}_{n,p}}\left(C\right)=U_{C}V_{C}^{\top}$, where $U_{C}\Sigma_{C}V_{C}^{\top}$ is the economic form of the singular value decomposition of $C$.

Since these $n\times n$ matrices $Q_{i}^{(k)}(i=1,\dotsc,d)$ are distributively maintained in $d$ agents, each agent is not able to independently solve subproblem (18). Fortunately, we only need to calculate

to solve this subproblem. In the distributed setting, the right hand side of (19) can be accomplished by calculating $Q_{i}^{(k)}Z^{(k)}$ in each agent and then invoking the all-reduce type of communication, where each agent just shares one $n\times p$ matrix. If the all-reduce type of communication is realized by the butterfly algorithm [43], the communication overhead per iteration is $O\left(np\log(d)\right)$. Furthermore, each local product $Q_{i}^{(k)}Z^{(k)}$ can be computed from

with a computational cost in the order of $O(np^{2})$. From the above formula, one observes that the $n\times n$ matrices $Q_{i}$ need not be stored explicitly.

Now we consider how to efficiently solve subproblem (18). By associating a multiplier $\Upsilon\in\mathbb{R}^{p\times p}$ to the linear equality constraint, the Lagrangian function of (18) can be written as:

We choose to apply the Uzawa method [5] to solve subproblem (18). At the $j$-th inner iteration, we first minimize the above Lagrangian function with respect to $D$ for a fixed $\Upsilon=\Upsilon(j)$:

where $D(j)$ and $\Upsilon(j)$ denote the $j$-th inner iterate of $D$ and $\Upsilon$, respectively. Here, we use $\mathrm{Prox}_{g}\left(X\right)$ to denote the proximal mapping of a given function $g:\mathbb{R}^{n\times p}\to\mathbb{R}$ at the point $X\in\mathbb{R}^{n\times p}$, which is defined by:

For the $\ell_{1}$-norm regularizer term $r(X)=\mu\left\|X\right\|_{1}$, the proximal mapping in (21) admits a closed-form solution:

where the subscript $[\,\cdot\,]_{ij}$ represents the $(i,j)$-th entry of a matrix. Then the multiplier is updated by a dual ascent step:

where $\tau>0$ is the step size. These two steps are repeated until convergence. The complete framework is summarized in Algorithm 1.

The Uzawa method can be viewed as a special case of the primal-dual hybrid gradient algorithm (PDHG) developed in [27] with the convergence rate $O(1/k)$ in the ergodic sense under mild conditions. Moreover, as an inner solver it can bring a higher overall efficiency than the SSN method used by ManPG [11] (see [56] for a recent study).

In this subsection, we focus on for the $i$-th local variable $X_{i}$, which can be rearranged as the following equivalent problem:

Here, $H_{i}^{(k)}$ is an $n\times n$ real symmetric matrix:

which is only related to the local data $A_{i}$. This is a standard eigenvalue problem where one needs to compute a $p$-dimensional dominant eigenspace of $H_{i}^{(k)}$.

As a subproblem, it is not necessary to solve (23) to high precision. In practice, we just need to find a point $X_{i}^{(k+1)}\in\mathcal{S}_{n,p}$ satisfying the following two conditions, which suffices to be a good inexact solution empirically, and to guarantee the global convergence of the whole algorithm. The first condition demands a sufficient decrease in function value:

where $c_{i}>0$ and $c_{i}^{\prime}>0$ are two constants independent of $\beta_{i}$. The second condition is a sufficient decrease in KKT violation:

where $\delta_{i}\in[0,1)$ is a constant independent of $\beta_{i}$. It turns out that these two rather weak termination conditions for subproblem (23) are sufficient for us to derive global convergence of our ADMM-like algorithm framework (12)-(15).

In practice, we can combine a warm-start strategy with a single iteration of SSI [44] to generate the next iterate $X_{i}^{(k+1)}=\mathrm{Proj}_{\mathcal{S}_{n,p}}\left(H_{i}^{(k)}X_{i}^{(k)}\right)$, i.e.,

which can be computed in the order of $O(np^{2})$ floating-point operations, given that the term $A_{i}A_{i}^{\top}X_{i}^{(k)}$ is inherited from the last iteration as a result of updating $W_{i}^{(k)}$, see (14).

We formally present the detailed algorithmic framework as Algorithm 2 below, named distributed subspace splitting algorithm with $\ell_{1}$ regularization and abbreviated to DSSAL1. In the distributed environment, all agents are initiated from the same point $Z^{(0)}\in\mathcal{S}_{n,p}$. And the initial guess of multipliers are computed by (15). After initialization, all agents first solve the common subproblem for $Z$ collaboratively by certain communication strategy. Then each agent solves its subproblem for $X_{i}$ and updates its multiplier $\Lambda_{i}$. These two steps only involve the local data privately stored at each agent, and hence can be carried out in $d$ agents concurrently. This procedure is repeated until convergence.

We claim that DSSAL1 can naturally protect the intrinsic privacy of local data. To form the global sum in (19), the shared information in DSSAL1 at iteration $k$ from the $i$-th agent is $S_{i}^{(k)}:=Q_{i}^{(k)}Z^{(k)}$ where $Z^{(k)}$ is known to all agents. However, the $n\times p$ system of equations, $Q_{i}^{(k)}Z^{(k)}=S_{i}^{(k)}$, is insufficient for obtaining the $n\times n$ mask matrix $Q_{i}^{(k)}$ which changes from iteration to iteration. Secondly, even if a few mask matrices $Q_{i}^{(k)}$ were unveiled, it would still be impossible to derive the local data matrix $A_{i}A_{i}^{\top}$ from these $Q_{i}^{(k)}$ without knowing corresponding $X_{i}^{(k)}$ (and $\beta_{i}$) which are always kept privately by the $i$-th agent. Finally, consider the ideal “converged” case where $X_{i}X_{i}^{\top}=ZZ^{\top}$ held at iteration $k$, and $\beta_{i}$ were known. In this case, $Q_{i}^{(k)}$ would be a known linear function of $A_{i}A_{i}^{\top}$ parameterized by $Z^{(k)}$. Still, the $n\times p$ system $Q_{i}^{(k)}Z^{(k)}=S_{i}^{(k)}$ would not be sufficient to uncover the $n\times n$ local data matrix $A_{i}A_{i}^{\top}$ (strictly speaking, one only needs to recover $n(n+1)/2$ entries since $A_{i}A_{i}^{\top}$ is symmetric). Based on this discussion, we call DSSAL1 a privacy-aware method.

We conclude this section by discussing the computational cost of our algorithm per iteration. We first compute the matrix multiplication $Q^{(k)}Z^{(k)}$ by (19) and (20), whose computational cost for each agent is $O(np^{2})$ as mentioned earlier. Then, at the center, the Uzawa method is applied to solving subproblem (18) which has a per-iteration complexity $O(np^{2})$. In practice, it usually takes very few iterations to generate $Z^{(k+1)}$. Next, each agent uses a single iteration of SSI to generate $X_{i}^{(k+1)}$ by (27) with the computational cost $O(np^{2})$ as discussed before. Finally, agent $i$ updates $W_{i}^{(k+1)}$ by (14) with the computational cost $4npm_{i}+O(np^{2})$ (which represents the multiplier matrix $\Lambda_{i}^{(k+1)}$ implicitly). Overall, for each agent, the computational cost of our algorithm is $4npm_{i}+O(np^{2})$ per iteration. At the center, the computational cost for approximately solving (18) is, empirically speaking, $O(np^{2})$.