The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
79,292 characters
Synthetic Nearest Neighbors: Extending Synthetic Controls for Matrix Completion with Missing Not at Random Data
\maketitle
\vspace{-15pt}
\begin{abstract}
We develop a causal framework for matrix completion under missing not at random (MNAR) data. Drawing on synthetic controls from the econometric panel data literature, our approach relaxes two assumptions common in MNAR matrix completion: positivity and independence of observation indicators.
Unlike traditional panel data models, which often require prescribed block-sparse geometries, our framework accommodates flexible, heterogeneous observation patterns through target-specific local information structures.
We propose synthetic nearest neighbors (\texttt{SNN}), a local synthetic-controls-inspired estimator, and establish finite-sample entrywise error bounds and consistency for mean recovery under suitable conditions.
We further derive asymptotic normality under heteroskedastic noise and develop feasible entrywise inference.
To estimate entry-specific noise variances, we apply the same local principle to squared outcomes, obtaining consistency under bounded noise and asymptotic unbiasedness under general subgaussian noise.
Simulation studies corroborate the theoretical findings across a range of missingness designs and observation patterns.
\end{abstract}
\medskip
\noindent
{\smaller
\textbf{Keywords:} heteroskedastic variance estimation; panel data; entrywise inference
}
\section{Introduction} \label{sec:intro}
Matrix completion seeks to recover an underlying matrix from partial, often noisy, observations.
Owing to its broad applicability, this field has witnessed rapid growth.
Traditionally, theoretical guarantees in matrix completion are established with respect to global error metrics, such as the Frobenius norm, which quantify the average performance across all entries.
While these metrics are important, they can obscure localized errors and fail to provide actionable insights in settings where individual predictions are consequential.
For instance, in personalized recommendations, accuracy at the level of specific entries---not just on average---is often of chief concern.
Moreover, these recovery guarantees typically rest on two core assumptions:
(i) the underlying matrix has a low-complexity structure (e.g., low-rank), and
(ii) the entries are {\em missing completely at random} (MCAR), i.e., each entry is observed independently with uniform probability $p > 0$, regardless of its value.
In many practical applications, however, the MCAR assumption is frequently violated.
One motivating example arises in recommender systems, a canonical application of matrix completion.
Here, user-item interactions are represented as a matrix, where rows correspond to users, columns to items,
and the $(i,j)$th entry records user $i$'s rating of item $j$.
Observations in such systems are often subject to {\em selection bias}.
Sports enthusiasts are more likely than non-fans to watch {\em The Last Dance}; meat-lovers are more likely than vegetarians to review steakhouses. Platform-side interventions can amplify these patterns: users searching for Grand Canyon trails are more likely to see ads for hiking boots than for wedding heels.
Across these cases, the user's preferences and the system's inferred beliefs about those preferences shape the sparsity pattern in such a way that the observed entries are influenced by the underlying values of the matrix itself---a phenomenon known as {\em missing not at random} (MNAR).
Recent works have made significant strides in addressing MNAR matrix completion, demonstrating that algorithms explicitly accounting for MNAR mechanisms often outperform standard methods designed under the MCAR assumption.
Yet, on the theoretical front, key questions remain open.
As highlighted by \cite{ma2019missing}, two major limitations persist in the existing literature:
(i) {\em positivity}---every entry must have a strictly positive probability of being observed;
(ii) {\em independence}---entries are assumed to be observed independently of one another.
This work aims to challenge both assumptions, which can prohibit applicability in many real-world settings.
To move beyond these constraints, we build on a growing line of research that connects matrix completion with the panel data framework in econometrics.
In panel data, rows index units (e.g., states) and columns index time-treatment pairs (e.g., economic policies over time) such that the $(i, (d,t))$th entry records the outcome of unit $i$ under treatment $d$ at time $t$.
Like recommender systems, panel data often exhibits MNAR patterns, as policymakers may selectively implement programs based on expected outcomes.
However, the structure of policy implementation introduces additional sparsity constraints: only one policy can be applied at a time, and policies, once introduced, generally cannot be removed.
These constraints clearly violate positivity and independence, but they also induce {\em block sparsity patterns} that are more structured and more restrictive than the irregular sparsity patterns that naturally arise in matrix completion problems.
Thus, the panel data literature offers frameworks that relax positivity and independence, but often under restrictive observation geometries. Standard matrix completion allows richer and more irregular observation patterns, but usually relies on stronger probabilistic assumptions.
This intersection motivates the central question of this paper:
{\em can we obtain entrywise guarantees for matrix completion under MNAR data, without relying on the limiting assumptions of positivity and independence, and for a broad class of observation patterns?}
\subsection{How the Missingness Mechanism Affects Matrix Recovery: a Teaser} \label{sec:teaser}
We begin with a simple illustration designed to isolate the role of the missingness mechanism. To remove the confounding effect of noise, we consider a noiseless setting and compare recovery under three observation designs.
\begin{enumerate} [label=(\alph*)]
\item {\bf MCAR}: Each entry is observed independently with probability $p=0.35$.
\item {\bf Conventional MNAR}: Each entry is observed with its own probability that is correlated with the underlying value. The standard assumptions of positivity and independence are maintained.
\item {\bf Panel-based MNAR}: Observation probabilities again depend on the underlying values, but now the two classical assumptions fail: some entries may have zero probability of being observed, and observations may be dependent across entries.
\end{enumerate}
Figure~\ref{fig:sparsity_patterns} reveals the observation patterns induced by these three mechanisms.
A recommender system interpretation provides useful intuition. Under MCAR, users rate movies entirely at random. Under conventional MNAR, users are more likely to rate movies they strongly liked or disliked and less likely to rate movies toward which they were indifferent.
Under panel-based MNAR, users rate a common set of core movies and, among the remaining titles, only those in their preferred genres. This represents a more extreme form of selection. Detailed constructions of designs (b) and (c) are provided in Appendix~\ref{sec:teaser.details} of the Supplementary Material.
As benchmarks, we consider two widely used matrix completion algorithms: universal singular value thresholding (\texttt{USVT}) of \cite{Chatterjee15}, a spectral method, and \texttt{softImpute}~of \cite{softimpute14}, an optimization-based method. Both were originally designed for MCAR data. To adapt them to MNAR settings, we apply standard debiasing techniques following \cite{bhattacharya2021matrix} and \cite{ma2019missing}. We also include a preview of our proposed method, {\em synthetic nearest neighbors} (\texttt{SNN}), which is introduced formally in Section~\ref{sec:estimator}.
The regularization hyperparameters for de-biased \texttt{softImpute}~and \texttt{SNN}~are selected using five-fold cross-validation.
Figure~\ref{fig:teaser.err} reports the entrywise root mean squared error (RMSE) for each method. Under MCAR, de-biased \texttt{softImpute}~and \texttt{SNN}~achieve substantially lower error than modified \texttt{USVT}. Under the conventional MNAR mechanism, \texttt{SNN}~remains accurate, whereas both benchmark methods incur considerably larger errors. Under panel-based MNAR, \texttt{SNN}~continues to have the lowest average error, although its performance is more dispersed across trials.
These results illustrate that the observation mechanism can materially affect matrix recovery. Even methods equipped with standard debiasing corrections may deteriorate when positivity or independence fails. This motivates the development of procedures that remain reliable across a broad range of MNAR mechanisms and observation patterns.
\begin{figure*}[t!]
\centering
\begin{subfigure}{0.22\linewidth}
\centering
\includegraphics[width=\textwidth]
{images/teaser/mcar.png}
\caption{MCAR.}
\label{fig:sparsity_MCAR}
\end{subfigure}
\quad
\begin{subfigure}{0.22\linewidth}
\centering
\includegraphics[width=\textwidth]
{images/teaser/c_mnar.png}
\caption{Conventional MNAR.}
\label{fig:sparsity_limited_MNAR}
\end{subfigure}
\quad
\begin{subfigure}{0.22\linewidth}
\centering
\includegraphics[width=\textwidth]
{images/teaser/p_mnar.png}
\caption{Panel-based MNAR.}
\label{fig:sparsity_general_MNAR}
\end{subfigure}
\caption{Sparsity patterns induced by different missingness mechanisms.}
\label{fig:sparsity_patterns}
\end{figure*}
\begin{figure}[t!]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]
{images/teaser/err_mcar.png}
\caption{MCAR.}
\label{fig:teaser_MCAR_USVT}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]
{images/teaser/err_c_mnar.png}
\caption{Conventional MNAR.}
\label{fig:teaser_MCAR_softimpute}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]
{images/teaser/err_p_mnar.png}
\caption{Panel-based MNAR.}
\label{fig:teaser_MCAR_SNN}
\end{subfigure}
\caption{Entrywise RMSE for modified \texttt{USVT}, de-biased \texttt{softImpute}, and \texttt{SNN}~under (a) MCAR, (b) conventional MNAR, and (c) panel-based MNAR.}
\label{fig:teaser.err}
\end{figure}
\subsection{Contributions}
This paper develops a causal framework for matrix completion under a broad class of MNAR mechanisms. The framework is inspired by the panel data literature, but is designed for the more irregular sparsity patterns that arise in matrix completion.
Within a low-rank factor model, we provide an identification principle for individual entries. We show that an entry can be identified from observed data even when its probability of observation may be zero, when observation indicators may be dependent across entries, and when the observation mechanism may be correlated with the underlying signal.
This directly relaxes the positivity and independence assumptions that underlie much of the existing MNAR matrix completion literature. Relative to standard panel data models, our approach also relaxes the requirement of a prescribed block-sparse geometry.
Instead, identification is target-specific: a given entry is recoverable whenever the observed mask contains a {\em sufficiently informative local pattern} around that entry (see the conditions of Theorem~\ref{thm:consistency} and the subsequent discussion).
To operationalize this framework, we introduce \texttt{SNN}, an estimator that blends the local nature of nearest neighbor methods with the reweighting logic of synthetic controls.
Under suitable assumptions, we establish finite-sample entrywise error bounds and consistency for mean
recovery (Theorem \ref{thm:consistency}).
We also show that the theory specializes to the MCAR setting as a benchmark, and we exhibit an observation design under which a shared local information pattern yields uniform entrywise recovery with transparent sample-complexity trade-offs.
Beyond point estimation, we prove asymptotic normality under heteroskedastic noise (Theorem \ref{thm:normality})
and develop feasible studentized inference (Theorem \ref{thm:var.est.asymp}).
The accompanying variance estimator applies the \texttt{SNN}~principle to the squared observed matrix.
We prove that this estimator is consistent under bounded noise and asymptotically unbiased under general subgaussian noise (Theorem \ref{thm:var.est}).
Simulation studies support the theoretical results across a range of missingness mechanisms and observation patterns.
\subsection{Related Work}
\noindent \textbf{Matrix completion.}
Low-rank matrix completion has been extensively studied under the MCAR assumption.
Within this regime, a wide array of methods have been proposed, including optimization-based approaches \citep{CandesTao10, Recht11, softimpute}, spectral techniques \citep{KeshavanMontanariOh10a, donoho14, Chatterjee15}, and collaborative filtering algorithms \citep{goldberg1992using, koren2015advances}.
For a comprehensive overview of standard assumptions and theoretical guarantees, see \cite{ieee_matrix_completion_overview}.
A more nuanced paradigm is {\em missing at random} (MAR), where the missingness is allowed to depend on observed covariates.
Conditional on the covariates, the outcomes and missingness mechanism are assumed to be independent, an assumption also known as unconfoundedness or selection on observables in the causal inference literature \citep{imbens_rubin_2015}.
This naturally leads to methods that estimate propensity scores based on covariates.
Key contributions in this regime include \cite{Liang2016, schnabelfwang16, WangNeurips2018} and \cite{wang2019}.
The most general and challenging regime is MNAR, which encompasses all mechanisms beyond MCAR and MAR.
As discussed in \cite{ma2019missing}, existing MNAR approaches hinge on positivity and independence.
A common strategy under these assumptions is to model the propensity score matrix as low-rank and estimate it using matrix completion techniques on the observation mask \citep{ma2019missing, bhattacharya2021matrix}.
This generalizes MAR approaches by removing the dependence on covariate information.
Additional works addressing MNAR mechanisms include \cite{sportisse2020estimation_PCA, causal_recommender_systems, foucart} and \cite{zhu2019high}.
Many prior works in the matrix completion literature do not explicitly distinguish MAR from MNAR; for a detailed discussion of this distinction, see \cite{little2019statistical}.
Most existing analyses focus on global estimation error metrics. A few works provide entrywise error bounds \citep{LeeLiShahSong16, chen_2020, chen_2021}, but their analyses rely heavily on independent sampling and therefore do not directly cover settings with structured missingness such as ours. Their bounds also scale inversely with a polynomial of the minimum propensity, making positivity indispensable.
Recently, \cite{kanxu25} investigates matrix completion in the context of matching markets, providing entrywise error guarantees for a nuclear-norm-regularized estimator.
Their framework accommodates dependent observations but also allows entries to be repeatedly observed.
This stands in contrast to our setting, where each entry is observed at most once.
\medskip
\noindent \textbf{Panel data.}
A growing body of work connects matrix completion with econometric panel data models, which often assume that the potential outcomes matrix is low rank \citep{amjad2018robust, amjad2019mrsc, athey2021matrix, tianyi_1, fernandez2020low, lihua25}.
This line of research is also closely related to work on factor models \citep{bai2019matrix, CAHAN2023113}.
\cite{amjad2018robust, amjad2019mrsc, athey2021matrix} provide an important contribution that connected matrix completion literature with panel data setting and establishing global error bounds under block-sparsity pattern.
\cite{Choi20092024} refine approach of \cite{athey2021matrix} by partitioning missing entries into smaller groups and applying convex estimation to each group, obtaining entrywise error bounds and asymptotic normality for certain statistics.
\cite{yan2024entrywiseinferencemissingpanel} build on this idea by developing a computationally efficient algorithm with improved entrywise inference for staggered adoption designs.
Our contribution is complementary to this line of work in scope rather than in dominance. Whereas panel-factor methods generally exploit global treatment-timing geometry, \texttt{SNN}~instead relies on a target-specific recovery condition.
Panel-specific methods may therefore yield sharper rates or stronger inference when their global geometry is credible.
Our framework instead targets settings where the observation geometry is too irregular to be captured by classical panel designs, but still contains enough local information for entrywise recovery.
\medskip
\noindent \textbf{Earlier versions.}
An earlier version of this work circulated as \cite{agarwal2021causalmatrixcompletion}, and a shorter extended abstract appeared at the {\em Conference on Learning Theory} in 2023. Those versions introduced the causal formulation, a preliminary \texttt{SNN}~estimator, and an oracle asymptotic normality result. The present article gives a substantially expanded treatment.
Methodologically, it introduces a simpler estimator for mean recovery, a heteroskedastic variance estimation procedure based on the same local principle, and a more tractable subroutine for finding informative local submatrices. Statistically, it develops a sharper asymptotic normality result, feasible studentized inference under heteroskedastic noise, and a theory for low-rank heteroskedastic variance estimation. It also situates \texttt{SNN}~more explicitly relative to subsequent work on panel-factor models.
\subsection{Paper Organization}
Section~\ref{sec:problem_setup} introduces our causal framework and formalizes the causal estimand.
Section~\ref{sec:estimator} presents the \texttt{SNN}~algorithm.
Section~\ref{sec:theoretical_results} establishes statistical guarantees for \texttt{SNN}~as an estimator of the mean matrix.
Section~\ref{sec:var.est} develops feasible entrywise inference and extends the analysis to heteroskedastic variance estimation.
Section~\ref{sec:sims} reports simulation studies.
Section~\ref{sec:conclusion} concludes.
\subsection{General Notation}
For a positive integer $n$, let $[n] = \{1, \dots, n\}$.
For index sets $\Ic \subseteq [m] $ and $\Jc \subseteq [n]$, let $\bM_{\Ic, \Jc}$ denote the $| \Ic | \times | \Jc |$ sub-matrix of $\bM \in \Rb^{m \times n}$ whose rows and columns are indexed by $\Ic$ and $\Jc$.
For a vector $\boldsymbol{v} \in \Rb^m$, define $\boldsymbol{v}_{\Ic}$ analogously.
We reserve $\boldsymbol{I}$ and $\boldsymbol{0}$ as the identity matrix and matrix/vector of zeros.
Let $\dagger$ denote the Moore-Penrose pseudoinverse and $\circ$ the entrywise (Hadamard) product.
For $p \in [1, \infty]$, let $\|\boldsymbol{v}\|_p$ denote the $\ell_p$-norm.
Denote the operator and Frobenius norms as $\| \bM \|_\text{op}$ and $\| \bM \|_F$.
The subgaussian and sub-exponential norms of a random vector $\boldsymbol{v}$ are denoted by $\| \boldsymbol{v} \|_{\psi_2}$ and $\|\boldsymbol{v}\|_{\psi_1}$.
Convergence in probability and distribution are denoted as $\xrightarrow{p}$ and $\rightsquigarrow$.
\section{A Causal Framework for Matrix Completion} \label{sec:problem_setup}
We adopt the potential outcomes framework \citep{neyman, rubin} to formulate matrix completion as a causal estimation problem.
Let $\widetilde{\boldsymbol{Y}} = [\widetilde{Y}_{ij}] \in \Rb^{m \times n}$ denote the matrix of potential outcomes, where $\widetilde{Y}_{ij}$ is the outcome that would be observed if unit $i \in [m]$ was exposed to treatment $j \in [n]$.
Let $\bD = [D_{ij}] \in \{0,1\}^{m \times n}$ denote the matrix of treatment assignments.
The observed outcomes matrix $\boldsymbol{Y} = [Y_{ij}] \in \{\Rb \cup \{\star\}\}^{m \times n}$, where $\star$ denotes a missing value, is defined entrywise as
\begin{align}\label{eq:sutva}
Y_{ij} \coloneqq \widetilde{Y}_{ij} \cdot \boldsymbol{1} \{D_{ij} = 1\} + \left( \star \right) \cdot \boldsymbol{1} \{D_{ij} = 0 \}.
\end{align}
The observation law of \eqref{eq:sutva}, known as the {\em stable unit treatment value assumption}, is a standard assumption in the causal inference literature that precludes interference \citep{imbens_rubin_2015}.
\subsection{Matrix Factor Model}
Without additional assumptions, the available observations offer no information about a missing entry.
To enable a faithful recovery of the underlying matrix, we impose a low-rank structure on the expected outcomes, as is standard in matrix completion.
\begin{assumption}[Low-rank factor model]\label{assump:lfm}
For every pair $(i,j)$, let $\widetilde{Y}_{ij} = \left\langle \boldsymbol{u}_i, \boldsymbol{v}_j \right \rangle + \varepsilon_{ij}$,
where $\boldsymbol{u}_i, \boldsymbol{v}_j \in \Rb^r$ are latent row and column factors, and $\varepsilon_{ij} \in \Rb$ is idiosyncratic noise.
\end{assumption}
With $\varepsilon_{ij}$ denoting (typically zero-mean) noise, $\langle \boldsymbol{u}_i, \boldsymbol{v}_j \rangle$ represents the underlying signal.
When $r \ll \min\{m, n\}$, the expected potential outcomes can be explained by significantly fewer factors than the dimensions of the matrix.
This renders the matrix recovery problem well-posed and identifiable under sensible conditions.
For these reasons, the signal-plus-noise factor model in Assumption~\ref{assump:lfm} is foundational in both matrix completion \citep{mc_survey} and econometrics factor models \citep{athey2021matrix}.
Empirical and theoretical work also supports the prevalence of approximate low-rank structure in data science applications \citep{udell2018big}, and smooth nonlinear latent variable models can often be well approximated by low-rank representations \citep{xu2017rates}.
\subsection{Causal Estimand} \label{sec:framework_causal_parameters}
Anchoring on Assumption~\ref{assump:lfm}, we define the causal estimand for entry $(i,j)$ as $$A_{ij} \coloneqq \mathbb{E}\left[\widetilde{Y}_{ij} \mid \boldsymbol{u}_i, \boldsymbol{v}_j \right],$$ the expected potential outcome for the $(i,j)$th entry, conditional on its latent factors.
The choice of $(i,j)$ is arbitrary: the analysis targets entrywise recovery, allowing the error associated with an individual missing entry to be studied directly.
Appendix~\ref{sec:add.params} of the Supplementary Material extends this framework to linear functionals of the mean matrix, including weighted averages and contrasts.
\subsubsection{Additional Notation}
Let $\bU \in \Rb^{m \times r}$ and $\bV \in \Rb^{n \times r}$ collect the latent row and column factors, and let $\bE = [\varepsilon_{ij}] \in \Rb^{m \times n}$ denote the noise matrix.
Define $\Ec = \{\bU, \bV, \bD\}$ as the information set, which we treat as fixed.
Further, define the neighborhood rows of column $j$ and neighborhood columns of row $i$ by $\textsf{NR}(j) = \{\ell \in [m]: D_{\ell j} = 1\}$ and $\textsf{NC}(i) = \{q \in [n]: D_{iq} = 1\}$.
We call $\textsf{AR} \subseteq \textsf{NR}(j)$ and $\textsf{AC} \subseteq \textsf{NC}(i)$ a set of {\em anchor rows} and {\em anchor columns} for the pair $(i,j)$ if $D_{\ell q} = 1$ for all $(\ell, q) \in \textsf{AR} \times \textsf{AC}$.
Thus, anchor rows are observed at the target column, anchor columns are observed for the target row, and the anchor submatrix $\normalfont{\textsf{AR}} \times \normalfont{\textsf{AC}}$ is fully observed.
Together, they form a local observed cross surrounding the missing entry $(i,j)$.
For a visualization, see Figure~\ref{fig:obs_pattern}.
\begin{figure} [!t]
\centering
\includegraphics[width=0.65\linewidth]
{images/obs_pattern.png}
\caption{Visualization of the observed matrix $\boldsymbol{Y}$, highlighting the neighborhood sets $(\textsf{NR}(j), \textsf{NC}(i))$ and anchor sets $(\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}})$ associated with the target $(i,j)$ pair. Unobserved entries are depicted in white.}
\label{fig:obs_pattern}
\end{figure}
\subsubsection{Assumptions for Estimation}
We state assumptions to estimate $A_{ij}$ from observed data.
\begin{assumption}[Selection on latent factors] \label{assump:mean_ind}
Let $\mathbb{E}[ \bE \mid \Ec] = \boldsymbol{0}$.
\end{assumption}
By the tower law, Assumption~\ref{assump:mean_ind} implies $\mathbb{E}[\bE \mid \bU, \bV] = \boldsymbol{0}$.
Coupled with Assumption~\ref{assump:lfm}, this implies $\bA = \mathbb{E}[\widetilde{\boldsymbol{Y}} \mid \bU, \bV]$ has rank $r$, with $r$ typically much smaller than $m$ and $n$.
Assumptions~\ref{assump:lfm} and \ref{assump:mean_ind} also jointly imply $\mathbb{E}[\widetilde{\boldsymbol{Y}} \mid \bU, \bV] = \mathbb{E}[\widetilde{\boldsymbol{Y}} \mid \Ec]$, so the potential outcomes are mean independent of the treatment assignments, conditioned on the latent factors.
Hence, the latent factors serve as latent confounders, playing an analogous role to observed covariates under the
classical selection on observables assumption.
This form of latent ignorability is conceptually aligned with assumptions explored in panel data models \citep{athey2021matrix, avi_eli_21, agarwal2021synthetic} and related contexts \citep{kallus2018causal}.
As with any assumption concerning unobserved confounding, Assumption~\ref{assump:mean_ind} is fundamentally untestable.
Its justification relies on domain expertise and context-specific reasoning about the data-generating process.
\begin{assumption}[Linear row span inclusion] \label{assump:linear_span}
Conditioned on $\Ec$ and given anchor rows $\normalfont{\textsf{AR}}$, we have $\boldsymbol{u}_i \in {\normalfont \text{span}}(\{ \boldsymbol{u}_\ell: \ell \in \normalfont{\textsf{AR}} \})$.
\end{assumption}
Assumption~\ref{assump:linear_span} requires the latent factor of the target row to be expressible as a linear combination of the latent factors associated with the anchor rows.
This condition arises naturally under the low-rank structure in Assumption~\ref{assump:lfm}. If $\text{span}(\{\boldsymbol{u}_\ell: \ell \in \textsf{AR}\}) = \Rb^r$, then Assumption~\ref{assump:linear_span} follows immediately because $\boldsymbol{u}_i \in \Rb^r$.
In practical terms, enough rows must be observed at column $j$ so that their latent factors span a rich enough subspace to include the target row factor. Although this condition cannot be verified directly because the factors are latent, its plausibility can be assessed empirically through training error diagnostics.
\subsubsection{Estimation Result}
With our assumptions in hand, we next establish that our causal estimand $A_{ij}$ can be expressed using quantities that are, in principle, estimable from observed data.
\begin{proposition}\label{thm:identification}
Let Assumptions~\ref{assump:lfm}--\ref{assump:linear_span} hold.
Then, there exists a $\boldsymbol{\beta} \in \Rb^{|\normalfont{\textsf{AR}} |}$ such that (a) $A_{ij} = \sum_{\ell \in \normalfont{\textsf{AR}}} \beta_\ell \cdot \mathbb{E}[Y_{\ell j} \mid \Ec]$, and (b) $\mathbb{E}[Y_{iq} \mid \mathcal{E}] = \sum_{\ell \in \normalfont{\textsf{AR}}} \beta_\ell \cdot \mathbb{E}[ Y_{\ell q} \mid \Ec]$ for every $q \in \normalfont{\textsf{AC}}$.
\end{proposition}
The first display states that the target estimand $A_{ij}$ can be recovered from the anchor rows once the coefficient vector $\boldsymbol{\beta}$ is known. The second display suggests a natural way to estimate $\boldsymbol{\beta}$: regress the observed entries of the target row over the anchor columns on the corresponding observed entries of the anchor rows. Section~\ref{sec:estimator} develops this idea into a concrete algorithm.
This framework accommodates a broad class of MNAR mechanisms. In particular, the assignment matrix $\bD$ may depend arbitrarily on the latent factors $\bU$ and $\bV$. Since $\bU \bV^\top$ is the signal matrix of expected potential outcomes, each observation indicator $D_{ij}$
may be a random or deterministic function of the underlying value $A_{ij}$.
Moreover, the indicators in $\bD$ may be dependent across entries, and positivity is not required.
The flexibility of the observation mechanism is supplied by the factor model assumptions. The framework is therefore outcome-model-based, rather than design-based: identification follows from the structure of the potential outcomes and the availability of an informative local anchor cross, not from probabilistic restrictions on the assignment matrix.
\section{An Algorithm for Matrix Completion with MNAR Data} \label{sec:estimator}
We introduce {\em synthetic nearest neighbors} (\texttt{SNN}), an algorithm for matrix completion under MNAR data.
Motivated by Proposition~\ref{thm:identification}, \texttt{SNN}~takes the following steps to estimate the $(i,j)$th entry.
\begin{enumerate}
\item {\bf Anchor set discovery}: Obtain anchor sets $(\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}})$ such that the anchor submatrix $\boldsymbol{Y}_{\textsf{AR}, \textsf{AC}} \coloneqq [Y_{\ell q}: (\ell, q) \in \textsf{AR} \times \textsf{AC}] \in \Rb^{|\textsf{AR}| \times |\textsf{AC}|}$, target-row vector $\boldsymbol{Y}_{i, \textsf{AC}} \coloneqq [Y_{iq}: q \in \textsf{AC}] \in \Rb^{|\textsf{AC}|}$, and target-column vector $\boldsymbol{Y}_{\normalfont{\textsf{AR}}, j} \coloneqq [Y_{qj}: q \in \textsf{AC}] \in \Rb^{|\textsf{AC}|}$ are all fully observed.
\item {\bf Spectral denoising}: Let the singular value decomposition of the anchor submatrix be $\boldsymbol{Y}_{\textsf{AR}, \textsf{AC}} = \sum_{\ell \ge 1} \widehat{s}_\ell \widehat{\boldsymbol{u}}_\ell \widehat{\boldsymbol{v}}_\ell^\top$, where $\widehat{s}_\ell \in \Rb$ are the singular values in decreasing order, and $\widehat{\boldsymbol{u}}_\ell \in \Rb^{|\textsf{AR}|}$, $\widehat{\boldsymbol{v}}_\ell \in \Rb^{|\textsf{AC}|}$ are the corresponding left and right singular vectors.
Given a spectral threshold $t \le \min\{|\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}|\}$, define the rank-$t$ approximation by $\boldsymbol{Y}^{(t)}_{\textsf{AR}, \textsf{AC}} \coloneqq \sum_{\ell = 1}^{t} \widehat{s}_\ell \widehat{\boldsymbol{u}}_\ell \widehat{\boldsymbol{v}}_\ell^\top$.
\item {\bf Parameter estimation}: Estimate the coefficients by principal component regression (\texttt{PCR})
\begin{align}
\widehat{\boldsymbol{\beta}}&\in \operatorname*{arg\,min}_{\boldsymbol{\beta} \in \Rb^{|\textsf{AR}|}} ~\left\| \boldsymbol{Y}_{i, \textsf{AC}} - \left( \boldsymbol{Y}^{(t)}_{\textsf{AR}, \textsf{AC}} \right)^\top \boldsymbol{\beta} \right\|_2^2
= \left( \boldsymbol{Y}^{(t)}_{\textsf{AR}, \textsf{AC}} \right)^{\top, \dagger} \cdot \boldsymbol{Y}_{i, \textsf{AC}}.
\label{eq:beta.hat}
\end{align}
\item {\bf Point estimation}: Estimate the target entry by $\hA_{ij} \coloneqq \sum_{\ell \in \textsf{AR}} \widehat{\beta}_\ell \cdot Y_{\ell j} = \left \langle \boldsymbol{Y}_{\textsf{AR}, j}, \widehat{\boldsymbol{\beta}} \right \rangle$.
\item {\bf Clipping:} If $A_{ij} \in [a, b]$ is known, define
$\hA^{\text{clip}}_{ij} \coloneqq \min\left\{ b, \max\left\{ a, \hA_{ij} \right\} \right\}$.
\end{enumerate}
\subsection{Interpretation} \label{sec:snn.intuition}
As illustrated by Figure~\ref{fig:obs_pattern}, \texttt{SNN}~reduces the original $m \times n$ completion problem to the local array
\begin{align}
\begin{pmatrix}
&\boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}~ & \boldsymbol{Y}_{\textsf{AR}, j}
\\
& \boldsymbol{Y}^\top_{i, \textsf{AC}}~ & \star
\end{pmatrix}.
\end{align}
Within this local submatrix, all entries are observed except the target entry $(i,j)$. These observed entries provide the information used to estimate $A_{ij}$.
\texttt{SNN}, however, does not operate directly on the raw entries.
Under Assumption~\ref{assump:lfm}, each outcome contains both signal and idiosyncratic noise, whereas Assumption~\ref{assump:mean_ind} implies that the conditional mean of the anchor block is low rank. The leading singular directions of $\boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}$ are therefore interpreted as carrying signal, while the smaller singular values primarily reflect noise. Accordingly, \texttt{SNN}~performs spectral denoising to extract the dominant components, $\boldsymbol{Y}^{(t)}_{\textsf{AR}, \textsf{AC}}$, which serves as an estimate of the low-rank signal $\mathbb{E}[\boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}} \mid \Ec]$.
After denoising, \texttt{SNN}~follows the intuition of nearest-neighbor collaborative filtering \citep{goldberg1992using, linden2003amazon}, but with a key modification. Standard nearest-neighbor methods search for one or more rows that are individually close to row $i$. \texttt{SNN}~instead constructs a synthetic neighbor for row $i$ by reweighting the rows in $\textsf{AR}$ according to $\widehat{\boldsymbol{\beta}}$. Thus, \texttt{SNN}~does not require any single anchor row to resemble the target row. Rather, it only requires that a weighted combination of anchor rows approximate the target row.
From this perspective, \texttt{SNN}~is conceptually aligned with synthetic controls \citep{abadie1, abadie2}. In its canonical form, synthetic controls constrains the elements of $\widehat{\boldsymbol{\beta}}$ to be nonnegative and sum to one. By contrast, \texttt{SNN}~constrains the fitted weights through \texttt{PCR}: $\widehat{\boldsymbol{\beta}}$ lies in the linear subspace spanned by the top $t$ left singular vectors of the anchor submatrix.
This restriction is naturally adapted to the low-rank factor model.
In short, \texttt{SNN}~localizes estimation to an observed anchor neighborhood around a target entry, as in nearest-neighbor methods, but replaces uniform or distance-based averaging with PCR-based synthetic-control weighting.
Although \texttt{SNN}~is written in row-side form, the same point estimate admits an equivalent column-side representation. Define the column-side \texttt{PCR}~estimate by
\begin{align}
\widehat{\boldsymbol{\alpha}} \coloneqq \left(\boldsymbol{Y}^{(t)}_{\normalfont{\textsf{AR}},\normalfont{\textsf{AC}}}\right)^\dagger \cdot \boldsymbol{Y}_{\normalfont{\textsf{AR}}, j}.
\label{eq:alpha.hat}
\end{align}
By \cite[Theorem 1 and Corollary 1]{sameroot}, $\hA_{ij} \coloneqq \big\langle \boldsymbol{Y}_{\textsf{AR}, j}, \widehat{\boldsymbol{\beta}} \big\rangle = \big\langle \boldsymbol{Y}_{i, \textsf{AC}}, \widehat{\boldsymbol{\alpha}} \big\rangle = \big\langle \widehat{\boldsymbol{\alpha}}, \big(\boldsymbol{Y}^{(t)}_{\normalfont{\textsf{AR}},\normalfont{\textsf{AC}}}\big)^\top \widehat{\boldsymbol{\beta}} \big\rangle$.
These identities give three equivalent interpretations of \texttt{SNN}: a row-side synthesis, a column-side synthesis, and a two-sided synthesis that uses both simultaneously through the denoised anchor block.
This symmetry recurs throughout the theory. Assumption~\ref{assump:subspace} in Section~\ref{sec:additional_assumptions} formalizes the column-side span condition, analogous to Assumption~\ref{assump:linear_span}; the inferential theory in Section~\ref{sec:normality} shows how both sides enter the first-order uncertainty of $\hA_{ij}$; and Section~\ref{sec:var.est} turns the resulting decomposition into feasible studentization.
For concreteness, we retain the row-side notation.
\subsection{Practical Considerations} \label{sec:fine.print}
\noindent \textbf{Discovering anchor sets.}
The first step of \texttt{SNN}~is to identify a valid anchor pair $(\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}})$.
Such pairs arise naturally in panel applications with block-sparse observation patterns, but similar structures also occur in less regular settings.
For instance, \cite{ma2019missing} documented pronounced block patterns in recommender systems datasets, such as the \texttt{Coat} \citep{schnabelfwang16} and \texttt{MovieLens-100k} \citep{movielens}.
Spectral biclustering \citep{spectral_bicluster} offers one computationally tractable way to locate coherent row-column blocks. Alternatively, the observation mask can be viewed as the incidence matrix of a bipartite graph, in which every valid anchor pair corresponds to a biclique.
Appendix~\ref{sec:max_biclique} in the Supplementary Material formalizes this reduction and reviews algorithms for finding maximal bicliques. We use these approaches in the examples of
Section~\ref{sec:teaser}. Corollary~\ref{cor:mcar} further shows that valid anchor pairs also arise with high probability under MCAR sampling when the observation probability is sufficiently large.
\medskip \noindent \textbf{Choosing the spectral threshold.}
The parameter $t$ controls the degree of spectral regularization. Choosing $t$ too small may discard relevant signal; choosing it too large may reintroduce noise.
A simple diagnostic is to inspect the singular-value spectrum and select $t$ near its elbow. A more data-driven alternative is cross-validation on held-out observed entries. One may also use a universal singular-value cutoff, retaining only those components above a prescribed threshold, as in \cite{donoho14} and \cite{Chatterjee15}.
\medskip \noindent \textbf{Ensemble \texttt{SNN}.} Appendix~\ref{sec:ensemble} of the Supplementary Material introduces an ensemble extension that constructs synthetic neighbors from multiple anchor crosses and aggregates the resulting estimates.
\section{Entrywise Guarantees for Mean Recovery} \label{sec:theoretical_results}
This section establishes statistical guarantees for \texttt{SNN}~as an estimator of the target mean $A_{ij}$.
Throughout, for any random object $\bW$, write $\bar{\bW} \coloneqq \mathbb{E}[\bW \mid \Ec]$.
\subsection{Additional Assumptions for Mean Recovery} \label{sec:additional_assumptions}
We state additional assumptions utilized in our analysis.
\begin{assumption}[Subgaussian noise] \label{assump:subg}
Conditioned on $\Ec$, $\varepsilon_{ij}$ are independent subgaussian random variables satisfying $\operatorname{Var}(\varepsilon_{ij})= \sigma_{ij}^2 \le \sigma^2$ and $\| \varepsilon_{ij} \|_{\psi_2} \le C_\varepsilon \sigma_{ij}$ for constants $C_\varepsilon, \sigma > 0$.
\end{assumption}
While the latent row and column factors may be arbitrarily correlated, Assumption~\ref{assump:subg} imposes independence and subgaussian tail behavior on the noise. Although this is a standard assumption in the matrix completion literature, it may be restrictive in certain applications. A formal treatment under more general noise models remains an important direction for future research.
\begin{assumption}[Bounded signal] \label{assump:bounded}
Let $\langle \boldsymbol{u}_i, \boldsymbol{v}_j \rangle \in [-1,1]$.
\end{assumption}
The bound $[-1,1]$ can be extended to $[a,b]$ for any $a, b \in \Rb$ with $a \le b$.
\begin{assumption}[Well-balanced spectra]\label{assump:spectra}
Given anchor sets $(\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}})$, the condition number $\kappa$ of $\bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}$ satisfies $\kappa^{-1} \ge c_y$, and $\| \bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}} \|_F^2 \ge c_y' \cdot |\normalfont{\textsf{AR}}| \cdot |\normalfont{\textsf{AC}}|$ for constants $c_y, c_y' > 0$.
\end{assumption}
Assumption~\ref{assump:spectra} requires the nonzero singular values of the anchor signal matrix to be neither too small nor too uneven. While potentially restrictive, it ensures the low-rank structure is visible above the noise floor. This condition parallels familiar pervasiveness, incoherence, and beta-min requirements from factor models, high-dimensional inference, and matrix completion \citep{chamberlainfactor, beta_min, fan2018eigenvector, bai_matrix, cai2021nonconvex}. Empirically, it can be assessed through the singular-value profile of $\boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}$, as discussed in Section~\ref{sec:fine.print}.
\begin{assumption}[Linear column span inclusion]\label{assump:subspace}
Conditioned on $\Ec$ and given anchor columns $\normalfont{\textsf{AC}}$, we have $\boldsymbol{v}_j \in {\normalfont \text{span}}(\{ \boldsymbol{v}_q: q \in \normalfont{\textsf{AC}} \})$.
\end{assumption}
Assumption~\ref{assump:subspace} is the column-side counterpart of Assumption~\ref{assump:linear_span}. Assumption~\ref{assump:linear_span} ensures that the target row $i$ can be synthesized from the anchor rows $\normalfont{\textsf{AR}}$; Assumption~\ref{assump:subspace} ensures that the target column $j$ can be extrapolated from the anchor columns $\normalfont{\textsf{AC}}$. Together, they make the local anchor cross informative enough to identify the target entry.
To see why Assumption~\ref{assump:subspace} is needed, recall that \texttt{SNN}~estimates the row-side weights by regressing $\boldsymbol{Y}_{i, \normalfont{\textsf{AC}}}$ on the denoised anchor block $\big(\boldsymbol{Y}^{(t)}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}\big)^\top$. If the population anchor block $\bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}$ is rank deficient, the coefficients $\boldsymbol{\beta}$ satisfying the row representation in Proposition~\ref{thm:identification} are not unique. \texttt{PCR}~instead identifies the projection of any valid representer onto the column space of $\bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}$, denoted as
\begin{align}
\boldsymbol{\beta}^* \coloneqq \bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}^{\top, \dagger} \bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}^\top \cdot \boldsymbol{\beta} = \bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}^{\top, \dagger} \bar{\boldsymbol{Y}}_{i, \normalfont{\textsf{AC}}}.
\label{eq:beta.star}
\end{align}
Assumption~\ref{assump:subspace} ensures that this projection preserves the prediction at the target column. Accordingly, Assumption~\ref{assump:linear_span} provides a representation of the target row, while Assumption~\ref{assump:subspace} validates its extrapolation.
Viewed more broadly, Assumption~\ref{assump:subspace} is a local generalization condition: the held-out factor $\boldsymbol{v}_j$ must lie in the span of the training-column factors $\{\boldsymbol{v}_q : q \in \textsf{AC}\}$.
This is a linear-algebraic analogue of covariate similarity assumptions in statistical learning.
\subsection{Non-asymptotic Error}
Below, we use $\lesssim$ to reflect inequalities that ignore constants factors independent of our model parameters.
For convenience, we denote $|\textsf{ARC}| \coloneqq |\textsf{AR}| \cdot |\textsf{AC}|$ and $L \coloneqq \log(|\normalfont{\textsf{ARC}}|)$.
\begin{theorem} \label{thm:consistency}
Let Assumptions~\ref{assump:lfm}--\ref{assump:subspace} hold with $t = \bar{r} \coloneqq \emph{rank}(\mathbb{E}[\boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}\mid\Ec])$. Suppose that
\begin{align}
|\normalfont{\textsf{ARC}} | \ge C \sigma^2 \bar{r} \left( \sqrt{|\normalfont{\textsf{AR}}|} + \sqrt{|\normalfont{\textsf{AC}}|} + \sqrt{L} \right)^2 \label{eq:rank.sep}
\end{align}
and $c \sigma^2 L \le \min\{|\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}|\}$ for sufficiently large constants $C, c > 0$.
Then, conditional on $\Ec$, with probability at least $1 - \mathcal{O}(|\normalfont{\textsf{ARC}}|^{-10})$,
\begin{align}
\left| \hA_{ij} - A_{ij} \right| &\lesssim \Phi \coloneqq \frac{\sigma \bar{r}}{\min \left\{ \sqrt{|\normalfont{\textsf{AR}}|}, \sqrt{|\normalfont{\textsf{AC}}|} \right\} }
+ \frac{\sqrt{\sigma \bar{r}} L^{1/4}}{\min\left\{ |\normalfont{\textsf{AC}}|, |\normalfont{\textsf{AR}}|^2 / \left(\sigma^2 L \right) \right\}^{1/4}}.
\label{eq:thm.hp.special}
\end{align}
\end{theorem}
\noindent \textbf{Rate implications.}
Theorem~\ref{thm:consistency} gives a non-asymptotic entrywise error bound that suggests \texttt{SNN}~is consistent whenever $\Phi = o(1)$.
In the balanced regime $|\normalfont{\textsf{AR}}| \asymp |\normalfont{\textsf{AC}}|$, the bound reduces to $\Phi = \Oc\left(\bar{r} |\normalfont{\textsf{AR}}|^{-1/2} + \sqrt{\bar{r}}L^{1/4} |\normalfont{\textsf{AR}}|^{-1/4} \right)$. Thus, a sufficient effective-rank condition is $\bar{r} = o(\sqrt{|\normalfont{\textsf{AR}}| / L} )$. In the favorable asymmetric regime $|\normalfont{\textsf{AC}}| \asymp |\normalfont{\textsf{AR}}|^2$, the rate becomes $\Phi = \Oc\left( |\normalfont{\textsf{AR}}|^{-1/2} \cdot (\bar{r} + \sqrt{\bar{r} L}) \right)$, yielding $\bar{r} = \tilde{o}(\sqrt{|\normalfont{\textsf{AR}}|})$. Hence, for fixed $\bar{r}$, \texttt{SNN}~approaches the near-parametric $|\normalfont{\textsf{AR}}|^{-1/2}$ scale.
The apparent advantage of the asymmetric regime should not be interpreted as an intrinsic benefit of having more anchor columns than anchor rows. It arises from the current perturbation analysis, which controls the row-side \texttt{PCR}~error using conservative operator-norm bounds. Since \texttt{SNN}~is symmetric in rows and columns, as highlighted in Section~\ref{sec:snn.intuition}, a sharper analysis may yield more symmetric sufficient conditions.
The auxiliary condition $c \sigma^2 L \le \min\{|\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}|\}$ is primarily used to simplify the displayed rate. The more substantive requirement is the signal-to-noise condition in \eqref{eq:rank.sep}, which ensures that the smallest nonzero singular value of the population anchor block dominates the operator-norm fluctuation of the noisy anchor block.
Finally, Theorem~\ref{thm:consistency} is stated for the oracle choice $t = \bar{r}$.
A formal analysis under misspecified rank is left for future work, although existing results suggest that overestimating the rank is generally less damaging than underestimating it \citep{pcr_aos}.
\medskip \noindent \textbf{Connection to MCAR.}
The next result specializes Theorem~\ref{thm:consistency} to the classical MCAR setting.
\begin{corollary} \label{cor:mcar}
Fix a target entry $(i,j)$, $\delta > 0$, and $0 < \eta < 1$. Let $m = n$, and let $d = o\big(\log(n)\big)$ with $d \rightarrow \infty$. Consider the MCAR regime with the observation probability $p$ satisfying
\begin{align}
p \ge \left[\frac{C\log(2/\eta)}{\binom{n-1}{d}^2}\right]^{1/[d(d+2)]}
\end{align}
for a sufficiently large constant $C > 0$.
Then, for all sufficiently large $n$, with probability at least $1 - \eta/2$, a valid $d \times d$ anchor pair $(\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}})$ exists.
Suppose further that a mask-measurable selection of such an anchor cross satisfies the conditions of Theorem~\ref{thm:consistency} and $d \ge c \sigma^2 \bar{r}^2 \max\{\delta^{-2}, \delta^{-4} \log(d)\}$ for a sufficiently large $c > 0$. Then, over both the MCAR mask and outcome noise,
\begin{align}
\Pb\left( \left|\hA_{ij}- A_{ij} \right| \le \delta \mid \bU, \bV \right) \ge 1 - \eta.
\end{align}
\end{corollary}
Corollary~\ref{cor:mcar} shows that valid anchor crosses arise even under fully random observation patterns, provided the observation probability is large enough (though still scaling as $o(1)$). This is not intended as a sharp MCAR guarantee. Rather, it serves as a bridge: the same local-anchor logic that drives \texttt{SNN}~under MNAR also recovers the familiar MCAR setting as a special case.
The comparison also clarifies the price of robustness. Suppressing logarithmic factors, classical MCAR matrix completion can exploit global random sampling to obtain a Frobenius error $(1/n) \| \bhA - \bA \|_F = \tilde{\Oc}(\delta)$ with $\tilde{\Oc}(nr \delta^{-2})$ samples.
Under MNAR sampling, such global sample-complexity statements are generally harder to formulate.
A more transparent comparison arises when every target can use a shared anchor pair $(\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}})$, as in Figure~\ref{fig:Lshape}. In a recommender system, the anchor rows may represent highly active viewers or critics, while the anchor columns may represent widely viewed canonical films, such as those used during onboarding.
Up to logarithmic factors, uniform versions of Theorem~\ref{thm:consistency} give $\Oc(\delta)$ entrywise recovery when $|\normalfont{\textsf{AR}}| = \tilde{\Omega}(\bar{r}^2 \delta^{-2})$ and $|\normalfont{\textsf{AC}}| = \tilde{\Omega}(\bar{r}^2 \delta^{-4})$ using $\tilde{\Oc}(n \bar{r}^2 \delta^{-4})$ samples.
Thus, when $\bar{r} \asymp r$, \texttt{SNN}~pays an additional factor of order $r \delta^{-2}$ relative to the MCAR benchmark, but provides entrywise guarantees under substantially more general missingness mechanisms.
A sharper account of the trade-offs between MCAR-tailored and MNAR-robust procedures is left for future work.
\begin{figure} [!t]
\centering
\includegraphics[width=0.3\linewidth]
{images/Lshape.png}
\caption{Proposed observation pattern to achieve entrywise matrix recovery for \texttt{SNN}.}
\label{fig:Lshape}
\end{figure}
\subsection{Asymptotic Normality} \label{sec:normality}
We next refine the finite-sample guarantee of Theorem~\ref{thm:consistency} into an asymptotic distributional statement.
Recall the row-side population representer $\boldsymbol{\beta}^*$ defined in \eqref{eq:beta.star}, and define its column counterpart by $\boldsymbol{\alpha}^* \coloneqq \bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}},\normalfont{\textsf{AC}}}^\dagger \bar{\boldsymbol{Y}}_{\normalfont{\textsf{AR}}, j}$.
To facilitate our presentation below, define
\begin{align}
\upsilon^2 \coloneqq \sum_{\ell \in \normalfont{\textsf{AR}}} \left(\beta^*_\ell\right)^2 \sigma^2_{\ell j}
+ \sum_{q \in \normalfont{\textsf{AC}}} \left(\alpha^*_q\right)^2 \sigma^2_{iq}
+ \sum_{\ell \in \normalfont{\textsf{AR}}} \sum_{q \in \normalfont{\textsf{AC}}} \left(\beta^*_\ell \alpha^*_q\right)^2 \sigma^2_{\ell q}.
\label{eq:var.asymp}
\end{align}
\begin{theorem} \label{thm:normality}
Let the setup of Theorem~\ref{thm:consistency} hold.
Conditional on $\Ec$, $\hA_{ij} - A_{ij} = G + R$.
If Assumption~\ref{assump:subg} is specialized to Gaussian noise, then $G / \upsilon \sim \mathcal{N}(0,1)$; for general subgaussian noise, if
\begin{align}
\frac{\max\left\{ \max_{\ell \in \normalfont{\textsf{AR}}} \left| \beta^*_\ell \right| \sigma_{\ell j}, ~\max_{q \in \normalfont{\textsf{AC}}} \left| \alpha^*_q \right| \sigma_{iq}, ~\max_{(\ell, q) \in \normalfont{\textsf{AR}} \times \normalfont{\textsf{AC}}} \left| \beta^*_\ell \alpha^*_q \right| \sigma_{\ell q} \right\}}{\upsilon} = o(1),
\label{eq:clt}
\end{align}
then $G / \upsilon \rightsquigarrow \mathcal{N}(0,1)$.
Moreover, with probability at least $1 - \Oc(|\normalfont{\textsf{ARC}}|^{-10})$,
\begin{align}
\left| R \right| \lesssim
\Psi \coloneqq
\sigma \Phi \cdot \left( \frac{ \max\left\{\sqrt{\bar{r}}, \sqrt{L}\right\}}{\min \left\{ \sqrt{|\normalfont{\textsf{AR}}|}, \sqrt{|\normalfont{\textsf{AC}}|} \right\}} + \frac{ \sqrt{\bar{r} L}}{\sqrt{|\normalfont{\textsf{ARC}}|}}
\right),
\label{eq:psi}
\end{align}
where $\Phi$ is defined in \eqref{eq:thm.hp.special}.
Consequently, if $\Psi / \upsilon = o(1)$, then as $|\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}| \rightarrow \infty$,
\begin{align}
\frac{\hA_{ij} - A_{ij}}{\upsilon} \rightsquigarrow \mathcal{N}(0,1).
\end{align}
\end{theorem}
\noindent \textbf{Decomposition interpretation.} The variance decomposition in \eqref{eq:var.asymp} makes the row-column symmetry of \texttt{SNN}~statistically consequential. Its first term reflects noise in the target column, weighted by the row representer $\boldsymbol{\beta}^*$; its second reflects noise in the target row, weighted by the column representer $\boldsymbol{\alpha}^*$; and its third reflects noise in the anchor block, weighted by their products. Thus, even though the point estimate $\hA_{ij}$ may be written from either side, its first-order uncertainty depends on both.
The remainder envelope $\Psi$ in \eqref{eq:psi} collects the higher-order effects of estimating the two representers and the anchor subspace.
\medskip \noindent \textbf{Rate implications.}
Two ingredients are needed for asymptotic normality.
First, the leading fluctuation must satisfy a central limit theorem; condition \eqref{eq:clt} is a Lyapunov-type notion that prevents any single weighted noise entry from dominating the leading term. Second, the remainder must vanish relative to the standard-error scale so that $\Psi / \upsilon = o(1)$.
For illustration, suppose the relevant variances are bounded below and the population representers are nondegenerate: $\| \boldsymbol{\beta}^* \|_2 \gtrsim |\normalfont{\textsf{AR}}|^{-1/2}$ and $\| \boldsymbol{\alpha}^* \|_2 \gtrsim |\normalfont{\textsf{AC}}|^{-1/2}$. In the balanced regime $|\normalfont{\textsf{AR}}| \asymp |\normalfont{\textsf{AC}}|$, $\Psi / \upsilon = \Oc\left( (\sqrt{\bar{r}} + \sqrt{L}) \Phi \right)$, yielding the condition $\bar{r} = o(|\normalfont{\textsf{AR}}|^{1/4} L^{-1/4} )$.
Thus, inference imposes a stricter effective-rank requirement than consistency: the point estimator must not only converge, but admit a first-order linear approximation on the standard-error scale. The product structure of $\Psi$---first-stage error $\Phi$ multiplied by a perturbation factor---makes this distinction explicit.
\subsection{Interpretation and Scope}
The guarantees in Theorems~\ref{thm:consistency} and \ref{thm:normality} are inherently local. Recovery of $A_{ij}$ depends on the existence of an informative anchor cross around that entry---not on positivity, independent sampling, staggered adoption, or any prescribed geometry of the full observation mask. Outside the anchor structure, the missingness pattern may be deterministic, dependent, and irregular, subject to our factor model.
This locality defines both the strength and the limitation of \texttt{SNN}. The method replaces global design assumptions with local outcome-model conditions: suitable row and column spans and a sufficiently strong anchor signal. When the observation pattern has exploitable global structure, procedures tailored to that structure may be sharper. When the global pattern is irregular but contains informative local crosses, \texttt{SNN}~provides an entrywise route to recovery and inference.
\section{Feasible Entrywise Inference} \label{sec:var.est}
Theorem~\ref{thm:normality} gives an oracle normal approximation for \texttt{SNN}~in terms of the unknown leading standard error $\upsilon$. We now convert that result into a feasible entrywise inference procedure. This requires two pieces: an estimator of $\upsilon^2$, and estimators of the heteroskedastic noise variances that enter $\upsilon^2$.
\subsection{Plug-in Studentization}
Recall the \texttt{PCR}~estimates $(\widehat{\boldsymbol{\beta}}, \widehat{\boldsymbol{\alpha}})$ from \eqref{eq:beta.hat} and \eqref{eq:alpha.hat}.
We estimate $\upsilon^2$ using the natural plug-in estimator
\begin{align}
\widehat{\upsilon}^2 \coloneqq \sum_{\ell \in \normalfont{\textsf{AR}}} \widehat{\beta}_\ell^2 \widehat{\sigma}^2_{\ell j}
+ \sum_{q \in \normalfont{\textsf{AC}}} \widehat{\alpha}_q^2 \widehat{\sigma}^2_{iq}
+ \sum_{\ell \in \normalfont{\textsf{AR}}} \sum_{q \in \normalfont{\textsf{AC}}} \widehat{\beta}^2_\ell \widehat{\alpha}_q^2 \widehat{\sigma}^2_{\ell q},
\label{eq:hvar.asymp}
\end{align}
where each $\widehat{\sigma}^2$ estimates the corresponding noise variance and is constrained to lie in $[0, \sigma^2_+]$.
\begin{theorem} \label{thm:var.est.asymp}
Let the setup of Theorem~\ref{thm:normality} hold. Assume further that there exists a deterministic sequence $\Lambda_\sigma$ such that, with probability at least $1-p_\sigma$,
\begin{align}
\max\left\{
\max_{\ell \in \normalfont{\textsf{AR}}} \left| \widehat{\sigma}^2_{\ell j} - \sigma^2_{\ell j} \right|,
\max_{q \in \normalfont{\textsf{AC}}} \left| \widehat{\sigma}^2_{iq} - \sigma^2_{iq} \right|,
\max_{\ell \in \normalfont{\textsf{AR}}, q \in \normalfont{\textsf{AC}}} \left| \widehat{\sigma}^2_{\ell q} - \sigma^2_{\ell q} \right|
\right\} \le \Lambda_\sigma.
\label{eq:var.bound}
\end{align}
Then, conditional on $\Ec$, with probability at least $1 - p_\sigma - \Oc(|\normalfont{\textsf{ARC}}|^{-10})$,
\begin{align}
\left|\widehat{\upsilon}^2 - \upsilon^2 \right|
\lesssim
\Gamma
&\coloneqq
\left(\sigma^2_+ + \Lambda_\sigma \right)
\left\{ \frac{\bar{r}\Phi + \Phi^2}{\min\left\{|\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}| \right\}}
\left( 1 + \frac{\bar{r} + \Phi^2}{\min\left\{|\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}|\right\}} \right)
\right\}
+ \frac{\bar{r} \Lambda_\sigma}{\min\{|\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}|\}}.
\label{eq:gamma}
\end{align}
Consequently, if $\Gamma / \upsilon^2 = o(1)$, then as $|\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}| \rightarrow \infty$, $\widehat{\upsilon} / \upsilon \xrightarrow{p} 1$ and thus,
\begin{align}
\frac{\hA_{ij} - A_{ij}}{\widehat{\upsilon}} \rightsquigarrow \mathcal{N}(0,1).
\end{align}
\end{theorem}
\noindent \textbf{Rate implications.}
Theorem~\ref{thm:var.est.asymp} states that $\widehat{\upsilon}^2$ consistently estimates $\upsilon^2$, and therefore permits feasible studentization, once two requirements are met. First, the local variance estimates must satisfy the uniform bound in \eqref{eq:var.bound} over the variance entries that appear in $\upsilon^2$. Second, the variance-perturbation envelope $\Gamma$ must be negligible relative to $\upsilon^2$.
The envelope $\Gamma$ has two sources. The term multiplied by $(\sigma^2_+ + \Lambda_\sigma)$ captures the effect of estimating the \texttt{PCR}~row and column weights. The final term measures the direct contribution of estimating the local heteroskedastic variances.
To interpret $\Gamma / \upsilon^2 = o(1)$, consider the same nondegenerate setting used after Theorem~\ref{thm:normality}. In the balanced regime $|\normalfont{\textsf{AC}}| \asymp |\normalfont{\textsf{AR}}|$, feasible studentization requires $\bar{r} = \tilde{o}(|\normalfont{\textsf{AR}}|^{1/6})$, which is stronger than the sufficient condition for oracle normality.
This reflects the additional need to estimate both representers and the local variance profile accurately on the variance scale.
\medskip
\noindent \textbf{An entrywise confidence interval.} For any significance level $\gamma \in (0,1)$, Theorem~\ref{thm:var.est.asymp} justifies the entrywise $(1-\gamma)100\%$ confidence interval for $A_{ij}$:
\begin{align}
\texttt{CI}(\gamma) \coloneqq \left[ \hA_{ij} \pm z_{1-\gamma/2} \cdot \widehat{\upsilon} \right], \label{eq:conf.iv}
\end{align}
where $z_{1-\gamma/2}$ is the $(1-\gamma/2)$ quantile of the standard normal distribution.
\subsection{Estimating Heteroskedastic Noise}
Theorem~\ref{thm:var.est.asymp} reduces feasible entrywise inference to estimating the variance entries that appear in the leading variance formula. We now provide such an estimation procedure.
\subsubsection{Moment Construction} \label{sec:b.strategy}
Theorem~\ref{thm:consistency} establishes recovery of the first moment $A_{ij}$. For variance estimation, we leverage the elementary identity $\sigma^2_{ij} = B_{ij} - A_{ij}^2$, where $B_{ij} \coloneqq \mathbb{E}[\widetilde{Y}^2_{ij}\mid\boldsymbol{u}_i, \boldsymbol{v}_j]$ is the conditional second moment. Let $\bSigma = [\sigma^2_{ij}]$ and $\bB = [B_{ij}]$. If both $\bSigma$ and $\bA$ are low-rank, then $\bB$ inherits a low-rank structure as well, since
\begin{align}
\rank(\bB) \le \rank(\bSigma) + \rank(\bA)^2. \label{eq:rankB}
\end{align}
This observation suggests the following \texttt{SNN}-based procedure.
\begin{enumerate}
\item {\bf First moment estimation}: Apply \texttt{SNN}~to $\boldsymbol{Y}$ and clip the resulting estimate to $\hA^\text{clip}_{ij} \in [-1,1]$.
\item {\bf Second moment estimation}: Apply \texttt{SNN}~to the elementwise squared matrix $\boldsymbol{Y} \circ \boldsymbol{Y}$, obtaining the raw estimate $\widehat{B}_{ij}$ and its clipped counterpart $\widehat{B}^\text{clip}_{ij} \in [0, 1 +\sigma^2]$, where $\sigma$ is the noise bound from Assumption~\ref{assump:subg}. Let $t_b$ denote the spectral threshold used for this second-moment regression.
\item {\bf Variance estimation}: Define
$\widehat{\sigma}^2_{ij} \coloneqq \widehat{B}_{ij} - \big( \hA^{\text{clip}}_{ij} \big)^2$ and $\widehat{\sigma}^{2, \text{clip}}_{ij} \coloneqq \widehat{B}^\text{clip}_{ij} - \big( \hA^{\text{clip}}_{ij} \big)^2$.
\end{enumerate}
This moment-subtraction strategy is not specific to \texttt{SNN}. In principle, it may be paired with any matrix completion method capable of estimating first and second moments under the relevant observation structure.
\subsubsection{Second-Moment Conditions} \label{sec:var.est.assump}
We impose second-moment analogues of the assumptions used for the mean matrix.
Write $B_{ij} \coloneqq \langle \boldsymbol{\theta}_i, \boldsymbol{\varphi}_j \rangle$, where $\boldsymbol{\theta}_i, \boldsymbol{\varphi}_j \in \Rb^\zeta$ are latent factors associated with row $i$ and column $j$, and $\zeta$ is compatible with the rank bound in \eqref{eq:rankB}.
\begin{assumption} \label{assump:linear.span.B}
Conditioned on $\Ec$ and given anchor sets $(\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}})$, let $\boldsymbol{\theta}_i \in {\normalfont \text{span}}(\{ \boldsymbol{\theta}_\ell: \ell \in \normalfont{\textsf{AR}} \})$ and $\boldsymbol{\varphi}_j \in {\normalfont \text{span}}(\{ \boldsymbol{\varphi}_q: q \in \normalfont{\textsf{AC}} \})$.
\end{assumption}
\begin{assumption} \label{assump:spectra.B}
Given anchor sets $(\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}})$, the condition number $\kappa$ of $\mathbb{E}[\boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}} \circ \boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}} \mid \Ec]$ satisfies $\kappa^{-1} \ge c$, and $\| \mathbb{E}[\boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}} \circ \boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}} \mid \Ec] \|_F^2 \ge c' \cdot | \normalfont{\textsf{AR}}| \cdot |\normalfont{\textsf{AC}}|$ for constants $c, c' > 0$.
\end{assumption}
Assumption~\ref{assump:linear.span.B} is the second-moment counterpart of the row- and column-span conditions in Assumptions~\ref{assump:linear_span} and \ref{assump:subspace}. Assumption~\ref{assump:spectra.B}, analogous to Assumption~\ref{assump:spectra}, is the corresponding spectral condition: it requires the second-moment anchor block to have a sufficiently strong low-rank signal.
\subsubsection{Variance-Recovery Guarantees.} \label{sec:var.est.results}
We now state guarantees for heteroskedastic variance estimation in two regimes: bounded noise and general subgaussian noise.
\begin{theorem} \label{thm:var.est}
Let the setup of Theorem~\ref{thm:consistency} hold, and suppose Assumptions~\ref{assump:linear.span.B} and \ref{assump:spectra.B} also hold.
Let $t_b = \bar{\lambda} \coloneqq \rank(\mathbb{E}[\boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}} \circ \boldsymbol{Y}_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}}\mid\Ec])$, where $\bar{\lambda} \le \bar{\tau} + \bar{r}^2$ and $\bar{\tau} \coloneqq \rank(\bSigma_{\normalfont{\textsf{AR}}, \normalfont{\textsf{AC}}})$.
\begin{enumerate} [label=(\roman*)]
\item \underline{Bounded noise.} Suppose the noise entries $\varepsilon_{ij}$ are bounded. Then, conditional on $\Ec$, with probability at least $1 - \Oc(|\normalfont{\textsf{ARC}}|^{-10})$
\begin{align}
\left| \widehat{\sigma}^{2,\emph{clip}}_{ij} - \sigma^2_{ij} \right|
\lesssim
\frac{\left(1+ \sigma^2\right) \left(\bar{r}^2 + \bar{\tau} \right)}{ \min \left\{ \sqrt{|\normalfont{\textsf{AR}}|}, \sqrt{|\normalfont{\textsf{AC}}|} \right\} }
+ \frac{(1 + \sqrt{\sigma})(\bar{r} + \sqrt{\bar{\tau}}) L^{1/4}}{\min\left\{ |\normalfont{\textsf{AC}}|, |\normalfont{\textsf{AR}}|^2 / \left((1+\sigma^2) L \right) \right\}^{1/4}}.
\label{eq:var.est.bounded}
\end{align}
\item \underline{General subgaussian noise.} Suppose the noise entries $\varepsilon_{ij}$ satisfy Assumption~\ref{assump:subg} and $|\widehat{B}_{ij} - B_{ij}\mid\le M$ for some deterministic $M > 0$. Then
\begin{align}
\left| \mathbb{E}\left[ \widehat{\sigma}^2_{ij} - \sigma^2_{ij} \mid \Ec \right] \right|
&\lesssim \frac{\left(1 + \sigma^4 \right) \left(\bar{r}^2 + \bar{\tau} \right) }{\min \left\{ \sqrt{|\normalfont{\textsf{AR}}|} L^{-3/2}, \sqrt{|\normalfont{\textsf{AC}}|} \right\}}
+ \frac{ \left(1 + \sigma^6 \right) \left(\bar{r}^2 + \bar{\tau} \right)^{3/2}}{\min \left\{ |\normalfont{\textsf{AR}}| L^{-3}, |\normalfont{\textsf{AC}}| \right\}}
+ \frac{M}{|\normalfont{\textsf{ARC}}|^{10}}.
\label{eq:var.est.general}
\end{align}
\end{enumerate}
\end{theorem}
\noindent \textbf{Bounded noise.}
When $\bar{r}^4 + \bar{\tau}^2 = \tilde{o}( \min \{ |\normalfont{\textsf{AR}}|, |\normalfont{\textsf{AC}}| \})$, the clipped variance estimator is consistent.
This result parallels Theorem~\ref{thm:consistency}: under bounded noise, the squared deviation $Y_{ij}^2 - B_{ij}$ remains subgaussian, allowing the techniques from Theorem~\ref{thm:consistency} to extend.
In this regime, Theorem~\ref{thm:var.est}(i) can be used inside Theorem~\ref{thm:var.est.asymp} by taking $\Lambda_\sigma$ to be the right-hand side of \eqref{eq:var.est.bounded}. A union bound over the variance entries appearing in \eqref{eq:var.bound} then gives $p_\sigma \coloneqq \Oc(|\normalfont{\textsf{ARC}}|^{-9})$.
\medskip \noindent \textbf{General subgaussian noise.} In this broader regime, squared observations are sub-exponential rather than subgaussian. This changes the concentration behavior of the anchor block and prevents a direct high-probability analogue of the bounded-noise result.
The theorem therefore establishes asymptotic unbiasedness, not consistency. In particular, it does not by itself provide the uniform variance estimation bound $\Lambda_\sigma$ required for Theorem~\ref{thm:var.est.asymp}.
A sufficient condition for asymptotic unbiasedness is $\bar{r}^2 + \bar{\tau} = \tilde{o}( \min \{ \sqrt{|\normalfont{\textsf{AR}}|}, \sqrt{|\normalfont{\textsf{AC}}|} \} )$ and $M = o\left(|\normalfont{\textsf{ARC}}|^{10} \right)$. In contrast to the bounded noise setting, the balanced regime $|\normalfont{\textsf{AC}}| \asymp |\normalfont{\textsf{AR}}|$ is most favorable for the stated unbiasedness guarantee.
\subsection{Scope of the Heteroskedastic Extension}
Recent panel-factor inference results have primarily emphasized homoskedastic, or effectively homoskedastic, noise. In this respect, Theorems~\ref{thm:var.est.asymp} and \ref{thm:var.est} provide a useful heteroskedastic extension. They allow uncertainty to be calibrated at the same local, target-specific level as the point estimator itself, which is meaningful when some rows, columns, or anchor blocks are intrinsically noisier than others. This refinement should not be read as a universal replacement for existing inference methods: when homoskedasticity is credible, or when the observation design has exploitable global structure, specialized panel methods may remain simpler or sharper. Moreover, feasible heteroskedastic inference carries a real cost: it requires estimating the local variance profile in addition to the two representers, and the strongest guarantees currently available rely on bounded noise.
Theorem~\ref{thm:var.est} may also be useful beyond its role in Theorem~\ref{thm:var.est.asymp}.
The variance matrix describes irreducible outcome uncertainty, identifies entries whose outcomes are intrinsically unstable, supports risk-sensitive decisions, and can serve as a nuisance input for other inferential procedures. From this perspective, Theorem~\ref{thm:var.est} addresses a distinct statistical problem: recovering a heteroskedastic variance matrix from a single partially observed matrix under MNAR sampling. Thus, beyond enabling studentized \texttt{SNN}~intervals, the heteroskedastic extension may be serviceable whenever entry-specific uncertainty is itself a substantive object of interest.
\section{Simulation Studies} \label{sec:sims}
We conduct two simulation studies to evaluate the finite-sample performance of \texttt{SNN}. The first examines the coverage and length of the feasible entrywise confidence intervals developed in Section~\ref{sec:var.est} under heteroskedastic noise. The second compares \texttt{SNN}~with the panel-factor estimator of \cite{yan2024entrywiseinferencemissingpanel} in a staggered adoption setting specifically tailored to their method.
\subsection{Coverage Study}
We first evaluate the finite-sample performance of the feasible \texttt{SNN}~confidence interval in \eqref{eq:conf.iv} under heteroskedastic noise.
\subsubsection{Setup} For each matrix dimension $n \in \{25, 50, 100, 250, 500\}$, we generate an $n \times n$ outcome matrix. We designate the lower-right entry $(n,n)$ as the target, mask that entry, and observe all remaining entries. The signal matrix is deterministic for each $n$ and is constructed from smooth Fourier row and column features. Define $\rho_i \coloneqq (i-1/2)/n$, $\theta_j \coloneqq (j-1/2)/n+0.173$, and
$f(t) \coloneqq \big[1,\sqrt{2}\sin(2\pi t),\sqrt{2}\cos(2\pi t),\sqrt{2}\sin(4\pi t),\sqrt{2}\cos(4\pi t)\big]^\top$.
Let the diagonal matrix of Fourier coefficients be $\boldsymbol{\Gamma} \coloneqq \operatorname{diag}(0.50,0.38,0.30,0.23,0.17)$ and define the unscaled signal by
$A^{\mathrm{raw}}_{ij} \coloneqq f(\rho_i)^\top \boldsymbol{\Gamma} f(\theta_j)$.
The coefficients place the greatest weight on the intercept and first harmonic, with smaller weights on the second harmonic. The phase shift in $\theta_j$ distinguishes the row and column factors while retaining the same smooth Fourier basis. We then multiply the entire matrix by a common scaling constant so that $|A_{ij}| \le 0.8$ for every $(i,j)$. By construction, $\rank(\bA) = 5$, and the target row and target column belong to the spans of the corresponding observed rows and columns.
For each observed entry, the outcome is generated as $Y_{ij} = A_{ij} + \varepsilon_{ij}$. We consider two noise distributions. In the Gaussian setting, $\varepsilon_{ij} \sim \mathcal{N}(0, \sigma^2_{ij})$, whereas in the bounded noise setting, $\varepsilon_{ij} \sim \texttt{Uniform}\big[-\sqrt{3} \sigma_{ij}, \sqrt{3} \sigma_{ij} \big]$. Both distributions have variance $\sigma^2_{ij}$. The heteroskedastic variance matrix is generated as a bounded, positive, low-rank Fourier surface and normalized to satisfy $(1/n^2) \sum_{i=1}^n \sum_{j=1}^n \sigma^2_{ij} = \sigma^2$, where $\sigma = 0.3$. Thus, the two experiments share the same entrywise variance profile and differ only in the shape of the noise distribution.
For each matrix dimension and noise distribution, we conduct 100 independent replications. \texttt{SNN}~is supplied with the oracle signal rank of five so that the experiment isolates the performance of the inferential procedure rather than the accuracy of rank selection. In each replication, we construct the nominal 90\% confidence interval $\big[ \hA_{nn} \pm z_{0.95} \widehat{\upsilon} \big]$, where $\widehat{\upsilon}$ is obtained from the plug-in variance estimator in \eqref{eq:hvar.asymp}. Its entrywise noise-variance inputs are estimated using the moment construction in Section~\ref{sec:b.strategy}, which applies \texttt{SNN}~to both the observed outcome matrix and its entrywise square.
\subsubsection{Results} Figure~\ref{fig:inference} reports three measures of finite-sample performance. Panel (a) shows the RMSE of the estimated leading variance, Panel (b) shows empirical coverage, and Panel (c) shows average interval length. The variance estimation error declines sharply as $n$ grows. At the same time, empirical coverage remains close to the nominal 90\% level under both noise distributions, while the the average interval length steadily decreases. Taken together, these results indicate that increasing the local sample size improves precision without producing visible undercoverage.
The uniform noise experiment lies within the bounded-noise setting of Theorem~\ref{thm:var.est}(i), and its results are consistent with the feasible inference guarantee in Theorem~\ref{thm:var.est.asymp}. The Gaussian experiment exhibits similar behavior but lies outside the scope of the current feasible inference theory. The Gaussian results therefore provide empirical evidence that the procedure may perform well beyond the bounded noise regime, rather than a formal verification of coverage in that broader setting.
\begin{figure}[t!]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]
{images/sims/snn_ci_lead_variance_rmse.png}
\caption{Lead variance estimation.}
\label{fig:lead.var}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]
{images/sims/snn_ci_coverage.png}
\caption{Empirical coverage.}
\label{fig:coverage}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]
{images/sims/snn_ci_length.png}
\caption{Average interval length.}
\label{fig:length}
\end{subfigure}
\caption{Finite-sample performance of feasible \texttt{SNN}~inference under heteroskedastic noise: (a) RMSE of the estimated leading variance, (b) empirical coverage of nominal 90\% confidence intervals, and (c) average confidence-interval length.}
\label{fig:inference}
\end{figure}
\subsection{Comparison with Related Panel-Factor Work}
We next compare \texttt{SNN}~with the estimator of \cite{yan2024entrywiseinferencemissingpanel}, which we refer to as the \texttt{Yan-Wainwright}~estimator. We reproduce their ``Scaling with Rank'' experiment from Section 5.2 of their paper, which uses a staggered adoption observation pattern.
\subsubsection{Setup} The simulation uses a $500 \times 500$ matrix, targets the final entry, and varies the signal rank $r$ from 1 to 30. Following \cite{yan2024entrywiseinferencemissingpanel}, we consider their two data-generating cases. Both estimators are supplied with the oracle rank $r$, so the comparison isolates estimation performance rather than rank-selection accuracy. We refer the reader to \cite{yan2024entrywiseinferencemissingpanel} for the remaining details of the simulation design.
\subsubsection{Results} Figure~\ref{fig:yw_snn} reports the MSEs of each estimator as a function of the signal rank.
In Case 1, the errors of \texttt{SNN}~and the \texttt{Yan-Wainwright}~estimator are closely aligned across the full range of ranks. In Case 2, \texttt{SNN}~has modestly higher MSE. This difference is unsurprising: the simulated observation pattern has exactly the staggered adoption structure for which the \texttt{Yan-Wainwright}~estimator is designed. \texttt{SNN}, by contrast, is intended for a broader class of observation masks and bases each estimate only on the information contained in a local anchor cross.
The comparison illustrates the trade-off between specialization and robustness. When a credible global structure, such as staggered adoption, is available, an estimator tailored to that structure may achieve better finite-sample performance. \texttt{SNN}~instead remains applicable when the global observation pattern is irregular, provided that sufficiently informative local anchor crosses are available, as demonstrated in Section~\ref{sec:teaser}. \texttt{SNN}~should therefore be viewed as complementary to, rather than a replacement for, panel-specific estimators.
\begin{figure}[t!]
\centering
\begin{subfigure}[b]{0.35\textwidth}
\centering
\includegraphics[width=\linewidth]
{images/sims/yw_snn_rank_case1.png}
\caption{Case 1.}
\label{fig:stagger.1}
\end{subfigure}
\quad
\begin{subfigure}[b]{0.35\textwidth}
\centering
\includegraphics[width=\linewidth]
{images/sims/yw_snn_rank_case2.png}
\caption{Case 2.}
\label{fig:stagger.2}
\end{subfigure}
\caption{MSE as a function of signal rank for \texttt{SNN}~and the estimator of \cite{yan2024entrywiseinferencemissingpanel}: (a) Case 1 and (b) Case 2.}
\label{fig:yw_snn}
\end{figure}
\section{Conclusion} \label{sec:conclusion}
This article develops a causal framework for matrix completion under MNAR sampling.
The framework draws on ideas from synthetic controls and panel data, but targets a different regime from much of the existing literature. Unlike standard MNAR matrix completion approaches, it does not require positivity or independent observation indicators; an entry may have zero probability of being observed, and the observation pattern may be dependent across entries. Unlike classical panel-factor methods, it does not require a prescribed global geometry such as staggered adoption, tall-wide panels, or rigid block sparsity. Instead, identification is local: a target entry can be recovered whenever the observed mask contains a sufficiently informative anchor structure around that entry.
To implement this principle, we introduce \texttt{SNN}, an estimator that combines the locality of nearest-neighbor methods with the reweighting logic of synthetic controls. We establish finite-sample entrywise guarantees for mean recovery, derive asymptotic normality under heteroskedastic noise, and develop feasible studentized inference. We also extend the same local-anchor logic to heteroskedastic variance estimation, proving consistency under bounded noise and asymptotic unbiasedness under general subgaussian noise. The simulation evidence is consistent with these theoretical results and illustrates the robustness of \texttt{SNN}~across missingness mechanisms and observation patterns.
Looking ahead, \texttt{SNN}’s position between classical matrix completion and panel-factor methods suggests two related directions. First, a sharper comparison across MCAR sampling, local MNAR designs, and structured panel geometries would clarify when local anchor methods are preferable to procedures tailored to global random sampling or panel-specific observation patterns. Second, it would be valuable to develop hybrid procedures that exploit global structure when it is reliable while retaining the local anchor logic that makes \texttt{SNN}~robust to irregular MNAR patterns. For example, global factor estimates could stabilize anchor discovery and spectral denoising, whereas local anchor crosses could deliver target-specific recovery and inference in regions where a global observation model is weak or misspecified. Such hybrids may provide a more adaptive bridge between classical matrix completion, panel-factor methods, and the local MNAR framework developed here.
\bibliographystyle{alpha}
\bibliography{bib}
\clearpage