\section{Introduction}
\label{sec:introduction}

When principal components possess a certain sparsity structure, standard principal component analysis (PCA) is not preferred due to poor interpretability and the inconsistency of solutions under high-dimensional settings \cite{paul2007asymptotics, nadler2008finite, johnstone2009consistency}.
To solve these issues, {\it sparse PCA} has been proposed, which enforces sparsity in the PCA solution so that dimension reduction and variable selection can be simultaneously performed. 
Theoretical and algorithmic research on sparse PCA has been actively conducted over the past few years \cite{zou2006sparse, amini2008high, journee2010generalized, berk2019certifiably, 
richtarik2021alternating}.


In this paper, we focus on the case that the data to which sparse PCA is applied, are not completely observed, but partially missing. 
Missing data frequently occurs in a wide range of machine learning problems, and sparse PCA is no exception. 
This has led to several works that offer reliable solutions to sparse PCA on missing data \cite{lounici2013sparse, kundu2015approximating, park2019sparse, lee2022support}.
However, these methods make a restrictive assumption that the entries are observed uniformly at random,
and analyze only the effect of the observation rate on the solvability of the problem.


In real world applications, data missing patterns are usually not at random.
Furthermore, when the study focuses on the case of uniformly random missing, one fails to recognize other important observation-pattern-specific factors (e.g., some topological properties of the realized observation pattern) that may affect the solvability of the problem or the performance of an algorithm.
Accordingly, the analysis under a deterministic missing scheme without any assumptions on its pattern is in demand.
With this motivation, there have been several works in other matrix-related problems which handle a deterministic and general sampling scheme \cite{lee2013matrix, heiman2014deterministic, bhojanapalli2014universal},
while sparse PCA has not received such attention.


Examining a general sampling scheme is particularly important for sparse PCA
because structural properties of the observed entries, not just observation rate,
do affect the solvability of the problem.
Imagine that the leading eigenvector $\pmb{u}$ of a symmetric matrix is sparse and denote its support by $J$.
The goal of our sparse PCA problem is to exactly recover the support $J$ from an incomplete data matrix.
Note that non-zero values in $\pmb{u}$ only affect the entries in the $J\times J$ sub-matrix.
Therefore, to recover the support, we need to observe a sufficient number of entries in the $J\times J$ sub-matrix.
In fact, if we do not observe any entry in a specific row of the $J\times J$ sub-matrix, we can never identify the corresponding index as an element of the support, and thus exact recovery fails.
This implies that we need to observe the entries in the $J\times J$ sub-matrix abundantly and evenly;
in other words, if we think of a graph having the observed entries as its edge set (we call this an observation graph), its $J\times J$ sub-graph needs to be well-connected and have similar node degrees.


With this intuition, 
we aim to explore a sparse PCA algorithm for incomplete data under a general non-random sampling scheme,
which would be successful when observed entries have good structural properties.
The algorithm we consider is a semidefinite relaxation of the $l_1$ regularized PCA (we call this the SDP algorithm).
This algorithm has been analyzed on complete data and has been shown to have theoretically good properties \cite{d2004direct, lei2015sparsistency}.
It has also been shown to work well for incomplete data, when the observation rate is sufficiently large under the uniform random sampling scheme \cite{lee2022support}.


Our main contribution is as follows:
we provide theoretical justification (i.e., \cref{thm:sufficient_conditions}) that we can exactly recover the true support $J$ with high probability by obtaining a unique solution of the SDP algorithm with incomplete and noisy observation, under proper conditions. 
The conditions involve the spectral gap between the largest and second-largest eigenvalues of the true matrix, the magnitude of noise, and especially, structural properties of the observation graph.
We mathematically articulate the structural graph properties by using two interpretable concepts: algebraic connectivity and irregularity
(see \cref{def:algebraic_connectivity} and \cref{def:irregularity}).
Through these graph properties, we demonstrate that the algorithm works well if the $J\times J$ sub-graph is well-connected and has similar node degrees.
We note that the graph properties apply to any type of undirected graph,
that is, our theorem is applicable to any general deterministic sampling scheme.
To the best of our knowledge, this is the first work on sparse PCA with incomplete data under the general non-random sampling assumption.


We empirically validate our theorem with synthetic and real data analysis in \cref{sec:experiments}.
In particular, our simulation results show that the performance of the SDP algorithm is solely determined by the properties we derive in our theorem, 
which is a strong justification of our theory.
In real data analysis, we show that the SDP algorithm outperforms several other sparse PCA approaches, and only the SDP algorithm benefits from the good structure of the observation graph.

As by-products of our analysis, we provide two theorems to handle the non-random sampling scheme: 
the tail bound for random matrices with independent sub-Gaussian values in a subset of entries, and the bound of the difference between complete and incomplete matrices under non-random missingness (see \cref{subsec:by_products}).
These theorems can be applied to other matrix-related problems, when matrix entries are missing with a non-random sampling scheme.




\textbf{Paper Organization.}
In \cref{sec:problem_definition}, we introduce the support recovery problem in sparse PCA with incomplete and noisy observation. 
We discuss several applications of the problem here.
In \cref{sec:methods}, we present the SDP algorithm, and define several graph properties used in the theorem.
We present our main theorem in \cref{sec:main_results}, along with by-products of our analysis and the method to choose a tuning parameter in the SDP algorithm.
In \cref{sec:experiments}, we illustrate results from our empirical study.





\textbf{Notation.} 
Matrices are bold capital (e.g., $\pmb{A}$), vectors are bold lowercase (e.g., $\pmb{a}$), and scalars or entries are not bold.
$A_{i,j}$ and $a_i$ represent the $(i,j)$-th and $i$-th entries of $\pmb{A}$ and $\pmb{a}$, respectively.
For any index sets $I$ and $J$, $\pmb{A}_{I,J}$ and $\pmb{a}_I$ denote
the $|I|\times|J|$-dimensional sub-matrix of $\pmb{A}$ consisting of rows in $I$ and columns in $J$, and
the $|I|$-dimensional sub-vector of $\pmb{a}$ consisting of the entries in $J$,
respectively.
For any positive integer $n$, we denote $[n]:=\{1,\dots,n\}$.

$\|\pmb{a}\|_1$ and $\|\pmb{a}\|_2$ represent the $l_1$ and $l_2$ norms of $\pmb{a}$.
$\|\pmb{A}\|_2$ and $\|\pmb{A}\|_*$ indicate the spectral and nuclear norms of $\pmb{A}$.
We let $\|\pmb{A}\|_{1,1} = \sum_{i}\sum_{j}|A_{i,j}|$ and $\|\pmb{A}\|_{\max} = \max_{i,j}|A_{i,j}|$.
$\lambda_i(\pmb{A})$ represents the $i$-th largest eigenvalue of $\pmb{A}$.
The trace of $\pmb{A}$ is denoted by $tr(\pmb{A})$, and the matrix inner product of $\pmb{A}$ and $\pmb{B}$ is denoted by $\langle \pmb{A}, \pmb{B} \rangle$.
$\pmb{A}\circ\pmb{B}$ represents the Hadamard product of $\pmb{A}$ and $\pmb{B}$.

$f(x) = O(g(x))$ means that there exists a positive constant $C$ such that $f(x)\leq C g(x)$ asymptotically.
$f(x) = \Omega(g(x))$ is equivalent to $g(x) = O(f(x))$.
$f(x) = \tilde{O}(g(x))$ is shorthand for $f(x) = O(g(x)\log^k x)$ for some $k$.




\section{Problem Definition}
\label{sec:problem_definition}

\textbf{Sparse Principal Component.}
Let $\pmb{M}^*\in\mathbb{R}^{d\times d}$ be an unknown symmetric matrix
and 
$\pmb{M}^* = \sum_{k\in [d]}\lambda_k(\pmb{M}^*) \pmb{u}_k \pmb{u}_k^\top$
be the spectral decomposition of $\pmb{M}^*$,
where $\lambda_1(\pmb{M}^*) > \lambda_2(\pmb{M}^*) \geq \cdots \geq \lambda_d(\pmb{M}^*)$ are its eigenvalues and $\pmb{u}_1,\dots,\pmb{u}_d \in \mathbb{R}^d$ are the corresponding eigenvectors.
For identifiability of the leading eigenvector, we consider that $\lambda_1(\pmb{M}^*)$ is strictly greater than $\lambda_2(\pmb{M}^*)$.
We assume that the leading eigenvector $\pmb{u}_1$ of $\pmb{M}^*$ is sparse, i.e., for some index set $J\in [d]$,
$$
\begin{cases}
u_{1,i} \neq 0 & \text{if } i\in J, \\
u_{1,i} = 0 & \text{otherwise.}
\end{cases}
$$	
With a notation $supp(\pmb{a}):=\{i\in[d] : a_{i}\neq 0\}$ for any vector $\pmb{a}\in\mathbb{R}^d$, we can write $J = supp(\pmb{u}_1)$.
Also, we denote the size of $J$ by $s$. That is, $s = |J|$.


\textbf{Incomplete and Noisy Observation.}
Suppose that we have only noisy observations of the entries of $\pmb{M}^*$ over a fixed sampling set $\Omega \subseteq [d]\times[d]$.
Specifically, we observe a symmetric matrix $\pmb{M}\in\mathbb{R}^{d\times d}$ such that 
$$
M_{i,j}
=
\begin{cases}
M^*_{i,j} + N_{i,j}  & \text{if } (i,j)\in \Omega, \\
0 & \text{otherwise}
\end{cases}
$$
for $i,j \in [d]$, where $N_{i,j}$ is the noise at location $(i,j)$.
We assume that $N_{i,j}$'s are symmetric about zero and follow a sub-Gaussian distribution independently, i.e., 
$\mathbb{E}e^{\theta N_{i,j}}\leq e^{\frac{\sigma^2 \theta^2}{2}}$ for any $\theta \geq 0$ and some $\sigma \geq 0$.


\textbf{Goal.}
In this paper, we aim to exactly recover the true support $J$ of the leading eigenvector $\pmb{u}_1$ of $\pmb{M}^*$ only from the incomplete and noisy observation $\pmb{M}$.



\textbf{Applications.}
Here, we point out several motivating applications of the problem.
\begin{itemize}
\item 
\textit{Single-cell RNA sequence (scRNA-seq) data analysis.}
In scRNA-seq data, the cells are divided into several distinct types which can be characterized with only a small number of genes among tens of thousands of genes \cite{park2019sparse}.
Sparse PCA can be effectively utilized here to reduce the dimension (from numerous cells to a few cell types) and to select a small number of genes that affect the reduced data. 
However, scRNA-seq data usually have many missing values due to technical and sampling issues, and a sparse PCA method designed for incomplete data needs to be applied.
\item 
\textit{Covariance analysis in vertical federated setting.}
In vertical federated setting, 
one aims to utilize local data sets with different features about the same set of subjects to train machine learning models \cite{yang2019federated}.
Consider the case that we want to analyze the covariance matrix of the whole data.
Due to communication restriction or privacy protection requirement, only the sample covariance matrix of a few number of features can be provided to the analyst.
In this case, the analyst can only have an incomplete (noise-injected) sample covariance matrix.
Our analysis is applied to the case that sparse PCA is conducted in such setting.
\item 
\textit{Anomaly detection in network.}
Sparse PCA is useful for anomaly detection in network data sets \cite{singh2011anomalous}.
Consider an edge-weighted network data, where the weight for each edge can represent various types of information, e.g., similarity of two nodes, traffic volume between connected locations, or social distance between individuals or groups.
By applying sparse PCA to the weight matrix of the network, we can find an anomalous collection of individuals or nodes.
However, since unconnected nodes cause missingness, analysis for incomplete data must be considered.
\end{itemize}




\section{Methods}
\label{sec:methods}

In this section, we introduce the algorithm used to solve the sparse PCA problem,
and define several graph properties which will be utilized in our main theorem.


\textbf{SDP Algorithm.}
For the support recovery of the leading eigenvector, an intuitive approach is imposing a regularization term on the PCA quadratic loss.
When using the $l_1$ regularizer, the optimization problem can be written as:
\begin{equation*}
\hat{\pmb{x}} = \underset{\|\pmb{x}\|_2 = 1}{\arg\max}~ \pmb{x}^\top \pmb{M}\pmb{x} - \rho \|\pmb{x}\|_1^2.
\end{equation*}
Here, the true support $J$ is estimated with $supp(\hat{\pmb{x}})$.
However, the objective is non-convex and difficult to solve.
Therefore, we suggest the following semidefinite relaxation as an alternative:
\begin{equation}
\label{eq:sdp_problem}
\hat{\pmb{X}} = \underset{\pmb{X}\succeq 0 \text{ and } tr(\pmb{X}) = 1}{\arg\max}~ 
\langle \pmb{M}, \pmb{X} \rangle - \rho \|\pmb{X}\|_{1,1},
\end{equation}
where we estimate $J$ by $\hat{J} = supp(diag(\hat{\pmb{X}}))$.
We call this the SDP algorithm.
Efficient scalable SDP solvers exist \cite{yurtsever2021scalable}, so the SDP algorithm is computationally friendly.
This approach has been shown to have good theoretical properties and work well in practice for both of complete and randomly missing data \cite{d2004direct, lei2015sparsistency, lee2022support},
so we focus on the SDP algorithm in this paper.


\begin{remark}
We note that in the implementation of the SDP algorithm, we use the matrix $\pmb{M}$ where zero is imputed in the missing entries, without applying any matrix completion or imputation methods.
This is because matrix completion can introduce unwanted bias under inappropriate conditions.
It is well-known that most of the matrix completion methods can be successful only under the low-rank assumption.
In this paper, however, we allow the true matrix $\pmb{M}^*$ to be not necessarily low-rank.
In \cref{subsec:real_data_analysis}, we provide experimental evidence  showing that the SDP algorithm with zero-imputed $\pmb{M}$ performs well when the observation has a good structural property,
while the result yielded from matrix completion does not achieve good performance overall.
\end{remark}




\textbf{Graph Properties.}
Before presenting our main theorem in the next section,
we define several graph terminologies and properties which are involved in the theorem.

We first introduce the observation graph $\mathcal{G} = (\mathcal{V},\mathcal{E})$,
which is an undirected graph associated with the fixed sampling set $\Omega$,
that is, $\mathcal{V} = [d]$ and $(i,j)\in \mathcal{E}$ if and only if $(i,j)\in \Omega$.
Note that $\mathcal{G}$ is allowed to contain loops.
We denote the adjacency matrix corresponding to $\mathcal{G}$ by $\pmb{A}_\mathcal{G}$.
With this notation, we can write $\pmb{M} =\pmb{A}_\mathcal{G}\circ(\pmb{M}^* + \pmb{N})$ where $\pmb{N}$ is the noise matrix whose $(i,j)$-th entry is $N_{i,j}$.

Below are several convenient notations about graphs.
\begin{itemize}
\item $\mathcal{G}_{J,J}$, $\mathcal{G}_{J,J^c}$, $\mathcal{G}_{J^c,J^c}$: 
For the observation graph $\mathcal{G}$,
we denote by $\mathcal{G}_{J,J}$, $\mathcal{G}_{J,J^c}$ and $\mathcal{G}_{J^c,J^c}$ the sub-graphs of $\mathcal{G}$ which consist of only the edges inside $J\times J$, $J\times J^c$, and $J^c\times J^c$, respectively. 
$\mathcal{G}_{J,J}$ and $\mathcal{G}_{J^c,J^c}$ are undirected graphs with vertex sets $J$ and $J^c$, respectively.
$\mathcal{G}_{J,J^c}$ is a bipartite graph with independent vertex sets $J$ and $J^c$.
\item $\overline{\mathcal{G}}$: 
For any graph $\mathcal{G}$, $\overline{\mathcal{G}}$ denotes a graph which has the same vertex set as $\mathcal{G}$, but whose edge set is the complement of that of $\mathcal{G}$.
That is, the edge sets of $\mathcal{G}$ and $\overline{\mathcal{G}}$ are disjoint and their union forms a complete graph.
\item $\Delta_{\max}(\mathcal{G})$, $\Delta_{\min}(\mathcal{G})$:
For any graph $\mathcal{G}$, $\Delta_{\max}(\mathcal{G})$ and $\Delta_{\min}(\mathcal{G})$ denote the maximum and minimum node degrees of $\mathcal{G}$, respectively.
\end{itemize}


Now, we define two important structural graph properties, {\it algebraic connectivity} and {\it irregularity}.
Both properties are crucial to explain the effect of the structure of the observation graph on the solvability of our support recovery problem.

Algebraic connectivity, the well-known concept to measure the graph connectivity, is defined as follows:
\begin{definition}[Algebraic Connectivity]
\label{def:algebraic_connectivity}
The algebraic connectivity of a graph $\mathcal{G}$, denoted by $\phi(\mathcal{G})$, is the second-smallest eigenvalue of the Laplacian matrix of $\mathcal{G}$.
The magnitude of $\phi(\mathcal{G})$ reflects how well connected the overall graph is.
\end{definition}

Next, we define the concept of irregularity, which is an uncommon graphical quantity that proves crucial in presenting our results.

\begin{definition}[Irregularity]
\label{def:irregularity}
For any undirected graph $\mathcal{G}$ such that $\Delta_{\max}(\mathcal{G}) \geq \phi(\mathcal{G})$ and $\Delta_{\max}(\overline{\mathcal{G}}) \geq \phi(\overline{\mathcal{G}})$,
the irregularity of $\mathcal{G}$ is defined as
$$
\psi(\mathcal{G}) := \max\big\{\Delta_{\max}(\mathcal{G}) - \phi(\mathcal{G}), \Delta_{\max}(\overline{\mathcal{G}}) - \phi(\overline{\mathcal{G}}) \big\}.$$
The magnitude of $\psi(\mathcal{G})$ reflects how different the node degrees of $\mathcal{G}$ are.
\end{definition}

To better interpret the above concept, we derive some lower and upper bound results below.
Let $\pmb{A}_\mathcal{G}$ be the adjacency matrix of $\mathcal{G}$.
Then we can derive the following inequality\footnotemark:
\begin{align*}
\max_{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}} \pmb{x}^\top \pmb{A}_\mathcal{G} \pmb{x}
\leq 
\Delta_{\max}(\mathcal{G}) - \phi(\mathcal{G})
\leq 
\max_{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}} \pmb{x}^\top \pmb{A}_\mathcal{G} \pmb{x}
+ \Delta_{\max}(\mathcal{G}) - \Delta_{\min}(\mathcal{G})
\end{align*}
where $\pmb{1} = (1, 1, \dots, 1)^\top$.

\begin{itemize}
\item $\max_{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}} \pmb{x}^\top \pmb{A}_\mathcal{G} \pmb{x}$: Among different $\pmb{A}_\mathcal{G}$'s having the same largest and second-largest eigenvalues, we can see that the one corresponding to a regular graph (a graph is regular when each node has the same degree) has the smallest magnitude of $\underset{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}}{\max} \pmb{x}^\top \pmb{A}_\mathcal{G} \pmb{x}$.
This is because a regular graph has a normalized vector of $\pmb{1}$ as its leading eigenvector.
\item $\Delta_{\max}(\mathcal{G}) - \Delta_{\min}(\mathcal{G})$: This quantity decreases as nodes of $\mathcal{G}$ have similar degrees.
\end{itemize}
Since $\overline{\mathcal{G}}$ is regular when $\mathcal{G}$ is regular,
and $\Delta_{\max}(\mathcal{G}) - \Delta_{\min}(\mathcal{G})
= \Delta_{\max}(\overline{\mathcal{G}}) - \Delta_{\min}(\overline{\mathcal{G}})$,
the same argument holds for $\Delta_{\max}(\overline{\mathcal{G}}) - \phi(\overline{\mathcal{G}})$.
Therefore, we can say that as $\mathcal{G}$ is closer to a regular graph, the value of $\psi(\mathcal{G})$ decreases.
This is the reason why we name this concept the ``irregularity".

\footnotetext{
The inequality is derived as follows:
let $\pmb{D}_\mathcal{G}$ be a diagonal matrix whose diagonal entries are the node degrees of $\pmb{A}_\mathcal{G}$. Then,
$
\Delta_{\max}(\mathcal{G}) - \phi(\mathcal{G}) 
= 
\underset{\|\pmb{x}\|_2 = 1}{\max} \pmb{x}^\top \pmb{D}_\mathcal{G} \pmb{x}
- \underset{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}}{\min} \pmb{x}^\top (\pmb{D}_\mathcal{G}-\pmb{A}_\mathcal{G}) \pmb{x}
\geq
\underset{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}}{\max} \pmb{x}^\top \pmb{A}_\mathcal{G} \pmb{x}
$.
Also,
$
\Delta_{\max}(\mathcal{G}) - \phi(\mathcal{G}) 
\leq
\underset{\|\pmb{x}\|_2 = 1}{\max} \pmb{x}^\top \pmb{D}_\mathcal{G} \pmb{x}
+\underset{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}}{\max} \pmb{x}^\top \pmb{A}_\mathcal{G} \pmb{x}
-\underset{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}}{\min} \pmb{x}^\top \pmb{D}_\mathcal{G} \pmb{x}
\leq
\underset{\|\pmb{x}\|_2 = 1, \pmb{x}\perp \pmb{1}}{\max} \pmb{x}^\top \pmb{A}_\mathcal{G} \pmb{x}
+ \Delta_{\max}(\mathcal{G}) - \Delta_{\min}(\mathcal{G}).
$
}



\begin{remark}
Note that $\Delta_{\max}(\mathcal{G}) \geq \phi(\mathcal{G})$ and $\Delta_{\max}(\overline{\mathcal{G}}) \geq \phi(\overline{\mathcal{G}})$, i.e., $\psi(\mathcal{G})\geq 0$ holds in most cases except for some extreme types of graphs, e.g., complete graph without loops.
\end{remark}


\section{Main Results}
\label{sec:main_results}


Now, we introduce our main theorem, which shows the sufficient condition for the SDP algorithm to exactly recover the true support $J$.

\begin{theorem}
\label{thm:sufficient_conditions}
Under the problem definition in \cref{sec:problem_definition},
assume that the following condition holds: with some constant $c>0$,
\begin{align}
&
\|\pmb{M}^*_{J,J}\|_2 \cdot \psi(\mathcal{G}_{J,J})
+ \sigma \sqrt{\Delta_{\max}(\mathcal{G}_{J,J})\log s}
+ s \| \pmb{M}^*_{J^c,J} \|_{2}
+ \frac{1}{\sqrt{s}}\| \pmb{M}^*_{J^c,J^c}\|_2
\nonumber \\&
+ \sigma s\sqrt{ \max\big\{ \Delta_{\max}(\mathcal{G}_{J,J^c}), \Delta_{\max}(\mathcal{G}_{J^c,J^c}) \big\} \log d}
\leq \frac{c\phi(\mathcal{G}_{J,J}) \bar{\lambda}(\pmb{M}^*) \cdot \min_{i \in J}|u_{1,i}|}{s}
\label{eq:suff_cond}
\end{align}
where $\bar{\lambda}(\pmb{M}^*) := \lambda_1(\pmb{M}^*)-\lambda_2(\pmb{M}^*)$ and $\psi(\mathcal{G}_{J,J}) \geq 0$.
If the tuning parameter $\rho$ is properly chosen,
the optimal solution $\hat{\pmb{X}}$ to the optimization problem (\ref{eq:sdp_problem}) is unique
and satisfies $supp(diag(\hat{\pmb{X}})) = J$
with probability at least $1-O(s^{-1})$.
\end{theorem}

In a nutshell, \cref{thm:sufficient_conditions} asserts that the SDP algorithm produces reliable solutions under certain conditions imposed on 
the spectral gap $\bar{\lambda}(\pmb{M}^*)$, the noise intensity parameter $\sigma$, the matrix norms,
and 
the graph properties of the sub-graphs $\mathcal{G}_{J,J}$, $\mathcal{G}_{J,J^c}$ and $\mathcal{G}_{J^c,J^c}$.

To better understand the condition, we consider the setting that $\|\pmb{M}^*_{J,J}\|_2 = O(\bar{\lambda}(\pmb{M}^*))$ and 
$\min_{i \in J}|u_{1,i}| = \Omega(s^{-\frac{1}{2}})$, for instance.
We note that the first inequality holds as long as
$\|\pmb{M}^*_{J,J}\|_2 = \lambda_1(\pmb{M}^*_{J,J})$ and
$\frac{\lambda_2(\pmb{M}^*)}{\lambda_1(\pmb{M}^*)}\leq c$ for some $c\in (0,1)$,
and the second inequality holds when all the non-zero entries of the sparse leading eigenvector are of the same level of magnitude.
In this case, we can rewrite \eqref{eq:suff_cond} as follows:
\begin{align*}
&\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}
= {O}\Big(\frac{1}{s\sqrt{s}} \Big),
\\&
\sigma 
=  
\tilde{O} \Big( \frac{\phi(\mathcal{G}_{J,J})}{\sqrt{\max\{ \Delta_{\max}(\mathcal{G}_{J,J^c}), \Delta_{\max}(\mathcal{G}_{J^c,J^c}) \}}}\cdot \frac{\bar{\lambda}(\pmb{M}^*)}{s^2 \sqrt{ s}} \Big),
\\&
\| \pmb{M}^*_{J^c,J} \|_{2} 
=  
{O} \Big( \frac{\phi(\mathcal{G}_{J,J}) \bar{\lambda}(\pmb{M}^*)}{s^2\sqrt{s}}  \Big),
\\&
\| \pmb{M}^*_{J^c,J^c}\|_2 
=  
{O} \Big( \frac{\phi(\mathcal{G}_{J,J}) \bar{\lambda}(\pmb{M}^*)}{s} \Big).
\end{align*}

The first condition about the structural graph properties of $\mathcal{G}_{J,J}$ states that
for the algorithm to be successful, 
the sub-graph $\mathcal{G}_{J,J}$ is desired to have sufficiently large connectivity $\phi(\mathcal{G}_{J,J})$ and small irregularity $\psi(\mathcal{G}_{J,J})$.
This implies that we need to observe the entries in the $J\times J$ sub-matrix of the true matrix abundantly and evenly.
This result fits well with our first intuition discussed in the introduction.

The second through fourth conditions mean that the noise and the norms of $\pmb{M}^*_{J^c,J}$ and $\pmb{M}^*_{J^c,J^c}$ need to be well-controlled for the success of the algorithm.
This is in accordance with our common sense, since large $\sigma$ or $\pmb{M}^*$ values outside $J\times J$ matrix will mask the true information.
These conditions are alleviated when the connectivity of $\mathcal{G}_{J,J}$ is large (especially in the second condition, larger than the number of observed entries outside $J\times J$) and when the spectral gap $\bar{\lambda}(\pmb{M}^*)$ is large.

It is worth mentioning that a sufficiently large spectral gap requirement is to ensure the uniqueness and identifiability of the projection matrix with respect to the principal subspace, which has been also discussed in \citet{lei2015sparsistency} and \citet{lee2022support}.




\textbf{Proof Technique.}
Detailed proof of \cref{thm:sufficient_conditions} is given in \cref{sec:proof_of_thm1}.
At a high level, we use the KKT conditions under the primal-dual witness framework,
and find the sufficient conditions which guarantee that the solution of \eqref{eq:sdp_problem} is unique and satisfies $supp(diag(\hat{\pmb{X}})) = J$ with high probability.
We apply several techniques in the derivation, including Davis Kahan sin$\Theta$ theorem and Weyl's inequality for principal subspace estimation.

One challenge in the proof is to handle a deterministic sampling scheme.
We note that traditional concentration inequalities, such as matrix Bernstein inequality, are not useful to derive tail bounds under non-random missingness.
To overcome this, we obtain and utilize two theorems: one is for the tail bound of the matrix whose subset of entries are random, and the other is for bounding the difference between complete and incomplete matrices under a deterministic sampling scheme.
These two results can be widely used in other matrix-related problems for non-random missing data as well.
We introduce these two theorems as by-products of our analysis in the next section.




\subsection{By-products}
\label{subsec:by_products}

Below is a tail bound for the matrix containing independent sub-Gaussian random values only in a fixed subset of entries. 

\begin{theorem}[Tail Bound for Random Matrix with Independent Sub-Gaussian Values in a Subset of Entries]
\label{thm:by_product_tail_bound}
Consider a random matrix $\pmb{Z}\in \mathbb{R}^{m\times n}$ whose subset of entries independently follow sub-Gaussian distributions which are symmetric about zero and have parameter $\sigma > 0$, while the other entries are fixed as zero.
That is, there exists an index set $S \subseteq \{(i,j)~|~ i \in [m], j \in [n]\}$ such that for $i \in [m]$ and $j \in [n]$,
$$
Z_{i,j} =  
\begin{cases}
N_{i,j} & \text{if } (i,j)\in S, \\
0 & \text{otherwise},
\end{cases}
$$
where each $N_{i,j}$ is symmetric about zero and satisfies $\mathbb{E}e^{\theta N_{i,j}}\leq e^{\frac{\sigma^2 \theta^2}{2}}$ for any $\theta >0$.
Then for any $t \geq 0$,
$$
\mathbb{P} [\| \pmb{Z} \|_2 \geq t ]
\leq 
2(m+n)\cdot \exp\Big(-\frac{t^2}{2\sigma^2 \Delta_{\max}(\mathcal{G}_S)}\Big),
$$
where $\mathcal{G}_S$ is a bipartite graph whose vertex and edge sets are $[m]\times[n]$ and $S$, respectively.
\end{theorem}

We defer the proof to \cref{sec:proof_of_tail_bound}.
In the derivation of \cref{thm:sufficient_conditions}, we use the above theorem to obtain the tail bounds of $\|(\pmb{A}_\mathcal{G})_{J,J}\circ \pmb{N}_{J,J}\|_2$, $\|(\pmb{A}_\mathcal{G})_{J,J^c}\circ \pmb{N}_{J,J^c}\|_2$ and $\|(\pmb{A}_\mathcal{G})_{J^c,J^c}\circ \pmb{N}_{J^c,J^c}\|_2$.





Next is the bound of the difference between complete and incomplete matrices under non-random missingness.
We note that this is an extended result of Theorem 4.1 in \citet{bhojanapalli2014universal}.
While the theorem of \citet{bhojanapalli2014universal} is limited to the case that the observation graph is regular,
our theorem applies to any general undirected observation graph.
For regular graphs, our bound coincides with that of \citet{bhojanapalli2014universal}, i.e., our result generalizes the result of \citet{bhojanapalli2014universal}.



\begin{theorem}
\label{thm:by_product_diff_bound}
Consider a symmetric matrix $\pmb{Y}$ with dimension $n$.
Let $\pmb{Y} = \sum_{k\in [r]} \lambda_k(\pmb{Y}) \pmb{v}_k \pmb{v}_k^\top$ be the spectral decomposition of $\pmb{Y}$, where $r$ is rank of $\pmb{Y}$.
Define $\tau := \max_{i\in[n]} \sum_{k\in [r]} {v}^2_{k,i}$.

Also, consider an undirected graph $\mathcal{G}$ with $n$ nodes and denote its adjacency matrix by $\pmb{A}_\mathcal{G}$.
Then,
$$
\|\pmb{Y} - \frac{n}{\phi(\mathcal{G})}\pmb{A}_\mathcal{G}\circ \pmb{Y} \|_2 \leq \frac{n\tau\psi(\mathcal{G})}{\phi(\mathcal{G})}\cdot  \|\pmb{Y}\|_2.
$$
\end{theorem}

The proof is given in \cref{sec:proof_of_diff_bound}.
We use the above theorem to bound $\| \pmb{M}^*_{J,J} - \frac{s}{\phi(\mathcal{G}_{J,J})} (\pmb{A}_\mathcal{G})_{J,J}\circ \pmb{M}^*_{J,J} \|_2$ in the proof of \cref{thm:sufficient_conditions}.



\subsection{Choice of the Tuning Parameter}
\label{subsec:choice_of_tuning_parameter}

The theoretical choice of $\rho$ \eqref{eq:theoretical_choice_of_rho} is useless in practice since it relies on unknown quantities. 
Therefore, certain tuning procedure over $\rho$ is necessary for the implementation of the SDP algorithm \eqref{eq:sdp_problem}.
Our suggestion is to find $\rho$ to maximize the following AIC type criterion (see also \citet{qi2013sparse}):
$$
C_\rho = (1-a)\frac{\langle \pmb{M}, \hat{\pmb{X}}_\rho \rangle}{\langle \pmb{M}, \hat{\pmb{X}}_0 \rangle}
+ a\Big(1-\frac{|supp(diag(\hat{\pmb{X}}_\rho))|}{d}\Big).
$$
Here, $\hat{\pmb{X}}_\rho$ and $\hat{\pmb{X}}_0$ refer to the solutions of the SDP algorithm where the tuning parameters are set to be $\rho$ and $0$, respectively.
$\langle \pmb{M}, \hat{\pmb{X}}_\rho \rangle$ and $\langle \pmb{M}, \hat{\pmb{X}}_0 \rangle$ represent the explained variances of $\pmb{M}$ by the solutions.

The first term of the criterion is a measure for the quality of the estimate, and the second term penalizes for the complexity of the solution.
$a\in (0,1)$ is the weight to be chosen by practitioners. 
As one needs a sparse principal component, a relatively large value of $a$ is suggested.
In the experiments, we find that $0.4\leq a \leq 0.6$ generally work well.



\section{Experiments}
\label{sec:experiments}

In this section, we present several empirical results on synthetic and real data sets which support our theoretical results.

\subsection{Simulations}
\label{subsec:simulations}

The goal of this simulation study is to demonstrate the effects of the structural properties of the observation graph, the spectral gap between the largest and second-largest eigenvalues of the true matrix, and the magnitude of the noise, on the success of the support recovery by the SDP algorithm.

In particular, we want to check whether the performance of the SDP algorithm is solely determined by the properties we derive in \cref{thm:sufficient_conditions}.
We utilize the following rescaled parameter for this:
\begin{equation}
\label{eq:rescaled}
Rescaled = \frac{\text{LHS of \eqref{eq:suff_cond}}}{\text{RHS of \eqref{eq:suff_cond} without constant } c}.
\end{equation}
If the performance of the algorithm versus this rescaled parameter is the same across different settings, then we can empirically justify that 
the performance is solely determined by the factors in the rescaled parameter.
This kind of approach has been used in \citet{wainwright2009sharp} for sparse linear regression.



\textbf{Setting.}
We use synthetic data sets generated in the following manner.
The orthonormal eigenvectors of $\pmb{M}^*$ are randomly selected, 
while the leading eigenvector $\pmb{u}_1$ is made to be sparse and have $s$ non-zero entries with the values of $\frac{1}{\sqrt{s}}$.
$\lambda_2(\pmb{M}^*), \dots, \lambda_d(\pmb{M}^*)$ are randomly selected from a normal distribution with mean $0$ and standard deviation $1$,
and $\lambda_1(\pmb{M}^*)$ is set to $\lambda_2(\pmb{M}^*)$ plus the spectral gap.
We set the matrix dimension $d$ to be $50$ and the support size $s$ to be $10$.

We generate the observation graph $\mathcal{G}$ to have $1250$ edges out of $2500$.
The value of $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ is set to be included in one of the ranges $0$ to $2$, $2$ to $4$, ..., or $16$ to $18$.
The entry-wise noise is randomly selected from a normal distribution with mean $0$ and standard deviation $\sigma$.

In each setting, we run the algorithm \eqref{eq:sdp_problem} and examine if the solution exactly recovers the true support $J$.
The tuning parameter $\rho \in \{0.025, 0.5, \dots, 1\}$ is selected by the method in \cref{subsec:choice_of_tuning_parameter} with $a=0.5$.
We repeat each experiment $100$ times with different random seeds, and calculate the rate of exact recovery in each setting.


\begin{figure}[t]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=0.7\columnwidth]{figures/simulation1.png}}
\caption{(a) Rate of exact recovery of support $J$ versus the graph parameter $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ for four different spectral gaps $\bar{\lambda}(\pmb{M}^*)$,
(b) Same simulation results with exact recovery rate plotted versus the rescaled parameter in \eqref{eq:rescaled}.}
\label{fig:simulation1}
\end{center}
\vskip -0.2in
\end{figure}


\textbf{Results.}
\cref{fig:simulation1} shows the experimental results where we fix the noise parameter $\sigma$ as $0$ (noiseless) and try different spectral gaps $\bar{\lambda}(\pmb{M}^*) \in \{1, 2, 5, 10\}$ to check the effect of the spectral gap.
From (a) in \cref{fig:simulation1}, we can observe that the exact recovery rate increases as the spectral gap increases and the value of $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ decreases, which is consistent with our theoretical findings.


\begin{figure}[t]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=0.7\columnwidth]{figures/simulation2.png}}
\caption{(a) Rate of exact recovery of support $J$ versus the graph parameter $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ for four different noise parameters $\sigma$,
(b) Same simulation results with exact recovery rate plotted versus the rescaled parameter in \eqref{eq:rescaled}.}
\label{fig:simulation2}
\end{center}
\vskip -0.2in
\end{figure}


\cref{fig:simulation2} shows the experimental results where we fix the spectral gap $\bar{\lambda}(\pmb{M}^*)$ as $20$ and try different noise parameters $\sigma \in \{0.1, 0.3, 0.5, 0.7\}$ to check the effect of the magnitude of noise.
From (a) in \cref{fig:simulation2}, we can observe that the exact recovery rate increases as the standard deviation of the noise decreases and the value of $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ decreases, which also supports our theorem.

In \cref{fig:simulation1} (b) and \cref{fig:simulation2} (b),
we can see that the curves of the exact recovery rate versus the rescaled parameter share almost the same pattern under different settings of $\bar{\lambda}(\pmb{M}^*)$ and $\sigma$.
This provides empirical justification of our theorem in the sense that the performance of the SDP algorithm is solely determined by the properties we derive in \cref{thm:sufficient_conditions}.










\subsection{Real Data Analysis}
\label{subsec:real_data_analysis}

The primary goal of this experimental study is to check if our SDP algorithm performs well compared to other sparse PCA algorithms on incomplete data.
We will also conduct a study to validate the selection criterion of the tuning parameter $\rho$ in \cref{subsec:choice_of_tuning_parameter}.

\textbf{Pitprops Data.}
The pitprops data, which stores 180 observations of 13 variables, has been a standard benchmark to evaluate algorithms for sparse PCA (see, e.g., 
\citet{zou2006sparse, shen2008sparse, journee2010generalized, qi2013sparse}).
We aim to compute a leading principal component of this data.
It has been revealed that on the complete data, a sparse solution with $6$ nonzero entries (with respect to the variables `topdiam', `length', `ringbut', `bowmax', `bowdist', `whorls') has a comparable explained variance with that of the dense solution from original PCA.
We focus on recovering these $6$ nonzero entries with incomplete data, by generating observation graphs synthetically.


\textbf{Setting.}
We impose missingness and noise on the complete $13\times 13$ covariance matrix in the following manner.
We generate the observation graph $\mathcal{G}$ to have $100$ edges out of $169$.
The value of $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ is set to be included in one of the ranges $0$ to $0.2$, $0.2$ to $0.4$, ..., or $2$ to $2.2$.
The entry-wise noise is randomly selected from a normal distribution with mean $0$ and standard deviation $\sigma=0.1$.


We compare our SDP algorithm with three different methods.
Firstly, we consider two popular sparse PCA algorithms: the diagonal thresholding sparse PCA (DTSPCA) by \citet{johnstone2009consistency} and
the iterative thresholding sparse PCA (ITSPCA) by \citet{ma2013sparse},
which are efficient in computation and have theoretically good properties on complete data.
When implementing these methods on incomplete data, we treat missing cells as zero.
Secondly, we consider the combination of matrix completion and the SDP algorithm.
We first estimate the missing entries of the incomplete matrix $\pmb{M}$ by using the following matrix completion method based on the nuclear norm minimization:
$$
\tilde{\pmb{M}} = \underset{\pmb{Y}= \pmb{Y}^\top,~ \pmb{A}_\mathcal{G}\circ \pmb{Y} = \pmb{A}_\mathcal{G}\circ \pmb{M}}{\arg\min} \| \pmb{Y}\|_*.
$$
Then we implemented the SDP algorithm with the completed matrix $\tilde{\pmb{M}}$.


We run each algorithm $100$ times with different random seeds in each setting, and calculate the rate of exact recovery of $J$.
In the SDP algorithm, the tuning parameter $\rho \in \{0.025, 0.5, \dots, 1\}$ is selected by the method in \cref{subsec:choice_of_tuning_parameter} with $a=0.4$.
For the DTSPCA and ITSPCA algorithms, there is no well-known method to choose tuning parameters, so we try multiple values of tuning parameters and choose one with the largest exact recovery rate.


\begin{figure}[t]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=0.7\columnwidth]{figures/tuning.png}}
\caption{Rate of exact recovery of support $J$ versus the graph parameter $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ for SDP algorithm with different values of the tuning parameter $\rho$.
The thick blue line indicates the result of the case that $\rho$ is selected by the criterion $C_\rho$ in \cref{subsec:choice_of_tuning_parameter}.}
\label{fig:tuning}
\end{center}
\vskip -0.2in
\end{figure}

\textbf{Results.}
In \cref{fig:tuning}, we can see that our method of selecting $\rho$ works well.
Here, we compare the result from our tuning method with those in the settings where the tuning parameter $\rho$ is fixed as a value among $\{0.05, 0.1, \dots, 1\}$ through all the repetitions.
We can see that the exact recovery rate from our tuning method is larger than most of the results where $\rho$ is fixed as one value.


\begin{figure}[t]
\vskip 0.2in
\begin{center}
\centerline{\includegraphics[width=0.5\columnwidth]{figures/comparison.png}}
\caption{Rate of exact recovery of support $J$ versus the graph parameter $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ for four different sparse PCA algorithms.}
\label{fig:comparison}
\end{center}
\vskip -0.2in
\end{figure}


\cref{fig:comparison} shows that our SDP method outperforms the other sparse PCA methods when the observation graph has a good structural property.
First, we observe that DTSPCA and ITSPCA algorithms do not work well at all (the exact recovery rates are always less than 0.2), even when the graph parameter $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ is small enough.
Unlike our SDP method, those two algorithms can not benefit from good structure of the observation graph.

In addition, the SDP method with matrix completion has the exact recovery rate of around $0.4$ overall,
while our algorithm with zero-imputation produces the exact recovery rates greater than $0.6$ when the value of $\frac{\psi(\mathcal{G}_{J,J})}{\phi(\mathcal{G}_{J,J})}$ is small enough.
About the failure of the matrix completion approach, we conjecture the following rationale:
while the matrix completion algorithm can be successful under the low-rank assumption, the pitprops data is not low-rank,
which yields an unsuccessful matrix completion result.
Accordingly, the imputed cells introduce more noise into the inference and the result of sparse PCA becomes even worse than that of simply using zero for the missing entries.
We note that unlike matrix completion, our algorithm does not require the low-rank assumption to be successful.
Therefore, our algorithm has good performance under a general condition on the true matrix.





\section{Concluding Remarks}
\label{sec:concluding_remarks}


This paper examines the support recovery problem in sparse PCA with incomplete and noisy data, under a general and deterministic sampling scheme.
The problem formulation spans various applications including gene expression data analysis, vertical federated learning, and anomaly detection in networks.
We suggest a practical algorithm based on a semidefinite relaxation of the $\ell_1$-regularized PCA problem,
and provide sufficient conditions where the algorithm can exactly recover the true support with high probability.
The conditions involve the spectral gap, the noise parameter, the matrix norms, and the structural properties of the observation graph.
We show that the algorithm works well if we observe the entries in the $J\times J$ sub-matrix abundantly and evenly.
We empirically justify our theorem with synthetic and real data analysis,
and show that our algorithm outperforms several other sparse PCA approaches especially when the observation graph has a good structural property.

Our work primarily focuses on the theoretical understanding of the problem and the algorithm,
but we briefly discuss a practical use of our results here.
Our results present that as we observe the entries well in the $J\times J$ sub-matrix, we can recover the support $J$.
Unfortunately, we do not know the true support $J$ in practice, so it is impossible to check if the observation structure satisfies the sufficient condition.
One possible alternative is to check if any sub-graphs of the observation graph have good algebraic connectivity and irregularity.
If there is a sub-graph with connectivity or irregularity that is extremely low or high, then we could conservatively suspect that the result from the algorithm can not be fully trusted.
Here, we will need prior knowledge about the size of the support.





