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.
136,135 characters
\begin{center}
{\bf{\LARGE{Entrywise Inference for Missing Panel Data: \\ A Simple
and Instance-Optimal Approach}}}
\vspace*{.2in}
{\large{
\begin{tabular}{ccc}
Yuling Yan$^{\star,\dagger}$ && Martin J. Wainwright$^{\star,\dagger,
\ddagger,+}$
\end{tabular}
}}
\vspace*{.2in}
\begin{tabular}{c}
Institute for Data, Systems, and Society$^\star$ \\
Laboratory for Information and Decision Systems$^\dagger$ \\
Department of Electrical Engineering and Computer
Sciences$^\ddagger$\\
Department of Mathematics$^+$ \\
Massachusetts Institute of Technology \\
\texttt{\{yulingy,mjwain\}@mit.edu}
\end{tabular}
\date{}
\vspace*{.2in}
\begin{abstract}
Longitudinal or panel data can be represented as a matrix with rows
indexed by units and columns indexed by time. We consider
inferential questions associated with the missing data version of
panel data induced by staggered adoption. We propose a
computationally efficient procedure for estimation, involving only
simple matrix algebra and singular value decomposition, and prove
non-asymptotic and high-probability bounds on its error in
estimating each missing entry. By controlling proximity to a
suitably scaled Gaussian variable, we develop and analyze a
data-driven procedure for constructing entrywise confidence
intervals with pre-specified coverage. Despite its simplicity, our
procedure turns out to be instance-optimal: we prove that the width
of our confidence intervals match a non-asymptotic instance-wise
lower bound derived via a Bayesian Cram\'{e}r--Rao argument. We
illustrate the sharpness of our theoretical characterization on a
variety of numerical examples. Our analysis is based on a general
inferential toolbox for SVD-based algorithm applied to the matrix
denoising model, which might be of independent interest.
\end{abstract}
\end{center}
\section{Introduction}
Longitudinal or panel data consists of a collection of observations of
units (e.g., individuals, companies, countries) that are collected
over time. In many applications, a subset of units are exposed to a
``treatment'' (e.g., drugs, regulations, or governmental
interventions) beginning at some time. Given data of this type, it is
frequently of interest to estimate counterfactual quantities, such as
what would have been their response if they had not been treated.
Such estimates underpin inferential methods for the treatment effect,
corresponding to the difference between the treated and untreated
responses. From the perspective of the untreated observations,
performing treatment can be modeled as inducing \emph{missing data}:
we no longer have observations of a given unit's untreated response at
all times after treatment. This problem and its variants, known as
(causal) \emph{inference with panel data}, find wide applications in
economics, social sciences, and biomedical research
(e.g.,~\citep{hu2008ownership,lewis2008economics,huang2008causal,imbens2015causal}).
In this paper, we assume that the treatment assignments follow the
staggered adoption
design~\citep{athey2022design,shaikh2021randomization}, meaning that
units may begin treatment at possibly different times, but that once
initiated, the treatment is irreversible. We can collect data for the
untreated unit/period pairs into a matrix, whose rows index units and
columns index periods. Performing treatment can be viewed as inducing
a form of missingness in this matrix; we refer the reader
to~\Cref{fig:intro} for an illustration of the induced missing
pattern.
Broadly speaking, in the literature on causal panel data, there are at
least two main approaches for imputing missing entries in panel data:
those based on unconfoundedness
(e.g.,~\citep{rosenbaum1983central,imbens2015causal}), and those based
on synthetic controls
(e.g.,~\citep{abadie2003economic,abadie2015comparative,abadie2021using}).
In order to understand how these methods work, suppose that we are
interested in estimating the missing outcome for unit $i$ at period
$t$. Approaches based on unconfoundedness seek to identify a subset of
untreated units whose outcomes (before unit $i$ was treated) are
similar to those of unit $i$, and use their observed outcomes at
period $t$ to estimate unit $i$'s missing outcome. On the other hand,
approaches based on synthetic controls estimate unit $i$'s missing
outcome at time $t$ by using a weighted average of all untreated
units' observed outcomes at that time, where the weights are typically
determined by solving a re problem using a set of predictors
that are unaffected by the treatment.
\newlength{\mywidth}
\setlength{\mywidth}{0.55cm}
\begin{figure}[t]
\centering \begin{tikzpicture}
\matrix (m) [matrix of math nodes,
nodes={draw, minimum height=\mywidth, minimum width=\mywidth, anchor=center},
column sep=-\pgflinewidth, row sep=-\pgflinewidth, font=\small]{
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \\
|[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & |[fill=gray!20]| \mathrm{C} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \mathrm{T} & \\
};
\node[above=-1mm, font=\small] at (m.north) {time periods};
\node[left=0mm, font=\small] at (m.west) {units};
\draw[decorate,decoration={brace,amplitude=6pt,raise=1pt}, font=\small]
(m-1-15.north east) -- (m-4-15.south east) node [black,midway,xshift=30pt,yshift=-8pt,anchor=south] {group 1};
\draw[decorate,decoration={brace,amplitude=6pt,raise=1pt}, font=\small]
(m-5-15.north east) -- (m-6-15.south east) node [black,midway,xshift=30pt,yshift=-8pt,anchor=south] {group 2};
\draw[decorate,decoration={brace,amplitude=6pt,raise=1pt}, font=\small]
(m-7-15.north east) -- (m-9-15.south east) node [black,midway,xshift=30pt,yshift=-8pt,anchor=south] {group 3};
\draw[decorate,decoration={brace,amplitude=6pt,raise=1pt}, font=\small]
(m-10-15.north east) -- (m-11-15.south east) node [black,midway,xshift=30pt,yshift=-8pt,anchor=south] {group 4};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}, font=\small]
(m-11-1.south west) -- (m-11-4.south east) node [black, midway, yshift=-15pt] {stage 1};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}, font=\small]
(m-11-5.south west) -- (m-11-8.south east) node [black, midway, yshift=-15pt] {stage 2};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}, font=\small]
(m-11-9.south west) -- (m-11-11.south east) node [black, midway, yshift=-15pt] {stage 3};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}, font=\small]
(m-11-12.south west) -- (m-11-15.south east) node [black, midway, yshift=-15pt] {stage 4};
\end{tikzpicture}
\caption{An illustration of a panel data with missingness induced by
staggered adoption design. The labels C and T refer to ``control''
(i.e., untreated) and ``treated'' respectively. From the
perspective of the control group, all entries marked with T
correspond to missing data, and a key problem is to impute these
entries.} \label{fig:intro}
\end{figure}
Athey et al.~\citep{athey2021matrix} observed that both the
unconfoundedness and synthetic control methods can be related to the
literature on low-rank matrix completion
(e.g.,~\citep{ExactMC09,Negahban2012restricted}). Based on this
connection, they proposed to impute the missing potential outcomes via
the standard convex relaxation for low-rank matrix
completion---namely, minimizing a least-squares objective with a
nuclear norm regularization. When the potential outcome matrix is
approximately low-rank, this approach can be equipped with attractive
guarantees. Following this general avenue, the past few years have
witnessed substantial progress in the development and analysis of
matrix completion algorithms tailored for panel data with missingness
(e.g.,~\citep{bai2021matrix,cahan2023factor,agarwal2023causal,choi2023matrix}).
Matrix completion methods based on nuclear norm relaxation are known
to be optimal when entries are missing completely at random
(e.g.,~\cite{wainwright2019high}), but there are currently no such
optimality guarantees for missingness induced by staggered adoption.
\subsection{Our contributions}
With this context, the main contribution of this paper is to propose a
simple procedure, substantially less computationally intensive than
nuclear norm relaxation, for inferring individual treatment effects in
panel data. Most importantly, we are able to show that our
procedure---despite its simplicity---is unimprovable in a sharp
instance-optimal sense. More precisely:
\begin{itemize}
\item \textbf{Entrywise guarantees:} We give non-asymptotic bounds and
distributional characterizations of the error in estimating each
matrix entry. Using the distributional characterization, we provide
a data-driven procedure for computing confidence intervals for the
unobserved entries with prescribed coverage.
\item \textbf{Sharp instance-wise optimality:} We demonstrate that the
length of our confidence intervals match the Bayesian Cram\'er-Rao
lower bound associated with a genie-aided version of the problem (in
which an oracle provides partial information about the unknown
truth).
\item \textbf{Computational simplicity:} Our algorithm only involves
basic matrix operations and singular value decompositions, making it
computationally faster than standard convex relaxations (which
typically involve solving semi-definite programs).
\end{itemize}
Our procedure and theory is built upon an inferential toolbox
developed for the SVD algorithm in the context of low-rank matrix
denoising, which might be of independent
interest. See~\Cref{appendix:proof-thm-denoising} for details.
\paragraph{A preview:}
\Cref{FigPreview} provides a high-level preview of some of the
consequences of our results, and the sharpness of our theoretical
predictions. In particular, let $\Mstar \in \real^{N \times T}$ be
the underlying matrix of outcomes for the untreated group.
\begin{figure}[h!]
\begin{center}
\begin{tabular}{ccc}
\widgraph{0.45\textwidth}{figs/fig_qqplot_entry1} &&
\widgraph{0.45\textwidth}{figs/fig_incoherence_case1} \\
(a) && (b)
\end{tabular}
\caption{(a) Illustration of the asymptotic normality of the
rescaled error $\ensuremath{R}_{i,t}$ from
equation~\eqref{EqnDefnZscore} for entry $(i,t) = (500, 500)$
from a $500 \times 500$ matrix. Shown is a Q-Q plot of the
empirical quantile of the suitably rescaled estimation error
versus the standard normal quantiles; consistent with our
theory, excellent agreement is shown. (b) Illustration of the
sharpness of our theoretical predictions: we construct a family
of matrix missing data problems, whose difficultly is calibrated
by an ``incoherence parameter''. Plots of actual mean-squared
error (MSE) of our procedure (red circles) to a lower bound
derived using a Bayesian Cram\'{e}r-Rao lower bound (blue solid
line). See~\Cref{sec:numerical} for the details of these
simulation studies.}
\label{FigPreview}
\end{center}
\end{figure}
Our methodological innovation is to propose an
estimate $\Mhat$ that is attractively simple: its computation involves
only SVDs and other basic matrix operations. For each unobserved
entry $(i, t)$, we define the rescaled error term
\begin{align}
\label{EqnDefnZscore}
\ensuremath{R}_{i,t} = ( \widehat{M}_{i,t} - M_{i,t}^\star ) /
\sqrt{\ensuremath{\widehat{\ensuremath{\gamma}}}_{i,t}},
\end{align}
where $\ensuremath{\widehat{\ensuremath{\gamma}}}_{i,t}$ is a data-dependent scaling factor. On the
achievable side, we provide a non-asymptotic decomposition of the
rescaled error $\ensuremath{R}_{i,t}$, a particular consequence of which is
that it converges in distribution to a standard $N(0,1)$ variate;
see~\Cref{thm:distribution,thm:distribution-general} for statements of
these results. Panel (a) in~\Cref{FigPreview} provides an empirical
demonstration of this predicted entrywise asymptotic normality.
Moreover, we prove that these estimation-theoretic and inferential
guarantees are \emph{optimal in a strong sense:} our variance
estimates $\ensuremath{\widehat{\ensuremath{\gamma}}}_{i,t}$ converge to a population quantity
$\ensuremath{\ensuremath{\gamma}^*}_{i,t}$ that is defined via a non-asymptotic Bayesian
Cram\'{e}r--Rao lower bound; see~\Cref{thm:CRLB} and the discussion
following~\Cref{thm:distribution-general} for these claims. Panel (b)
provides empirical confirmation of the instance optimality of our
method; we construct an ensemble of problems whose difficulty is
indexed by an ``incoherence'' parameter. Panel (b) compares Monte
Carlo estimates of the mean-squared error obtained by our procedure
with the lower bound predicted by~\Cref{thm:CRLB}; note the excellent
agreement. We refer the reader to~\Cref{sec:numerical} for full
details on the simulation set-ups that underlie these empirical
studies.
\subsection{Related work}
So as to put our contributions in context, let us discuss some related
lines of work. Low-rank matrix completion has been studied in great
detail when entries are assumed to be missing completely at random;
the behavior of the standard convex relation based on the nuclear norm
minimization is now very well-understood. Initial
investigations~\citep{ExactMC09,recht2010guaranteed,CanTao10}
primarily focused on the noiseless scenario and examined the minimal
sample size required for exact recovery. In the context of noisy
matrix completion (where the observed entries are corrupted by random
noise), the last decade has seen the establishment of optimal
estimation guarantees for convex
relaxation~\citep{Negahban2012restricted,klopp2014noisy,chen2019noisy}. A
more recent line of work~\cite{chen2019inference,xia2021statistical}
has focused on debiasing techniques for convex relaxation, showing how
it is possible to construct confidence intervals for each entry of the
unknown matrix. We refer to the reader to the
survey~\cite{chi2018nonconvex} for an overview of other non-convex
approaches to matrix completion
(e.g.,~\citep{burer2003nonlinear,KesMonSew2010,srebro2004learning}).
Athey et al.~\cite{athey2021matrix} pioneered the use of convex matrix
completion for imputing potential outcomes in panel data. They
provided Frobenius norm bounds on the error in estimating the full
matrix; such a bound can be viewed as providing estimation guarantees
for averaged treatment effects. More recently, Choi and
Yuan~\cite{choi2023matrix} proposed a more refined approach that first
divides the missing entries into smaller groups, and then applies the
convex relaxation to each group. They derived entrywise error bounds
that are sharper than the Frobenius norm error bound from the
paper~\cite{athey2021matrix}, and also established asymptotic
normality for some statistics of interest. Agarwal et
al.~\cite{agarwal2023causal} proposed an algorithm based on synthetic
nearest-neighbors, for which they established error bounds and
asymptotic normality. There is also a line of
work~\cite{bai2021matrix,cahan2023factor} on models satisfying certain
factor constraints, to which we compare in more detail
in~\Cref{sec:comparison}. Other work on related but distinct
estimation problems for causal panel data include the
papers~\cite{farias2021learning,lei2023estimating,choi2023inference,agarwal2020synthetic}.
Our algorithm involves spectral techniques, and it is natural that our
analysis of it has connections to past work on spectral
methods~\cite{abbe2017entrywise,cai2019subspace,yan2021inference,zhou2023deflated}.
Notably, some past
work~\citep{el2015impact,sur2017likelihood,ma2017implicit,chen2020bridging}
has made effective use of ``leave-one-out'' methods; see the
monograph~\citep{chen2020spectral} for an overview. Among the
technical contributions of this paper is a natural generalization of
this idea, which we refer to as ``leave-one-block-out''. As will be
clarified by our analysis, this extension is essential in providing a
sharp characterization of the subspace perturbation error that arises
in the panel data model.
\paragraph{Notation:} For a positive integer $n$, we adopt the
shorthand $[n] \coloneqq \{ 1 , \ldots , n \}$, For a matrix $\bm{A}
\in \mathbb{R}^{n_1 \times n_2}$ and subsets $\mathcal{I} = \{ i_1,
i_2, \ldots, i_I\} \subseteq [n_1]$ and $\mathcal{J} = \{ j_1, j_2,
\ldots, j_J \} \subseteq [n_2]$, define $\bm{A}_{\mathcal{I} ,
\mathcal{J}} \in \mathbb{R} ^ { I \times J}$ to be the submatrix of
$\bm{A}$ comprising the intersection of its $I$ rows indexed by
$\mathcal{I}$ and $J$ columns indexed by $\mathcal{J}$, specifically $
(\bm{A}_{\mathcal{I}, \mathcal{J}})_{k, l}=A_{i_k, j_l} $ for any $k
\in [ I ]$ and $l \in [ J ]$. In addition, we abbreviate
$\bm{A}_{\mathcal{I}, \cdot} \coloneqq \bm{A}_{ \mathcal{I} , [ n_2 ]
}$ and $\bm{A}_{ \cdot, \mathcal{J} } \coloneqq \bm{A}_{ [n_1] ,
\mathcal{J}}$. If a subset $\mathcal{I} = \{ i \}$ is a singleton,
we may directly use the index $i$ to replace $\mathcal{I}$ for
simplicity. We use $\Vert \bm{A} \Vert$, $\Vert \bm{A}
\Vert_{\mathrm{F}}$ and $\Vert \bm{A} \Vert_{\infty}$ (respectively)
to denote the spectral norm, Frobenius norm and entrywise
$\ell_{\infty}$ norms of a matrix $\bm{A}$. For any symmetric matrices
$\bm{A}$ and $\bm{B} \in \mathbb{R}^{n\times n}$, the relation $\bm{A}
\succeq\bm{B}$ (resp.~$\bm{A}\preceq \bm{B}$) means that
$\bm{A}-\bm{B}$ (resp.~$\bm{B}-\bm{A}$) is positive semidefinite. For
any invertible matrix $\bm{H} \in \mathbb{R}^{r\times r}$ with SVD
$\bm{U} \bm{\Sigma} \bm{V}^\top$, define its sign matrix $\mathsf{sgn}
( \bm{H} ) \coloneqq \bm{U} \bm{V}^\top$. We also define the function
$\mathsf{svds}(\bm{A},r)$ to output the truncated rank-$r$ SVD
$(\bm{U}, \bm{\Sigma}, \bm{V})$ of $\bm{A}$, where the columns of
$\bm{U} \in \mathbb{R}^{n_1 \times r}$ and $\bm{V} \in \mathbb{R}^{n_2
\times r}$ are the top-$r$ left/right singular vectors, and
$\bm{\Sigma} \in \mathbb{R}^{r \times r}$ is a diagonal matrix
containing the top-$r$ singular values.
\section{Problem set-up}
In this section, we provide a more precise description of the
observation model, and its reformulation in terms of matrix estimation
with missing entries.
\subsection{Basic observation model}
\label{sec:model}
We consider a panel data setting in which there are $\Nunit$ units
over $\Time$ periods. For each unit $i \in [N]$ and at each time
period $t \in [T]$, we observe a response $Y_{i,t}$. The
interpretation of this observation depends on whether or not unit $i$
has undergone treatment by time $t$. We assume that treatment follows
the \emph{staggered adoption design}, meaning that each unit $i \in
[\Nunit]$ may differ in the time that it is first exposed to
treatment, and the treatment is irreversible
(e.g.,~\cite{athey2022design,athey2021matrix,shaikh2021randomization}).
More precisely, for each unit $i$, we let the integer $t_i$ denote the
time at which unit $i$ was first exposed to the treatment; for
completeness, we set $t_i = \infty$ if unit $i$ never undergoes the
treatment. Using the notation of potential outcomes, let $Y_{i,t}(0)$
be the mean outcome of unit $i$ at time $t$ \emph{if} the unit is
never exposed to the treatment. For a given unit $i$ and for all
times $t = 1, 2, \ldots, t_i-1$, we assume that the observed response
$Y_{i,t}$ is a noisy version of this mean outcome---that is
\begin{align}
Y_{i,t} & = Y_{i,t} (0) + E_{i,t} \qquad \mbox{for each $t = 1,
\ldots, t_i-1$,}
\end{align}
where $\{ E_{i,t} \}$ are independent noise variables. On the other
hand, for any time index $t \geq t_i$, we do not model any joint
structure between the observation $Y_{i,t}$ and $Y_{i,t}(0)$, because
we expect that the former might also depend on the adoption time
$t_i$.
For each unit $i$ and each time $t \geq t_i$, our goal is to estimate
the individual treatment effect (ITE)
\begin{align}
\tau_{i,t} \coloneqq Y_{i,t} - Y_{i,t}(0),
\end{align}
corresponding to the difference between the treated response $Y_{i,t}$
(that we observe) and and the (unobserved) untreated mean
$Y_{i,t}(0)$. Without further assumptions, none of our observations
give any information about $Y_{i,t}(0)$ for all $t \geq t_i$; thus,
the ITEs are not identifiable in general. Additional assumptions are
required to ensure identifiability, and following a line of previous
work, out approach is to assume that the untreated mean outcomes have
low-rank structure.
In order to state this assumption more precisely, we introduce an
$\Nunit \times \Time$ matrix that is populated by the untreated mean
outcomes---that is, matrix $\bm{M}^\star \in \RR^{N \times T}$ with
entries $M_{i,t}^\star = Y_{i,t}(0)$. Given the previously described
observation model, we observe
\begin{align}
Y_{i,t} & = M_{i,t}^\star + E_{i,t} \qquad \mbox{for each $i \in [N]$
and $t = 1, 2, \ldots, t_i-1$.}
\end{align}
Consequently, the problem of estimating each individual treatment
effect $\tau_{i,t}$ can reduced to estimating the missing entries
$\{M_{i,t}^\star \mid i\in[N], t \geq t_i\}$. As in a line of past
work~\cite{athey2021matrix,choi2023matrix}, we adopt the assumption
that $\bm{M}^\star$ is a low-rank matrix (i.e., has rank much smaller
than $\min \{N, T \}$), under which it becomes possible to estimate
$\bm{M}^\star$ even in the missing data setting given here.
\subsection{Formulation as low-rank matrix completion}
We now give a precise formulation in terms of low-rank matrix
completion.
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\begin{tikzpicture}
\matrix (m) [matrix of math nodes,
nodes={draw, minimum height=0.8cm, minimum width=1.2cm, anchor=center, inner sep=0.2pt},
column sep=-\pgflinewidth, row sep=-\pgflinewidth]{
\bm{1} & \bm{1} & \cdots & \bm{1} & \bm{1} & \bm{1} \\
\bm{1} & \bm{1} & \cdots & \bm{1} & \bm{1} & \bm{0} \\
\bm{1} & \bm{1} & \cdots & \bm{1} & \bm{0} & \bm{0} \\
\operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\ddots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \\
\bm{1} & \bm{0} & \cdots & \bm{0} & \bm{0} & \bm{0} \\
};
\node[above, font=\footnotesize] at (m-1-1.north) {$T_1$};
\node[above, font=\footnotesize] at (m-1-2.north) {$T_2$};
\node[above, font=\footnotesize] at (m-1-4.north) {$T_{k-2}$};
\node[above, font=\footnotesize] at (m-1-5.north) {$T_{k-1}$};
\node[above, font=\footnotesize] at (m-1-6.north) {$T_{k}$};
\node[left, font=\footnotesize] at (m-1-1.west) {$N_{1}$};
\node[left, font=\footnotesize] at (m-2-1.west) {$N_{2}$};
\node[left, font=\footnotesize] at (m-3-1.west) {$N_{3}$};
\node[left, font=\footnotesize] at (m-5-1.west) {$N_{k}$};
\node[below=1mm] at (m.south) {$\bm{\Omega}$};
\end{tikzpicture} &
\begin{tikzpicture}
\matrix (m) [matrix of math nodes,
nodes={draw, minimum height=0.8cm, minimum width=1.2cm, anchor=center, font=\footnotesize, inner sep=0.2pt},
column sep=-\pgflinewidth, row sep=-\pgflinewidth]{
|[fill=gray!20]| \bm{M}_{1,1}^\star & |[fill=gray!20]| \bm{M}_{1,2}^\star & |[fill=gray!20]| \cdots & |[fill=gray!20]| \bm{M}_{1,k-2}^\star & |[fill=gray!20]| \bm{M}_{1,k-1}^\star & |[fill=gray!20]| \bm{M}_{1,k}^\star \\
|[fill=gray!20]| \bm{M}_{2,1}^\star & |[fill=gray!20]| \bm{M}_{2,2}^\star & |[fill=gray!20]| \cdots & |[fill=gray!20]| \bm{M}_{2,k-2}^\star & |[fill=gray!20]| \bm{M}_{2,k-1}^\star & \bm{M}_{2,k}^\star \\
|[fill=gray!20]| \bm{M}_{3,1}^\star & |[fill=gray!20]| \bm{M}_{3,2}^\star & |[fill=gray!20]| \cdots & |[fill=gray!20]| \bm{M}_{3,k-2}^\star & \bm{M}_{3,k-1}^\star & \bm{M}_{3,k}^\star \\
\operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\ddots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \operatorname{\vphantom{\int^{0^{0^0}}}\smash[t]{\vdots}} & \\
|[fill=gray!20]| \bm{M}_{k,1}^\star & \bm{M}_{k,2}^\star & \cdots & \bm{M}_{k,k-2}^\star & \bm{M}_{k,k-1}^\star & \bm{M}_{k,k}^\star \\
};
\node[above, font=\footnotesize] at (m-1-1.north) {$T_1$};
\node[above, font=\footnotesize] at (m-1-2.north) {$T_2$};
\node[above, font=\footnotesize] at (m-1-4.north) {$T_{k-2}$};
\node[above, font=\footnotesize] at (m-1-5.north) {$T_{k-1}$};
\node[above, font=\footnotesize] at (m-1-6.north) {$T_{k}$};
\node[left, font=\footnotesize] at (m-1-1.west) {$N_{1}$};
\node[left, font=\footnotesize] at (m-2-1.west) {$N_{2}$};
\node[left, font=\footnotesize] at (m-3-1.west) {$N_{3}$};
\node[left, font=\footnotesize] at (m-5-1.west) {$N_{k}$};
\node[below=1mm] at (m.south) {$\bm{M}^\star$};
\end{tikzpicture}
\end{tabular}
\caption{An illustration of the the indicator matrix $\bm{\Omega}$
and the potential outcome matrix $\bm{M}^\star$ under staggered
adoption design. Here $\bm{1}$ (resp.~$\bm{0}$) denotes a matrix
whose entries are all one (resp.~zero). The gray blocks in the
right panel are the potential outcomes associated with the
untreated unit/period pairs.}
\label{fig:setup}
\end{center}
\end{figure}
We assume that the matrix $\MatTre{}$ has rank $r \ll \min
\{\Nunit,\Time \}$. Define the subset $\Omega \subseteq [N] \times
[T]$ of indices where $(i, t) \in \Omega$ if and only if unit
$i$ has not yet been treated at time $t$. In terms of this
notation, we can write our observation model in the compact
form
\begin{align}
\label{EqnObservationModel}
\bm{M} \mydefn \mathcal{P}_{ \Omega } \left ( \MatTre{} + \bm{E}
\right ),
\end{align}
where $\mathcal{P}_{\Omega}(\bm{X})$ denotes the Euclidean projection
of a matrix $\bm{X}$ onto the subspace of matrices supported on
$\Omega$, and $\bm{E}$ is a random noise matrix with
i.i.d.~$\mathcal{N} \big(0, \ensuremath{\omega}^{2}/(NT) \big)$ entries. Our
scaling of the noise variance ensures that $\bm{E}$ has Frobenius norm
close to $\ensuremath{\omega}$ with high probability, and facilitates
interpretation of results in the sequel.
Under the staggered adoption design, we can sort the $N$ units
according to the time $t_{i}$ at which they were first exposed to the
treatment. With this sorting, we have $t_{1} \geq t_{2} \geq \cdots
\geq t_{N}$, and we are guaranteed the existence of an integer $k \geq
1$ and partitions $N_{1}, \ldots N_{k}$ and $T_{1}, \ldots, T_{k}$
with $\sum_{i = 1}^{k}N_{i} = N$ and $\sum_{j = 1}^{k}T_{j} = T$ such
that the Boolean matrix
\begin{align}
\label{EqnDefnOmega}
\bm{\Omega} \in \{0,1\}^{N \times T} \quad \mbox{with entries
$\bm{\Omega}_{i,t} \mydefn \operatorname{\mathds{1}}_{(i, t) \in \Omega}$,}
\end{align}
associated with the subset $\Omega$ has the block structure as in the
left panel of~\Cref{fig:setup}. We also partition $\MatTre{}$
according to the pattern in $\bm{\Omega}$ to get the submatrices
$\bm{M}_{i, j}^{\star} \in \mathbb{R}^{N_{i} \times T_{j}}$, as shown
in the right panel of~\Cref{fig:setup}. For each $i\in[k]$ and $1 \leq
j \leq k-i$ define $\bm{M}_{i, j} = \bm{M}_{i, j}^{\star} + \bm{E}_{i,
j}$, where $\bm{E}_{i,j}$ is the corresponding submatrix of
$\bm{E}$.
Our goal is to design an algorithm that provides point estimates,
along with confidence intervals, for each unseen entry $M_{i,t}^\star$
as well as corresponding ITE $\tau_{i,t} = Y_{i,t} - M_{i,t}^\star$,
where $(i,t) \notin \Omega$. Since $Y_{i,t}$ is known, the estimation
and inference for these two sets of quantities are equivalent, hence
we will only present algorithms and results for $M_{i,t}^\star$ in the
following sections.
\subsection{A simple four-block structure}
Our estimator for the general case can be understood most easily by
understanding how it applies to a special four-block structure. More
precisely, suppose that a subset of $N_{1}$ units \emph{never} receive
the treatment, while the other $N_{2} = N - N_{1}$ units are exposed
to the irreversible treatment at time $T_{1} + 1$. Then our general
observation model can be written in the four-block form
\begin{align}
\MatTre{} = \left [ \begin{array}{cc} \MatTre{a} & \MatTre{b}
\\ \MatTre{c} & \MatTre{d}
\end{array} \right ], \qquad \bm{E} = \left [ \begin{array}{cc}
\bm{E}_{a} & \bm{E}_{b} \\
\label{eq:four-block-structure}
\bm{E}_{c} & \text{N/A}
\end{array} \right ].
\end{align}
Here the subscript $a$ denotes matrices of size $N_{1}$ by $T_{1}$;
the subscript $b$ denotes matrices of size $N_1$ by $T - T_1$; and so
on for the subscripts $c$ and $d$. As shown in the sequel, it
suffices to understand how to perform estimation and inference in this
scenario because our algorithm for the general staggered design first
reduces the problem to a simpler problem with four-block structure. In
view of the discussion at the end of the previous section, our goal is
to estimate and construct confidence intervals for the entries of
$\bm{M}_{d}^{\star}$, since they correspond to the counterfactual
outcomes associated with the control group.
\section{Estimation algorithms}
\label{sec:alg}
In this section, we first design an algorithm that is applicable to
the simpler four-block model~\eqref{eq:four-block-structure}. This
algorithm is relatively easy to describe, and is the basic building
block for our algorithm that applies to the general staggered design.
\subsection{Algorithm for the four-block design}
\begin{algorithm}[t]
\DontPrintSemicolon \SetNoFillComment
\textbf{Input:} Data matrix $ \bm{M}$, rank $r$ \\
\tcc{Step 1: Subspace Estimation} Compute the truncated rank-$r$ SVD
$(\bm{U}_{\mathsf{left}}, \bm{\Sigma}_{\mathsf{left}},
\bm{V}_{\mathsf{left}})$ of $\bm{M}_{\mathsf{left}}$. \\
Partition $\bm{U}_{\mathsf{left}}$ into two submatrices $\bm{U}_{1}$
and $\bm{U}_{2}$, where $\bm{U}_{1} \in \mathbb{R}^{N_{1}\times r}$
consists of its top $N_{1}$ rows and $\bm{U}_{2} \in
\mathbb{R}^{N_{2}\times r}$ consists of its bottom $N_{2}$
rows. \\
\tcc{Step 2: Matrix Denoising} Compute the truncated rank-$r$ SVD $(
\bm{U}_{\mathsf{upper}}, \bm{\Sigma}_{\mathsf{upper}},
\bm{V}_{\mathsf{upper}}) $ of $ \bm{M}_{\mathsf{upper}}$ \\
Partition $\bm{V}_{\mathsf{upper}}$ into two submatrices $\bm{V}_{1}$
and $\bm{V}_{2}$, where $\bm{V}_{1} \in \mathbb{R}^{T_{1}\times r}$
consists of its top $T_{1}$ rows and $\bm{V}_{2} \in
\mathbb{R}^{T_{2}\times r}$ consists of its bottom $T_{2}$ rows. \\
Compute the estimate $\widehat{ \bm{M}}_{b} \mydefn
\bm{U}_{\mathsf{upper}} \bm{\Sigma}_{\mathsf{upper}}
\bm{V}_{2}^{\top}$ of $\bm{M}_{b}^{\star}$. \\
\tcc{Step 3: Linear Regression} Compute the matrix $\widehat{
\bm{M}}_{d} \mydefn \bm{U}_{2}( \bm{U}_{1}^{\top} \bm{U}_{1})^{-1}
\bm{U}_{1}^{\top}\widehat{ \bm{M}}_{b}$. \\
\textbf{Output:} $\widehat{ \bm{M}}_{d}$ as estimate of
$\bm{M}_{d}^{\star}$
\caption{Estimating counterfactual outcomes: four-block
design \label{alg:4-blocks-Md}}
\end{algorithm}
We first define some necessary notation for the four-block model. We
begin by writing our observation matrix $\bm{M}$ in the form
\begin{align*}
\bm{M} = \left[\begin{array}{cc} \bm{M}_{a} & \bm{M}_{b}\\ \bm{M}_{c}
& \bm{?}
\end{array} \right] = \left[\begin{array}{cc}
\bm{M}_{a}^{\star} + \bm{E}_{a} & \bm{M}_{b}^{\star} +
\bm{E}_{b}\\ \bm{M}_{c}^{\star} + \bm{E}_{c} & \bm{?}
\end{array} \right].
\end{align*}
The unknown matrix $\bm{M}^{\star}$ has a (compact) singular value
decomposition, which we write in the form \mbox{$\bm{M}^\star =
\bm{U}^{\star} \bm{\Sigma}^{\star} \bm{V}^{\star\top}$,} where
$\bm{U}^{\star} \in \mathbb{R}^{N\times r}$ and $\bm{V}^{\star} \in
\mathbb{R}^{T\times r}$ have orthonormal columns, and
\mbox{$\bm{\Sigma}^{\star} = \mathsf{diag}(\sigma_{1}^{\star}, \ldots,
\sigma_{r}^{\star})$} is a diagonal matrix consisting of the ordered
singular values \mbox{$\sigma_{1}^{\star} \geq \sigma_{2}^{\star} \geq
\cdots \geq \sigma_{r}^{\star}$.} We partition the matrices of
singular vectors $\bm{U}^{\star} \in \real^{N \times r}$ and
$\bm{V}^{\star} \in \real^{T \times r}$ into two blocks as
\begin{align}
\bm{U}^{\star} = \left[\begin{array}{c} \bm{U}_{1}^{\star} \in
\mathbb{R}^{N_1 \times r} \\ \bm{U}_{2}^{\star} \in \real^{N_2
\times r}
\end{array} \right] \qquad
\text{and} \qquad \bm{V}^{\star} = \left[\begin{array}{c}
\bm{V}_{1}^{\star} \in \real^{T_1 \times r} \\ \bm{V}_{2}^{\star}
\in \real^{T_2 \times r}
\end{array} \right]. \label{eq:UV1-star}
\end{align}
It is also convenient to define
\begin{align*}
\bm{M}_{\mathsf{left}} = \left[\begin{array}{c} \bm{M}_{a}
\\ \bm{M}_{c}
\end{array} \right], \qquad
\bm{M}_{\mathsf{left}}^{\star} = \left[\begin{array}{c}
\bm{M}_{a}^{\star} \\
\bm{M}_{c}^{\star}
\end{array} \right],
\qquad
\bm{M}_{\mathsf{upper}} = \left[\, \bm{M}_{a}\;\; \bm{M}_{b}\,
\right], \qquad \bm{M}_{\mathsf{upper}}^{\star} = \left[\,
\bm{M}_{a}^{\star}\;\; \bm{M}_{b}^{\star}\, \right].
\end{align*}
To motivate our algorithm, consider the ideal setting where there is
no noise (i.e.~$ \bm{E} = \bm{0}$). Given that $\bm{M}_{a}^{\star}$ is
rank-$r$, we can first compute the matrix $\bm{U}^{\star}$ (up to a
rotation) by using the SVD of the left submatrix
$\bm{M}_{\mathsf{left}}^{\star}$, and then recover
$\bm{M}_{d}^{\star}$ by
\begin{align*}
\bm{M}_{d}^{\star} = \bm{U}_{2}^{\star}( \bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top} \bm{M}_{b}^{\star}.
\end{align*}
The idea behind the above formula is that, each column of $
\bm{M}^{\star}$ can be written uniquely as a linear combination of the
columns of $ \bm{U}^{\star}$. When $ \bm{U}_{1}^{\star}$ has full
rank, we can find the coefficients by regressing each column of $
\bm{M}_{b}^{\star}$ over the columns of $ \bm{U}_{1}^{\star}$. These
coefficients in turn allow reconstruction of $\bm{M}_{d}^{\star}$
using the columns of $ \bm{U}_{2}^{\star}$.
In the presence of random noise, the na\"{i}ve procedure needs to be
adjusted. Doing so leads to the following three-step procedure:
\begin{itemize}
\item (\textbf{Subspace estimation}) The goal of this step is to
estimate $\bm{U}^{\star}$. We compute the truncated rank-$r$ SVD of
$\bm{M}_{\mathsf{left}}$ as
\begin{align*}
\left( \bm{U}_{\mathsf{left}}, \bm{\Sigma}_{\mathsf{left}},
\bm{V}_{\mathsf{left}} \right)\;\longleftarrow\;\mathsf{svds} \left(
\bm{M}_{\mathsf{left}}, r \right),
\end{align*}
and use the top-$r$ column subspace $\bm{U}_{\mathsf{left}}$ as an
estimate of $\bm{U}^{\star}$ up to rotations. We partition $
\bm{U}_{\mathsf{left}}$ into two submatrices $ \bm{U}_{1}$ and $
\bm{U}_{2}$, where $ \bm{U}_{1} \in \mathbb{R}^{N_{1}\times r}$
consists of its top $N_{1}$ rows and $ \bm{U}_{2} \in
\mathbb{R}^{N_{2}\times r}$ consists of its bottom $N_{2}$ rows.
\item (\textbf{Matrix denoising}) Then we need an estimate of $
\bm{M}_{b}^{\star}$. The most natural idea is to use $ \bm{M}_{b}$,
which is an unbiased estimator for $ \bm{M}_{b}^{\star}$. However,
due to the noise in the observation, this naive estimate will incur
much larger estimation error, leading to suboptimal statistical
efficiency. Instead we compute the truncated rank-$r$ SVD of
$\bm{M}_{\mathsf{upper}}$, given by
\begin{align*}
\left( \bm{U}_{\mathsf{upper}}, \bm{\Sigma}_{\mathsf{upper}},
\bm{V}_{\mathsf{upper}} \right)\;\longleftarrow\;\mathsf{svds} \left(
\bm{M}_{\mathsf{upper}}, r \right),
\end{align*}
We partition $ \bm{V}_{\mathsf{upper}}$ into two submatrices
$\bm{V}_{1}$ and $\bm{V}_{2}$, where $ \bm{V}_{1} \in
\mathbb{R}^{T_{1}\times r}$ consists of its top $T_{1}$ rows and $
\bm{V}_{2} \in \mathbb{R}^{T_{2}\times r}$ consists of its bottom
$T_{2}$ rows. We use $\widehat{ \bm{M}}_{b} \mydefn
\bm{U}_{\mathsf{upper}} \bm{\Sigma}_{\mathsf{upper}}
\bm{V}_{2}^{\top}$ as an estimate of $ \bm{M}_{b}^{\star}$.
\item (\textbf{Linear regression}) For each $1 \leq t \leq T_{2}$, we
use the OLS solution to estimate the coefficient of the linear
combination of the $t$-th column of $ \bm{M}_{b}^{\star}$ with
respect to the basis $ \bm{U}_{1}^{\star}$
\begin{align*}
\widehat{ \bm{\beta}}_{t} \mydefn \mathop{\arg\min}_{ \bm{\beta} \in
\mathbb{R}^{r}}\big \Vert(\widehat{ \bm{M}}_{b})_{\cdot,
t}- \bm{U}_{1} \bm{\beta}\big
\Vert_{2}^{2} =( \bm{U}_{1}^{\top} \bm{U}_{1})^{-1} \bm{U}_{1}^{\top}(\widehat{ \bm{M}}_{b})_{\cdot,
t} \quad \text{for all }1 \leq t \leq T_{2}.
\end{align*}
This allows us to estimate the $t$-th column of $ \bm{M}_{d}^{\star}$
by $ \bm{U}_{2}\widehat{ \bm{\beta}}_{t}$. In matrix form, our final
estimate for $ \bm{M}_{d}^{\star}$ is $\widehat{ \bm{M}}_{d} \mydefn \bm{U}_{2}( \bm{U}_{1}^{\top} \bm{U}_{1})^{-1} \bm{U}_{1}^{\top}\widehat{ \bm{M}}_{b}$.
\end{itemize}
This procedure is summarized in~\Cref{alg:4-blocks-Md}.
As we will see in the following sections, under mild assumptions the
estimate $\widehat{ \bm{M}}_{d}$ achieves full statistical efficiency
for estimating $ \bm{M}_{d}^{\star}$ (even including the preconstant).
\subsection{Algorithm for the staggered adoption design}
\begin{figure}[t]
\centering
\begin{tabular}{cc}
\begin{tikzpicture}
\matrix (m) [matrix of math nodes,
nodes ={draw, minimum height =0.75cm, minimum width =1.1cm, anchor =center},
column sep =-\pgflinewidth, row sep =-\pgflinewidth]{
|[fill =red!20]| \bm{M}_{1, 1} & |[fill =red!20]| \bm{M}_{1, 2} & |[fill =blue!20]| \bm{M}_{1, 3} & |[fill =blue!20]| \bm{M}_{1, 4} & |[fill =gray!20]| \bm{M}_{1, 5} & |[fill =gray!20]| \bm{M}_{1, 6} \\
|[fill =red!20]| \bm{M}_{2, 1} & |[fill =red!20]| \bm{M}_{2, 2} & |[fill =blue!20]| \bm{M}_{2, 3} & |[fill =blue!20]| \bm{M}_{2, 4} & |[fill =gray!20]| \bm{M}_{2, 5} & \, \\
|[fill =red!20]| \bm{M}_{3, 1} & |[fill =red!20]| \bm{M}_{3, 2} & |[fill =blue!20]| \bm{M}_{3, 3} & |[fill =blue!20]| \bm{M}_{3, 4} & \, & \, \\
|[fill =yellow!20]| \bm{M}_{4, 1} & |[fill =yellow!20]| \bm{M}_{4, 2} & |[fill =gray!20]| \bm{M}_{4, 3} & \, & \, & \, \\
|[fill =yellow!20]| \bm{M}_{5, 1} & |[fill =yellow!20]| \bm{M}_{5, 2} & \, & \bm{?} & \, & \, \\
|[fill =gray!20]| \bm{M}_{6, 1} & \, & \, & \, & \, & \, \\
};
\foreach \n in {1, ..., 6} {
\node[above] at (m-1-\n.north) {$T_{\n}$};
}
\foreach \n in {1, ..., 6} {
\node[left] at (m-\n-1.west) {$N_{\n}$};
}
\draw[line width =0.5mm] (m-1-1.north west) rectangle (m-5-4.south east);
\draw[line width =0.5mm] (m-1-1.north west) rectangle (m-3-4.south east);
\draw[line width =0.5mm] (m-1-1.north west) rectangle (m-3-2.south east);
\draw[line width =0.5mm] (m-1-1.north west) rectangle (m-5-2.south east);
\end{tikzpicture} &
\begin{tikzpicture}
\matrix (m) [matrix of math nodes,
nodes ={draw, minimum height =0.75cm, minimum width =1.1cm, anchor =center},
column sep =-\pgflinewidth, row sep =-\pgflinewidth]{
|[fill =red!20]| \bm{M}_{1, 1} & |[fill =blue!20]| \bm{M}_{1, 2} & |[fill =blue!20]| \bm{M}_{1, 3} & |[fill =blue!20]| \bm{M}_{1, 4} & |[fill =blue!20]| \bm{M}_{1, 5} & |[fill =blue!20]| \bm{M}_{1, 6} \\
|[fill =yellow!20]| \bm{M}_{2, 1} & |[fill =gray!20]| \bm{M}_{2, 2} & |[fill =gray!20]| \bm{M}_{2, 3} & |[fill =gray!20]| \bm{M}_{2, 4} & |[fill =gray!20]| \bm{M}_{2, 5} & \, \\
|[fill =yellow!20]| \bm{M}_{3, 1} & |[fill =gray!20]| \bm{M}_{3, 2} & |[fill =gray!20]| \bm{M}_{3, 3} & |[fill =gray!20]| \bm{M}_{3, 4} & \, & \, \\
|[fill =yellow!20]| \bm{M}_{4, 1} & |[fill =gray!20]| \bm{M}_{4, 2} & |[fill =gray!20]| \bm{M}_{4, 3} & \, & \, & \, \\
|[fill =yellow!20]| \bm{M}_{5, 1} & |[fill =gray!20]| \bm{M}_{5, 2} & \, & \, & \, & \, \\
|[fill =yellow!20]| \bm{M}_{6, 1} & \, & \, & \, & \, & \bm{?} \\
};
\foreach \n in {1, ..., 6} {
\node[above] at (m-1-\n.north) {$T_{\n}$};
}
\foreach \n in {1, ..., 6} {
\node[left] at (m-\n-1.west) {$N_{\n}$};
}
\draw[line width =0.5mm] (m-1-1.north west) rectangle (m-6-6.south east);
\draw[line width =0.5mm] (m-1-1.north west) rectangle (m-6-1.south east);
\draw[line width =0.5mm] (m-1-1.north west) rectangle (m-1-6.south east);
\end{tikzpicture}
\tabularnewline
$ \quad \;\;$(a) & $ \quad \;\;\;$(b)\tabularnewline
\end{tabular}
\caption{Two examples on how to construct a four-block data matrix $ \bm{M}^{(i_0, j_0)}$ to estimate $ \bm{M}^\star_{i_0, j_0}$. We consider $k =6$, and the left and right panel correspond to $(i_0, j_0) =(5, 4)$ and $(6, 6)$ respectively. We use the question mark to denote the unobserved block that we want to estimate, and we use the bold line to single out the corresponding four-block design. For each observed block, we use different colors to distinguish their roles: red, blue and yellow blocks constitute $ \bm{M}_a$, $ \bm{M}_b$ and $ \bm{M}_c$ (cf.~(\ref{eq:reduction})), while gray blocks are those unimportant data that we discard while estimating $ \bm{M}^\star_{i_0, j_0}$.} \label{fig:alg-general}
\end{figure}
We move on to discuss the general staggered adoption design. At a high
level, our procedure is based on reducing to the four-block case.
Concretely, for any block $(i_{0}, j_{0})$ in~\Cref{fig:setup}, we
define an associated four-block estimation problem. The solution to
this problem---which can be obtained using the algorithm from the
preceding section---allows us to the block $\bm{M}_{i_{0},
j_{0}}^{\star}$. We obtain the four-block problem by removing a
subset of data so as to simplify the observation structure;
importantly, our theory to be described in the next section shows that
the left-out data is ultimately not of significant statistical
utility.
Let us now describe how to estimate a block $ \bm{M}_{i_{0},
j_{0}}^{\star}$, for some index $(i_{0}, j_{0})$ such that
\mbox{$j_{0} > k + 1-i_{0}$}, so that it is not observed. Defining the
indices $k_{1} =k + 1 - j_{0}$ and $k_{2} = k + 1 - i_{0}$, we
construct a four-block data matrix $\bm{M}^{(i_{0}, j_{0})}$ as
follows
\begin{align}
\label{eq:reduction}
\bm{M}^{(i_{0}, j_{0})} \mydefn \left[\begin{array}{cc} \bm{M}_{a} &
\bm{M}_{b} \\
\bm{M}_{c} & \bm{?}
\end{array} \right], \qquad \text{where} \quad
\begin{cases}
\bm{M}_{a} = \left[ \bm{M}_{i, j} \right]_{1 \leq i \leq k_{1}, 1
\leq j \leq k_{2}}, \\
\bm{M}_{b} = \left[ \bm{M}_{i, j} \right]_{1 \leq i \leq k_{1},
k_{2} < j \leq j_{0}}, \\
\bm{M}_{c} = \left[ \bm{M}_{i, j} \right]_{k_{1} < i \leq i_{0}, 1
\leq j \leq k_{2}}.
\end{cases}
\end{align}
Next we call~\Cref{alg:4-blocks-Md} with input $\bm{M}^{(i_{0},
j_{0})}$ to estimate its bottom-right block. Although the output
provides estimation jointly for all blocks $ \bm{M}_{i, j}^{\star}$
where $k_{1} < i \leq i_{0}$ and $k_{2} < j \leq j_{0}$, we only use
it to estimate the block $ \bm{M}_{i_{0}, j_{0}}^{\star}$. (For other
blocks, this estimate may not be statistically efficient.)
\Cref{fig:alg-general} and its caption provide two examples to help
illustrate how we construct the four-block data matrix for estimating
each unobserved block. The complete estimation procedure is summarized
in~\Cref{alg:general}.
As noted above, we show in the sequel that this algorithm achieves
full statistical efficiency (including the constant pre-factors).
Consequently, the data discarded while estimating each block are
indeed unimportant for estimating that block, due to the lack of
increase in statistical error.
\begin{algorithm}[t]
\DontPrintSemicolon \SetNoFillComment
\textbf{Input:} Data matrix $\bm{M}$, Boolean matrix $\bm{\Omega}$,
rank $r$. \\
Extract the dimension information $\{N_i\}_{1 \leq i \leq k}$ and
$\{T_j\}_{1 \leq j \leq k}$ from $ \bm{\Omega}$. \\
\For{$ i_0 = 1$ \KwTo $k$}{ \For{$j_0 = k + 1 - i_0$ \KwTo $k$}{
Construct the data matrix $ \bm{M}^{(i_0, j_0)}$ via
equation~\eqref{eq:reduction}. \\
Call~\Cref{alg:4-blocks-Md} with input $ \bm{M}^{(i_0, j_0)}$ and rank
$r$ to compute an estimate $\widehat{ \bm{M}}_d$ of its bottom left
block. \\
Extract the submatrix $\widehat{ \bm{M}}_{i_0, j_0}$ from $\widehat{
\bm{M}}_d$ corresponding to $\bm{M}_{i_0, j_0}^\star$. \\ } }
\textbf{Output:} $\widehat{ \bm{M}}_{i_0, j_0}$ as an estimate of
$\bm{M}_{i_0, j_0}^{\star}$ for each $(i_0, j_0)$
\caption{Estimating counterfactual outcomes: staggered adoption
design \label{alg:general}}
\end{algorithm}
\section{Main results}
\label{sec:main-results}
We are now positioned to present our main theoretical guarantees for
the algorithms described in~\Cref{sec:alg}. We begin with theory for
the four-block design in~\Cref{sec:theory-four}, before extending it
to the general case in~\Cref{sec:theory-general}.
\subsection{Theory for four-block design}
\label{sec:theory-four}
As noted previously, low-rank matrix recovery is an under-determined
problem, and certain regularity conditions are required to provide
theoretical guarantees. We begin by stating the conditions used in
our analysis; see the discussion in~\Cref{SecInterpret} for some
interpretation and intuition.
\paragraph{Sub-block conditioning:}
Our first set of conditions involve the top sub-blocks $\ensuremath{\bm{U}^\star}_1
\in\real^{N_1 \times r}$ and $\ensuremath{\bm{V}^\star}_1 \in \real^{T_1 \times r}$ of
$\ensuremath{\bm{U}^\star}$ and $\ensuremath{\bm{V}^\star}$, respectively. In particular, we assume that
there are constants $0 < \ensuremath{c_\ell} \leq \ensuremath{c_u} < \infty$ such that
\begin{align}
\label{EqnSubmatrix}
\ensuremath{c_\ell} \tfrac{\Nunit_{1}}{N} \bm{I}_{r} \preceq \bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star} \preceq \ensuremath{c_u} \tfrac{\Nunit_{1}}{N} \bm{I}_{r},
\qquad \text{and} \qquad \ensuremath{c_\ell} \tfrac{T_{1}}{T} \bm{I}_{r} \preceq
\bm{V}_{1}^{\star \top} \bm{V}_{1}^{\star} \preceq \ensuremath{c_u}
\tfrac{T_{1}}{T} \bm{I}_{r}.
\end{align}
\paragraph{Noise level:}
Throughout, so as to avoid degenerate corner cases, we assume that
$\min \{\Nunit_{1}, T_{1}\} \geq c \log (N+ T)$. In addition, some of
our results require an upper bound on the noise level, as measured by
the ratio $\ensuremath{\omega}/\sigma^*_r$. In particular, we define the noise
functional
\begin{align}
\label{eq:noise-level}
\ensuremath{\ensuremath{\rho}_{N,T}} \coloneqq \frac{\ensuremath{\omega}}{\sigma^*_r} \frac{1}{ \sqrt{ \min
\{\Nunit_1, \Time_1\}}}
\end{align}
For a given target error level $\delta > 0$, we require that $\min
\{\Nunit_1, \Time_1\}$ is sufficiently large so as to ensure that
\begin{align}
\label{eq:noise-condition-est}
\ensuremath{\ensuremath{\rho}_{N,T}} \; \sqrt{r+\log(N+T)} & \leq c_{\mathsf{noise}} \delta
\end{align}
for sufficiently small constant $c_{\mathsf{noise}} > 0$.
\paragraph{Local incoherence:}
For each $i \in [N]$ and $t \in [T]$, we define the local incoherence
parameters
\begin{subequations}
\label{eq:incoherence}
\begin{align}
\mu_{i} \coloneqq \sqrt{ \tfrac{N}{r} }\left\Vert \bm{U}_{i,
\cdot}^{\star} \right \Vert _{2} \qquad \text{and} \qquad \nu_{t}
= \sqrt{\tfrac{T}{r}} \left \Vert \bm{V}_{t,
\cdot}^{\star}\right\Vert _{2},
\end{align}
along with the noise levels $\ensuremath{\ensuremath{\rho}_N} \coloneqq
\frac{\sigma}{\sigma_{r}^{\star}} \; \sqrt{\frac{1}{ \Nunit_1 }}$ and
$\ensuremath{\ensuremath{\rho}_T} \coloneqq \frac{\sigma}{\sigma_{r}^{\star}} \; \sqrt{\frac{1}{
\Time_1 }}$. For any given target error level $\delta > 0$, in
order to provide recovery guarantees for entry $(i, t)$, we require
that
\begin{align}
\label{eq:incoherence-est}
\max \left\{ \mu_{i} \ensuremath{\ensuremath{\rho}_T}, \nu_{t} \ensuremath{\ensuremath{\rho}_N} \right\} \sqrt{r} \leq
\ensuremath{c_{\mathsf{inc}}} \delta \qquad \text{and} \qquad \min \left\{
\tfrac{\mu_{i}}{\sqrt{N_1}} , \tfrac{\nu_{t} }{\sqrt{T_1}} \right\}
\; \sqrt{r \log(N+T)} \leq \ensuremath{c_{\mathsf{inc}}} \delta
\end{align}
for some sufficiently small constant $\ensuremath{c_{\mathsf{inc}}} > 0$. For
proving proximity to a Gaussian variable (as stated in the main text),
we also require that
\begin{align}
\label{eq:signal-lb}
\min \left \{ \tfrac{\ensuremath{\ensuremath{\rho}_T}}{\mu_i}, \tfrac{\ensuremath{\ensuremath{\rho}_N}}{\nu_t} \right \} \; \Big
\{ \sqrt{\log(N+T)} + \tfrac{\log(N +T)}{\sqrt{r}} \Big \} & \leq
\ensuremath{c_{\mathsf{inc}}} \delta.
\end{align}
\end{subequations}
We refer the reader to~\Cref{lem:inference} for some analysis not
requiring this last condition.
\subsubsection{An achievable result and its consequences}
With these conditions in hand, we are now ready to state an achievable
result for our estimator. It provides a non-asymptotic bound, for any
unobserved index $(i,t)$, on the rescaled error of our estimator in
terms of a standard Gaussian variable $G_{i,t} \sim \ensuremath{\mathcal{N}}(0,1)$.
\begin{theorem}[\textsf{Non-asymptotics and distributional theory}]
\label{thm:distribution}
Under the sub-matrix condition~\eqref{EqnSubmatrix}, consider some
$\delta > 0$ for which the noise
condition~\eqref{eq:noise-condition-est} and incoherence
conditions~\eqref{eq:incoherence} hold. Then with probability at least
\mbox{$1 - O((N + T)^{-10})$,} \Cref{alg:4-blocks-Md} produces an
estimate $\widehat{\bm{M}}_d$ that for any unobserved index $(i,t)$,
we have
\begin{subequations}
\begin{align}
\label{EqnMainFourBound}
\Big| \tfrac{1}{ \smash[b]{\sqrt{\ensuremath{\ensuremath{\gamma}^*}_{i,t}}} }
\big(\widehat{\bm{M}}_{d} - \bm{M}_{d}^{\star}\big)_{i, t} - G_{i, t}
\Big| & \leq \delta,
\end{align}
and the variance takes the form
\begin{align}
\label{eq:variance-defn}
\ensuremath{\ensuremath{\gamma}^*}_{i, t} & \mydefn \frac{\ensuremath{\omega}^{2}}{N T} \Big \{ \bm{U}_{i,
\cdot}^{\star} (\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1}
\bm{U}_{i,\cdot}^{\star \top} + \bm{V}_{t, \cdot}^{\star}
(\bm{V}_{1}^{\star\top} \bm{V}_{1}^{\star})^{-1} \bm{V}_{t,
\cdot}^{\star\top} \Big \}.
\end{align}
\end{subequations}
\end{theorem}
Let us explore some consequences of~\Cref{thm:distribution}, beginning
with its uses in constructing confidence intervals. Observe that the
bound~\eqref{EqnMainFourBound} guarantees that the rescaled error is
close to a standard Gaussian variable. The scaling factor
$\ensuremath{\ensuremath{\gamma}^*}_{i,t}$, corresponding to the error variance, depends on
unknown population quantities. However, it can be estimated
accurately by standard plug-in approach, leading to data-driven
construction of entrywise confidence intervals. In particular, we
define
\begin{align*}
\ensuremath{\widehat{\ensuremath{\gamma}}}_{i,t} & \coloneqq \frac{\widehat{\ensuremath{\omega}}^{2}}{NT} \Big \{
\bm{U}_{i, \cdot}(\bm{U}_{1}^{\top} \bm{U}_{1})^{-1} \bm{U}_{i,
\cdot}^{\top} + \bm{V}_{t,\cdot} (\bm{V}_{1}^{\top} \bm{V}_{1})^{-1}
\bm{V}_{t, \cdot}^{\top} \Big \} \quad \text{where} \quad
\widehat{\ensuremath{\omega}}^{2} = \frac{N}{N_1} \big\Vert \bm{M}_{\mathsf{upper}}
- \widehat{\bm{M}}_{\mathsf{upper}} \big\Vert_{\mathrm{F}}^{2} .
\end{align*}
Here all of the matrices $\bm{U} = \bm{U}_{\mathsf{left}}$, $\bm{V} =
\bm{V}_{\mathsf{left}}$ and $\widehat{\bm{M}}_{\mathsf{upper}} =
\bm{U}_{\mathsf{upper}} \bm{\Sigma}_{\mathsf{upper}}
\bm{V}_{\mathsf{upper}}^\top$ can be computed from
\Cref{alg:4-blocks-Md}. Letting $\Phi$ denote the CDF of standard
Gaussian distribution, we then define the interval
\begin{subequations}
\begin{align}
\mathsf{CI}_{i, t}^{1-\alpha} \coloneqq \big[ \widehat{M}_{i,t} \pm
\Phi^{-1} (1-\alpha/2) \sqrt{\ensuremath{\widehat{\ensuremath{\gamma}}}_{i,t}} \big],
\end{align}
With this definition, we have the following guarantee:
\begin{proposition}
\label{prop:CI}
In addition to the conditions of \Cref{thm:distribution}, suppose that
$\min\{ N_1, T \} \geq \delta^{-2} r \log( N + T)$. Then the interval
$\mathsf{CI}_{i,t}^{1-\alpha}$ has the coverage guarantee
\begin{align}
\label{EqnCI}
\mathbb{P} \Big( \mathsf{CI}_{i,t}^{1-\alpha} \ni M_{i,t}^\star
\Big) = 1 - \alpha + O \big( \delta + (N+T)^{-10} \big).
\end{align}
\end{proposition}
\end{subequations}
\noindent See~\Cref{sec:proof-CI} for the proof. \\
Second, while~\Cref{thm:distribution} is stated in a way that
facilitates its inferential use, it also allows us to derive entrywise
error bounds for $\widehat{\bm{M}}_d$. If the conditions
of~\Cref{thm:distribution} hold with some $\delta$ decreasing to zero
as $\min \{N_1, T_1 \}$ grows, then it can be verified that
\begin{subequations}
\begin{align}
\mathbb{E}\big[(\widehat{M}_{i,t} - M_{i,t}^\star)^2 \big] = (1 +
o(1)) \; \ensuremath{\ensuremath{\gamma}^*}_{i, t}
\lesssim \frac{\ensuremath{\omega}^2}{N_1 T} \Vert \bm{U}_{i,\cdot}^\star \Vert_2^2
+ \frac{\ensuremath{\omega}^2}{N T_1} \Vert \bm{V}_{t,\cdot}^\star \Vert_2^2 .
\end{align}
To provide intuition for the convergence rate, suppose that the local
incoherence parameters $\mu_i$ and $\nu_t$ are viewed as quantities of
constant order; in this case, we have
\begin{align}
\mathbb{E}\big[(\widehat{M}_{i,t} - M_{i,t}^\star)^2 \big] \lesssim
\frac{\ensuremath{\omega}^2 r}{NT \min\{ N_1, T_1\}} \asymp \frac{r}{N T} \;
\ensuremath{\ensuremath{\rho}_{N,T}}^2.
\end{align}
\end{subequations}
In fact, even without the incoherence
conditions~\eqref{eq:incoherence} and with a relaxation of the noise
condition~\eqref{eq:noise-condition-est}, our analysis still provides
bounds on the estimation error; see~\Cref{lem:master} for details.
\subsubsection{Matching lower bound}
Is it possible to improve the variance term~\eqref{eq:variance-defn}
in our achievable result? In this section, we provide a negative
answer to this question by showing that it matches with a Bayesian
Cram\'er--Rao lower bound.
We decompose the unknown matrix as $\bm{M}^{\star} = \bm{X}^{\star}
\bm{Y}^{ \star \top }$, where $\bm{X}^{\star} \coloneqq
\bm{U}^{\star}( \bm{\Sigma}^{\star} )^{1/2}$ and $\bm{Y}^{\star}
\coloneqq \bm{V}^{\star}( \bm{\Sigma}^{\star} )^{1/2}$ are $N \times
r$ and $T \times r$ matrices, respectively. Now suppose that there is
a genie that reveals all rows of $\bm{X}^{\star}$ \emph{except} for
the $i$-th one $\xstar \coloneqq \bm{X}_{i, \cdot}^{\star}$, and also
all rows of $\bm{Y}^{\star}$ \emph{except} for its $t$-th one $\ystar
\coloneqq \bm{Y}_{t, \cdot}^{\star}$. Given this information, the
statistician is left with the following problem:
\begin{enumerate}
\item[(a)] The only unknown parameters are the $r$-dimensional row
vectors $\bm{x}^\star \coloneqq \bm{X}_{i, \cdot}^{\star}$ and
$\bm{y}^\star \coloneqq \bm{Y}_{t, \cdot}^{\star}$, and the goal is
to estimate their inner product $(\bm{M}_{d}^{\star})_{i, t} \equiv
\inprod{\bm{x}^\star}{\bm{y}^{\star}}$.
\item[(b)] The remaining observations of any use are the
$\Nunit_{1}+T_{1}$ linear measurements
\begin{align}
\label{eq:lb-samples}
\big\{ \bm{M}_{i, s} = \inprod{\xstar}{\bm{Y}_{s, \cdot}^{\star}} +
E_{i, s} \big\}_{s=1}^{T_1} \qquad \text{and} \quad \big\{ \bm{M}_{k,
t} = \inprod{\bm{X}_{k, \cdot}^{\star}}{\ystar} + E_{k, t}
\big\}_{k=1}^{\Nunit_1}.
\end{align}
All other observations are uninformative given the information
revealed by the genie.
\end{enumerate}
Thus, with the aid from the genie, the causal panel data model reduces
to a linear regression model with dimension $2 r$.
We use this genie-aided problem to compute lower bounds for the
estimation task---both the Cram\'er-Rao lower bound as well as a local
minimax risk. The local minimax risk at scale $\varepsilon > 0$ is
given by
\begin{align}
\label{eq:local-minimax-risk}
R_{\mathsf{loc}} (\varepsilon) \coloneqq \inf_{\widehat{M}_{i,t}} \;
\; \sup_{ \substack{\bm{x} \in \ensuremath{\mathcal{B}}_\infty (\bm{x}^\star ,
\varepsilon ) \\ \bm{y} \in \ensuremath{\mathcal{B}}_\infty (\bm{y}^\star ,
\varepsilon )} } \mathbb{E} \left[ \big (\widehat{M}_{i,t} -
\inprod{\bm{x}}{\bm{y}} \big)^2 \right].
\end{align}
Here the infimum is taken over all estimators of the scalar
$\inprod{\bm{x}}{\bm{y}}$ based on samples from the
model~\eqref{eq:lb-samples}, except that the true parameters are
$\bm{x}$ and $\bm{y}$ instead of $\bm{x}^\star$ and $\bm{y}^\star$.
\begin{theorem}[\textsf{Optimality}]
\label{thm:CRLB}
In the genie-aided setting, the Cram\'er-Rao lower bound for
estimating $(\bm{M}_d^\star)_{i, t}$ is given by
$\gamma_{i,t}^\star$. More generally, under the sub-block
condition~\eqref{EqnSubmatrix}, the local minimax risk at scale
$\varepsilon_{N,T} = \frac{\sqrt{\sigma_r^\star}}{\sqrt{ \min \{N, T
\} }}$ is lower bounded by
\begin{align}
\label{EqnLocalMinimax}
R_{\mathsf{loc}} ( \varepsilon_{N,T} ) \geq \Big( 1 -
\frac{\ensuremath{c_u} \pi^2}{\ensuremath{c_\ell}^2} \ensuremath{\ensuremath{\rho}_{N,T}}^2 \Big ) \; \gamma_{i,t}^\star,
\end{align}
where $\ensuremath{\ensuremath{\rho}_{N,T}}$ is the noise level previously
defined~\eqref{eq:noise-level}.
\end{theorem}
\noindent See~\Cref{sec:proof-thm-crlb} for the proof. \\
Under the noise assumption~\eqref{eq:noise-condition-est}, we have
$\ensuremath{\ensuremath{\rho}_{N,T}}^2 \leq \frac{c^2_{\mathsf{noise}} \delta^2}{r + \log(N +
T)}$, whence $\tfrac{R_{\mathsf{loc}}
(\varepsilon_{N,T})}{\gamma_{i,t}^\star} \geq 1 - O(\delta^2)$ as
$\delta$ goes to zero. Thus, we see that our our estimator
(cf.~\Cref{alg:4-blocks-Md}) matches the local minimax risk aided by
the genie, including a sharp pre-factor. Thus, our procedure is
optimal in a strong instance-dependent sense.
\subsubsection{Interpretation of structural conditions}
\label{SecInterpret}
Let us interpret the structural conditions---namely, the sub-block and
local incoherence conditions---that underlie~\Cref{thm:distribution}.
\paragraph{Sub-block condition:}
Since $\bm{M}^{\star}_a = \ensuremath{\bm{U}^\star}_1 \ensuremath{\bm{\Sigma}^\star} \ensuremath{\bm{V}^\star}_1$, the existence
of the constant $\ensuremath{c_\ell} > 0$ ensures that $\MatTre{a}$ has full rank.
This condition is necessary, because otherwise it would be impossible
to recover the missing block $\MatTre{d}$ even in the noiseless
setting.\footnote{For example, if the row rank of $\MatTre{a}$ is less
than $r$, then we may find a row of $\MatTre{b}$ that is not in the
subspace spanned by the rows of $\MatTre{a}$. In this case, the
corresponding row of $\MatTre{d}$ is unidentifiable.}. Otherwise, we
observe that $\ensuremath{\bm{U}^\star}_1$ is an $N_1 \times T_1$ sub-matrix of the
orthonormal matrix $\ensuremath{\bm{U}^\star} \in \real^{N \times T}$. Since $\ensuremath{\bm{U}^\star}$
has orthonormal columns, we expect that the typical size of each entry
is $1/\sqrt{N}$. Thus, the typical squared norm of each column of the
sub-matrix $\ensuremath{\bm{U}^\star}_1$ is of the order $N_1/N$, which accounts for the
normalization in the condition~\eqref{EqnSubmatrix} on $\ensuremath{\bm{U}^\star}_1$.
Similar comments apply to the condition on $\ensuremath{\bm{V}^\star}_1$.
\paragraph{Local incoherence conditions:}
Since the Frobenius norm of $\bm{U}^{\star}$ and $\bm{V}^{\star}$ are
both $\sqrt{r}$, the ``typical size'' of their rows (in $\ell_{2}$
norm) are $\sqrt{r/N}$ and $\sqrt{r/T}$ respectively. We can see that
$\mu_{i}$ (resp.~$\nu_{i}$) evaluates the ratio between $\Vert
\bm{U}_{i, \cdot}^{\star} \Vert_{2}$ (resp.~$\Vert \bm{V}_{t,
\cdot}^{\star} \Vert_{2}$) and its typical size. This definition
connects to the (global) incoherence condition that is widely assumed
in the low-rank matrix completion literature~\cite{ExactMC09,
chen2015incoherence, chi2018nonconvex}, where the (global)
incoherence parameters $\mu \coloneqq\max_{1 \leq i \leq N} \mu_{i}$
and $\nu \coloneqq \max_{1\leq t\leq T}\nu_{t}$ defined therein cannot
be too large in order for the recovery to be possible. In contrast,
our theory does not assume any global incoherence condition to hold.
If our goal is to estimate and infer the counterfactual outcome or
treatment effect of unit $i$ at time $t$, we just need the local
incoherence parameters $\mu_{i}$ and $\nu_{t}$ to be reasonably
well-behaved.
Our constraints involving the local incoherence
conditions~\eqref{eq:incoherence-est} are reasonably mild. When
$\min\{\Nunit_{1}, T_{1}\}\gg r$ and the noise condition
(\ref{eq:noise-condition-est}) is satisfied, we can allow both
$\Vert\bm{U}_{i, \cdot}^{\star}\Vert_{2}$ and $\Vert \bm{V} _ {t ,
\cdot}^{\star} \Vert_{2}$ to be much larger than their typical
size. The other condition \eqref{eq:signal-lb} requires one of
$\Vert\bm{U}_{i, \cdot}^{\star}\Vert_{2}$ and $\Vert\bm{V}_{t,
\cdot}^{\star}\Vert_{2}$ to be not vanishingly small compared to
their typical size. In fact, even when (\ref{eq:signal-lb}) does not
hold, we can still characterize the entrywise distribution of
$\widehat{\bm{M}}_{d}-\bm{M}_{d}^{\star}$, although it is not
approximately Gaussian; we refer interested readers to
\Cref{lem:inference} for this general result.
\subsection{Theory for staggered adoption design}
\label{sec:theory-general}
In this section, we extend our theoretical analysis to the general
(non-four-block) setting to which \Cref{alg:general} applies. Without
loss of generality, we focus on estimating and inferring entries of an
unobserved submatrix $\bm{M}_{i_{0}, j_{0}}^{\star}$. In
\Cref{sec:alg}, we described how to construct a data matrix
$\bm{M}^{(i_{0}, j_{0})}$ for estimating $\bm{M}_{i_{0},
j_{0}}^{\star}$ that admits the four-block design; see
equation~\eqref{eq:reduction} for details. \Cref{fig:theory} provides
a simple example so as to help visualize the set-up and to facilitate
later discussion.
At the left of~\Cref{fig:theory}, we depict a panel data problem with
units separated into $k = 6$ subsets, as indicated by the large
blocks. The $i$-th subset has $N_i$ units that all received the
treatment at time $1+\sum_{j=1}^{i} T_j$; recall that a time larger
than $T$ means that the unit never received the treatment. All the
colored blocks are observed, while the uncolored blocks are
missing. Suppose we want to estimate the sub-matrix
$\bm{M}_{i_0,j_0}^\star$ indexed by the pair $(i_0, j_0)=(5,4)$, as
indicated by a the question mark in the figure. The data matrix
$\bm{M}_{(i_{0}, j_{0})}^\star$ consists of the $i_0=5$ by $j_0=4$
block matrix on the top left corner of $\bm{M}$. The gray blocks,
although observed, are discarded to enforce the four-block
structure. We use three colors (red, blue and yellow) to show the
sub-matrices that play the role of $\bm{M}_a^\star$, $\bm{M}_b^\star$
and $\bm{M}_c^\star$ in the four-block model. The red component
consists of the $k_1=3$ by $k_2=2$ block matrix on the top left
corner; in general, these two integers satisfy the relations $k_1 = k
+ 1 - j_0$ and $k_2 = k + 1 - i_0$.
Now we define the dimension parameters $\bar{N}_{1}$, $\bar{N}_{2}$,
$\bar{T}_{1}$ and $\bar{T}_{2}$ with respect to the four-block
structure, given by
\begin{align*}
\bar{N}_{1} \coloneqq \sum_{i=1}^{k_{1}}N_{i}, \quad\bar{N}_{2}
\coloneqq \sum_{i=k_{1}+1}^{i_{0}}N_{i}, \quad\bar{T}_{1} \coloneqq
\sum_{j=1}^{k_{2}} T_{j}, \quad \mbox{and} \quad \bar{T}_{2} \coloneqq
\sum_{j=k_{2}+1}^{j_{0}} T_{j},
\end{align*}
We also define $\bar{N} \coloneqq \bar{N}_{1}+\bar{N}_{2}$ and
$\bar{T} \coloneqq \bar{T}_{1}+\bar{T}_{2}$. Let $\bm{U}_{1}^{\star}
\in \mathbb{R}^{\bar{N}_{1}\times r}$ and $\bm{U}_{2}^{\star} \in
\mathbb{R}^{\bar{N}_{2}\times r}$ be submatrices of $\bm{U}^{\star}
\in \mathbb{R}^{N\times r}$, consisting of its first $\bar{N}_{1}$
rows and the next $\bar{N}_{2}$ rows respectively. Similarly, let
$\bm{V}_{1}^{\star} \in \mathbb{R}^{\bar{T}_{1}\times r}$ and
$\bm{V}_{2}^{\star} \in \mathbb{R}^{\bar{T}_{2}\times r}$ be
submatrices of $\bm{V}^{\star} \in \mathbb{R}^{T\times r}$, consisting
of its first $\bar{T}_{1}$ rows and the next $\bar{T}_{2}$ rows
respectively. These dimension parameters are also illustrated
in~\Cref{fig:theory}.
\begin{figure}[h]
\centering \begin{tabular}{ccc}
\begin{tikzpicture}
\matrix (m) [matrix of math nodes,
nodes={draw, minimum height=0.75cm, minimum width=0.9cm, anchor=center, font=\small},
column sep=-\pgflinewidth, row sep=-\pgflinewidth]{
|[fill=red!20]| \, & |[fill=red!20]| \, & |[fill=blue!20]| \, & |[fill=blue!20]| \, & |[fill=gray!20]| \, & |[fill=gray!20]| \, \\
|[fill=red!20]| \, & |[fill=red!20]| \, & |[fill=blue!20]| \, & |[fill=blue!20]| \, & |[fill=gray!20]| \, & \, \\
|[fill=red!20]| \, & |[fill=red!20]| \, & |[fill=blue!20]| \, & |[fill=blue!20]| \, & \, & \, \\
|[fill=yellow!20]| \, & |[fill=yellow!20]| \, & |[fill=gray!20]| \, & \, & \, & \, \\
|[fill=yellow!20]| \, & |[fill=yellow!20]| \, & \, & \bm{?} & \, & \, \\
|[fill=gray!20]| \, & \, & \, & \, & \, & \, \\
};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt}]
(m-1-1.north west) -- (m-1-2.north east) node [black, midway, yshift=13pt] {$\bar{T}_1$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt}]
(m-1-3.north west) -- (m-1-4.north east) node [black, midway, yshift=13pt] {$\bar{T}_2$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}]
(m-6-1.south west) -- (m-6-4.south east) node [black, midway, yshift=-13pt] {$\bar{T}$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}]
(m-1-1.north west) -- (m-3-1.south west) node [black, midway, xshift=-15pt, yshift=-8pt, anchor=south] {$\bar{N}_1$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}]
(m-4-1.north west) -- (m-5-1.south west) node [black, midway, xshift=-15pt, yshift=-8pt, anchor=south] {$\bar{N}_2$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt}]
(m-1-6.north east) -- (m-5-6.south east) node [black, midway, xshift=13pt, yshift=-8pt, anchor=south] {$\bar{N}$};
\draw[line width=0.5mm] (m-1-1.north west) rectangle (m-5-4.south east);
\draw[line width=0.5mm] (m-1-1.north west) rectangle (m-3-4.south east);
\draw[line width=0.5mm] (m-1-1.north west) rectangle (m-3-2.south east);
\draw[line width=0.5mm] (m-1-1.north west) rectangle (m-5-2.south east);
\end{tikzpicture}
&
\begin{adjustbox}{valign=t, raise=56mm}
\begin{tikzpicture}
\matrix (m) [matrix of math nodes,
nodes={draw, minimum width=0.7cm, anchor=center, font=\small},
column sep=-\pgflinewidth, row sep=-\pgflinewidth]{
|[minimum height=3*0.75cm, fill=red!20]| \bm{U}_{1}^\star \\
|[minimum height=2*0.75cm, fill=yellow!20]| \bm{U}_{2}^\star \\
|[minimum height=0.75cm, fill=gray!20]| \, \\
};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt}]
(m-1-1.north west) -- (m-1-1.north east) node [black, midway, yshift=13pt] {$r$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}]
(m-1-1.north west) -- (m-1-1.south west) node [black, midway, xshift=-15pt, yshift=-8pt, anchor=south] {$\bar{N}_1$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}]
(m-2-1.north west) -- (m-2-1.south west) node [black, midway, xshift=-15pt, yshift=-8pt, anchor=south] {$\bar{N}_2$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt}]
(m-1-1.north east) -- (m-2-1.south east) node [black, midway, xshift=13pt, yshift=-8pt, anchor=south] {$\bar{N}$};
\draw[line width=0.5mm] (m-1-1.north west) rectangle (m-1-1.south east);
\draw[line width=0.5mm] (m-1-1.north west) rectangle (m-2-1.south east);
\end{tikzpicture}
\end{adjustbox}
&
\begin{adjustbox}{valign=t, raise=37mm}
\begin{tikzpicture}
\matrix (m) [matrix of math nodes,
nodes={draw, minimum height=0.7cm, anchor=center, font=\small},
column sep=-\pgflinewidth, row sep=-\pgflinewidth]{
|[minimum width=2*0.9cm, fill=red!20]| \bm{V}_{1}^{\star\top} &
|[minimum width=2*0.9cm, fill=blue!20]| \bm{V}_{2}^{\star\top} &
|[minimum width=2*0.9cm, fill=gray!20]| \, \\
};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}]
(m-1-1.north west) -- (m-1-1.south west) node [black, midway, xshift=-11pt, yshift=-6pt, anchor=south] {$r$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt}]
(m-1-1.north west) -- (m-1-1.north east) node [black, midway, yshift=13pt] {$\bar{T}_1$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt}]
(m-1-2.north west) -- (m-1-2.north east) node [black, midway, yshift=13pt] {$\bar{T}_2$};
\draw[decorate, decoration={brace, amplitude=6pt, raise=1pt, mirror}]
(m-1-1.south west) -- (m-1-2.south east) node [black, midway, yshift=-13pt] {$\bar{T}$};
\draw[line width=0.5mm] (m-1-1.north west) rectangle (m-1-1.south east);
\draw[line width=0.5mm] (m-1-1.north west) rectangle (m-1-2.south east);
\end{tikzpicture}
\end{adjustbox}
\tabularnewline
$\;\;$ $\bm{M}$ & $\;\;$ $\bm{U}^\star$ & $\quad\quad$ $\bm{V}^{\star\top}$ \tabularnewline
\end{tabular}
\caption{An illustration of the dimension parameters and the
partition used in~\Cref{sec:theory-general}. }
\label{fig:theory}
\end{figure}
In analogy with the four-block design, we first introduce the
conditions required for establishing theoretical guarantees in the
general staggered adoption design.
\paragraph{Sub-block conditioning:}
We assume that there are constants $0 < \ensuremath{c_\ell} \leq \ensuremath{c_u} < \infty$
such that
\begin{subequations}
\label{EqnSubmatrix-general}
\begin{align}
\ensuremath{c_\ell} \frac{\bar{N}_{1}}{N} \bm{I}_{r} \preceq \bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star} \preceq \ensuremath{c_u} \frac{\bar{N}_{1}}{N} \bm{I}_{r}, &
\qquad \ensuremath{c_\ell} \frac{\bar{N}}{N} \bm{I}_{r} \preceq
\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star} + \bm{U}_{2}^{\star\top}
\bm{U}_{2}^{\star} \preceq \ensuremath{c_u} \frac{\bar{N}}{N} \bm{I}_{r}, \\
\ensuremath{c_\ell} \frac{\bar{T}_{1}}{T} \bm{I}_{r} \preceq \bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star} \preceq \ensuremath{c_u} \frac{\bar{T}_{1}}{T} \bm{I}_{r}, &
\qquad \ensuremath{c_\ell} \frac{\bar{T}}{T} \bm{I}_{r} \preceq
\bm{V}_{1}^{\star\top} \bm{V}_{1}^{\star} + \bm{V}_{2}^{\star\top}
\bm{V}_{2}^{\star} \preceq \ensuremath{c_u} \frac{\bar{T}}{T} \bm{I}_{r}.
\end{align}
\end{subequations}
Recall the sub-block condition \eqref{EqnSubmatrix} for the four-block
design and the discussion in \Cref{SecInterpret}, here we are assuming
in addition that $\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star} +
\bm{U}_{2}^{\star\top} \bm{U}_{2}^{\star}$ and $\bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star} + \bm{V}_{2}^{\star\top} \bm{V}_{2}^{\star}$ do not
deviate too much from their typical size. This additional requirement
connects the spectrum of the ground truth of $\bm{M}^{(i_{0}, j_{0})}$
(which is a submatrix of $\bm{M}^{\star}$) with that of the full
ground truth $\bm{M}^{\star}$, which allows us to present cleaner
results that depend on the spectrum of $\bm{M}^{\star}$ instead of its
submatrices.
\paragraph{Noise level:}
The noise level in this general setting takes the following form:
\begin{align}
\label{eq:noise-level-general}
\ensuremath{\ensuremath{\bar{\rho}}_{N,T}} \coloneqq \max \{ \ensuremath{\ensuremath{\bar{\rho}}_N} , \ensuremath{\ensuremath{\bar{\rho}}_T} \}, \qquad \text{where}
\qquad \ensuremath{\ensuremath{\bar{\rho}}_N} \coloneqq \frac{\ensuremath{\omega}}{\sigma_{r}^{\star}} \;
\sqrt{\frac{1}{ \bar{\Nunit}_1 }} \qquad \text{and} \qquad \ensuremath{\ensuremath{\bar{\rho}}_T}
\coloneqq \frac{\ensuremath{\omega}}{\sigma_{r}^{\star}} \; \sqrt{\frac{1}{
\bar{\Time}_1 }}.
\end{align}
For any prescribed target error level $\delta > 0$, we assume that it
is upper bounded by
\begin{align}
\label{eq:noise-condition-est-general}
\ensuremath{\ensuremath{\bar{\rho}}_{N,T}} \sqrt{r+\log(N+T)} \leq c_{\mathsf{noise}} \delta ,
\end{align}
for some sufficiently small constant $c_{\mathsf{noise}} > 0$.
\paragraph{Local incoherence:}
For each $i \in [N]$ and $t \in [T]$, the local incoherence parameters
$\mu_i$ and $\nu_t$ are defined as in equation~\eqref{eq:incoherence}.
For any given target error level $\delta>0$, we assume that there is
some sufficiently small constant $\ensuremath{c_{\mathsf{inc}}}>0$ such that,
\begin{subequations}
\label{eq:incoherence-general}
\begin{align}
\label{eq:incoherence-est-general}
\max \left\{ \mu_{i} \ensuremath{\ensuremath{\bar{\rho}}_T}, \nu_{t} \ensuremath{\ensuremath{\bar{\rho}}_N} \right\} \sqrt{r} \leq
\ensuremath{c_{\mathsf{inc}}} \delta \qquad \text{and} \qquad \min \left\{
\tfrac{\mu_{i}}{\sqrt{\bar{N}_1}} , \tfrac{\nu_{t}
}{\sqrt{\bar{T}_1}} \right\} \; \sqrt{r \log(N+T)} \leq
\ensuremath{c_{\mathsf{inc}}} \delta
\end{align}
and that at least one of the following conditions holds
\begin{align}
\label{eq:signal-lb-1}
\min \left \{ \tfrac{\ensuremath{\ensuremath{\bar{\rho}}_T}}{\mu_i}, \tfrac{\ensuremath{\ensuremath{\bar{\rho}}_N}}{\nu_t} \right \}
\; \Big \{ \sqrt{\log(N+T)} + \tfrac{\log(N +T)}{\sqrt{r}} \Big \} &
\leq \ensuremath{c_{\mathsf{inc}}} \delta.
\end{align}
\end{subequations}
With this set-up, we let us describe how our theory applies to this
more general setting, in particular by showing proximity of the
rescaled error to a standard normal $G_{i,t} \sim \ensuremath{\mathcal{N}}(0,1)$.
\begin{theorem}
[\textsf{Non-asymptotics and distribution: General setting}]
\label{thm:distribution-general}
For a given $\delta > 0$, suppose that the noise
condition~\eqref{eq:noise-condition-est-general} and the incoherence
conditions~\eqref{eq:incoherence-general} hold. For any index $(i,t)$
belonging to any unobserved submatrix $\bm{M}_{i_0,j_0}^\star$, with
probability at least \mbox{$1 - O((N + T)^{-10})$,} \Cref{alg:general}
produces an estimate $\widehat{\bm{M}}_{i_0,j_0}$ such that
\begin{subequations}
\begin{align}
\label{EqnMainStaggeredBound}
\Big| \tfrac{1}{ \smash[b]{\sqrt{\ensuremath{\ensuremath{\gamma}^*}_{i,t}}} }
\big(\widehat{\bm{M}}_{i_0,j_0} - \bm{M}_{i_0,j_0}^{\star}\big)_{i, t}
- G_{i, t} \Big| & \leq \delta,
\end{align}
where the variance is given by
\begin{align}
\label{eq:variance-defn-general}
\ensuremath{\ensuremath{\gamma}^*}_{i, t} & \mydefn \frac{\ensuremath{\omega}^{2}}{NT} \Big \{ \bm{U}_{i,
\cdot}^{\star} (\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1}
\bm{U}_{i}^{\star\top} + \bm{V}_{t}^{\star} (\bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1} \bm{V}_{t, \cdot}^{\star\top} \Big \}.
\end{align}
\end{subequations}
\end{theorem}
\noindent
See~\Cref{sec:proof-general-distribution} for the proof. \\
Observe that this result is the natural generalization of our previous
achievable result (\Cref{thm:distribution}) for the simpler four-block
setting. This correspondence makes sense, since our algorithm for the
general case involves ``reducing'' the problem to instances of the
simpler four-block case. Of course, one might wonder whether or not
this reduction---which involves discarding some data that might be
relevant---leads to a sub-optimal guarantee.
Interestingly, our lower bound---establishing the optimality of the
variance $\gamma^\star_{i,t}$, as stated in~\Cref{thm:CRLB}---also
extends to this multi-block setting. (We do not state this lower
bound formally here, but note that our proof is actually given for the
multi-block setting; see~\Cref{sec:proof-thm-crlb} for details.) As a
consequence, under the conditions of~\Cref{thm:distribution-general},
there is \emph{no appreciable loss} of accuracy in solving the problem
by successive reduction to the four-block setting.
\subsection{Comparison with previous work}
\label{sec:comparison}
Having given precise statements of our results, let us now provide a
detailed comparison with related results established in past
work~\cite{athey2021matrix,bai2021matrix,agarwal2023causal,choi2023matrix}.
To simplify discussion, we ignore any log factors in making these
comparisons, and use $\lesssim$ to reflect inequalities that ignore
constants and such log factors. For making comparisons, we record
here the estimation error bound that can can be deduced
from~\Cref{thm:distribution-general}: for any entry $(i,t)$, the error
is upper bounded as
\begin{align}
\label{eq:our-bound}
\big | \widehat{M}_{i , t} - M_{i , t}^\star \big | \lesssim \ensuremath{\omega}
\sqrt{ \frac{1}{\bar{N}_1 T} } \Vert \bm{U}_{i, \cdot}^\star \Vert_2
+ \ensuremath{\omega} \sqrt{ \frac{1}{N \bar{T}_1} } \Vert \bm{V}_{t,
\cdot}^\star \Vert_2
\end{align}
with high probability.
\begin{itemize}
\item Athey et al.~\citep{athey2021matrix} proposed to use
convex relaxation (nuclear norm minimization) to estimate
the missing potential outcomes. They showed that with high
probability
\begin{equation}
\label{eq:athey-bound}
\frac{1}{ \sqrt{ N T } } \big \Vert \widehat{ \bm{M}
}^{\mathsf{ABD}} - \bm{M}^\star \big
\Vert_{\mathrm{F}} \lesssim \frac{ N^{1/4} \Vert
\bm{M}^\star \Vert_\infty }{ \sqrt{N_1} } +
\frac{ \ensuremath{\omega} }{ N_1 T } \sqrt{\frac{ r }{ \min \{ N, T \} }
} .
\end{equation}
This is a Frobenius norm bound, which measures the average error over
the full matrix and does not provide entrywise guarantees like
(\ref{eq:our-bound}). In addition, their bound (\ref{eq:athey-bound})
does not vanish as the noise level $\sigma$ goes to zero, so that even
when there is no noise, it does not ensure exact recovery.
\item In recent work, Choi and Yuan~\citep{choi2023matrix} built upon
the paper~\cite{athey2021matrix}, providing an algorithm with
superior estimation guarantees. In particular, they suggested a
more sophisticated and multi-round use of convex relaxation, based
on dividing the missing entries into smaller groups, and then
solving a nuclear-norm regularized problem for each of these groups.
They showed that with high probability, their estimator
$\widehat{\bm{M}}^{\mathsf{CY}}$ satisfies the entrywise
$\ell_{\infty}$-bound
\begin{align}
\label{eq:choi-bound}
\big \Vert \widehat{ \bm{M} }^{\mathsf{CY}} - \bm{M}^\star \big
\Vert_{\infty} \lesssim \frac{\ensuremath{\omega}}{\sqrt{NT}} \sqrt{\frac{ \kappa^5 \mu^2 r^3 }{
\min \{ \bar{N}_1, \bar{T}_1 \} } }.
\end{align}
Here $\kappa=\sigma_1^\star/\sigma_r^\star$ is the condition number of
$\bm{M}^\star$, and $\mu$ is the \emph{global} incoherence parameter
(i.e., the maximum over all the local parameters $\mu_i$ and $\nu_t$
defined in equation~\eqref{eq:incoherence}). In comparison, our
bound~\eqref{eq:our-bound}) is a localized improvement of their
bound~\eqref{eq:choi-bound}; in particular, it has better dependency
on the rank $r$, and no dependence on the condition number. In
addition, we also provide data-driven confidence intervals for each
entry of $\bm{M}^\star$ with optimal width. It is also worth
mentioning that our spectral algorithm is more computationally
efficient than the SDP-based algorithm in \cite{choi2023matrix}.
\item Agarwal et al.~\citep{agarwal2023causal} proposed a procedure
called ``synthetic nearest neighbors'' (SNN), for which they proved
entrywise error bound and established asymptotic normality. Their
algorithm works for a wider range of missing mechanism and is more
general than our algorithm. However, their error rate is not optimal
under the staggered adoption design. For example, if $N=T$ and there
is only one missing entry, their estimation error scales as
$N^{-1/4}$, which is slower than the $N^{-1/2}$ rate for our
algorithm and convex relaxation.
\item Another line of work~\citep{bai2021matrix,cahan2023factor} is
based on assuming a factor structure for the full panel data matrix,
and applying spectral algorithms that share similar spirit with
ours. In general, the analyses in these papers impose stronger
assumptions than those used in our analysis. First, they model the
factors and factor loadings---roughly speaking, these correspond to
the rows of $\bm{U}^\star$ and $\bm{V}^\star$--- as being random,
and required to satisfy certain asymptotic moment conditions.
Second, they require that the observation sizes satisfy the
constraint $\min\{ \bar{N}_1, \bar{T}_1 \} \gg \sqrt{ \max \{ N, T
\}}$, as well as certain eigengap conditions. Neither of these
conditions are imposed in our analysis. Third, their underlying
noise condition (used in the asymptotic analysis) is
$\sqrt{\min\{\bar{N}_1, \bar{T}_1\}}$ times more stringent than our
noise condition~\eqref{eq:noise-condition-est}. On the other hand,
we note that their set-up is general enough to handle weakly
correlated noise, which is not covered by our framework.
\end{itemize}
\section{Numerical experiments}
\label{sec:numerical}
In this section, we report the results of some numerical experiments
to corroborate the computational simplicity of our algorithm, and to
study the sharpness of our theoretical predictions. For certain
random ensembles, we draw matrices $\bm{U} \in \mathbb{R}^{N\times r}$
from the Haar distribution over the Grassmann manifold; for
simplicity, we say that any such $\bm{U}$ is a random orthonormal
matrix.
\subsection{Scaling with incoherence}
All experiments in this section were based on matrices with dimensions
$N = T = 500$, rank $r = 3$, and using the staggered adoption design
with $k = 3$ groups with dimension parameters $N_{1} = T_{1} = 200$,
$N_2 = T_2 = 200$ and $N_3 = T_3 = 100$. We used the fixed noise
level $\omega = 1$.
Our goal in this set of experiments was to examine the necessity and
role of incoherence
conditions~\eqref{eq:incoherence-general}. Focusing on estimation of
the entry $(i, t) =(500, 500)$ for concreteness, we generated random
problems of the following type. For an incoherence parameter $\mu \in
(0, \sqrt{N/r})$, we generated an orthonormal matrix as
\begin{subequations}
\label{eq:ensemble}
\begin{align}
\label{eq:U-mu}
\bm{U}^\star (\mu) = \left[\begin{matrix} \sqrt{1 - \mu^2 \frac{r}{N}}
\, \bm{U}^{-} \\ \mu \sqrt{ \frac{r}{N}} \, \bm{O}_1
\end{matrix}\right]
\end{align}
where $\bm{U}^{-} \in \mathbb{R}^{(N-r) \times r}$ and $\bm{O}_1 \in
\mathbb{R}^{r\times r}$ are random orthonormal matrices. For any $\mu
\in (0,\sqrt{N/r})$, it is straightforward to check that $\bm{U}^\star
(\mu) \in \mathbb{R}^{N\times r}$ is an orthonormal matrix, and it has
associated local incoherence parameter $\mu_i \equiv \sqrt{N/r} \Vert
\bm{e}_i^\top \bm{U}^\star (\mu) \Vert_2 = \mu$. Similarly, for any
$\nu \in (0,\sqrt{T/r})$, we can construct an orthonormal matrix
$\bm{V}^\star (\nu)$ with a given $\nu_t = \nu$ by
\begin{align}
\label{eq:V-nu}
\bm{V}^\star(\nu) = \left[\begin{matrix} \sqrt{1 - \nu^2 \frac{r}{T}}
\, \bm{V}^{-} \\ \nu \sqrt{ \frac{r}{T}} \, \bm{O}_2
\end{matrix}\right]
\end{align}
\end{subequations}
where $\bm{V}^{-} \in \mathbb{R}^{(T-r) \times r}$ and $\bm{O}_2 \in
\mathbb{R}^{r\times r}$ are random orthonormal matrices. We consider
two cases:
\begin{itemize}
\item {\bf Case 1:} We generate $\bm{M}^{\star}=\bm{U}^{\star} (\mu)
\bm{V}^{\star\top}$ where $\bm{U}^\star(\mu) \in \mathbb{R}^{N\times
r}$ is defined in equation~\eqref{eq:U-mu} with
$\mu$ varying between $10^{-4}$ and $10^2$, while $\bm{V}^\star \in
\mathbb{R}^{T\times r}$ is a random orthonormal matrix.
\item {\bf Case 2:} We generate $\bm{M}^{\star}=\bm{U}^{\star} (\mu)
\bm{V}^{\star}(\nu)^\top$ where $\bm{U}^\star (\mu) \in
\mathbb{R}^{N\times r}$ and $\bm{V}^\star (\nu) \in \mathbb{R}^{T
\times r}$ are defined in equations~\eqref{eq:U-mu}
and~\eqref{eq:V-nu}, respectively, with $\mu=\nu$ varying between
$10^{-4}$ and $10^2$.
\end{itemize}
\begin{figure}[h]
\begin{center}
\begin{tabular}{ccc}
\widgraph{0.45\textwidth}{figs/fig_incoherence_case1} &&
\widgraph{0.45\textwidth}{figs/fig_incoherence_case2} \\
(a) && (b)
\end{tabular}
\caption{Red circles: Mean-squared error (MSE) of the estimate
$\widehat{M}_{i,t}$ of entry $(i,t) = (500, 500)$ returned by
Algorithm~\ref{alg:general} versus the the incoherence parameters.
Blue lines: theoretically predicted scaling $\gamma_{i,t}^\star$
of the MSE versus incoherence. Empirical MSE was obtained as an
average of $100$ Monte Carlo trials. (a) Results for the matrix
ensemble in Case 1. We see good agreement between simulation and
theory across the full range of incoherence. (b) Results for Case
2. Simulation/theory agreement is good for incoherence parameters
above $10^{-2}$. Theory breaks down below this level because
higher-order terms start to dominate the MSE for a fixed matrix
dimension.}
\label{fig:incoherence}
\end{center}
\end{figure}
For each case, we performed a total of $100$ Monte Carlo trials, each
time computing our estimate and forming a Monte Carlo approximation of
the mean-squared error
$\mathbb{E}[(\widehat{M}_{i,t}-M_{i,t}^\star)^2]$ associated with
estimating entry $(i,t)$ of $M_{i,t}^\star$. Panels (a) and (b),
respectively, of~\Cref{fig:incoherence} give plots of these empirical
MSEs versus the incoherence parameter for Cases 1 and 2, respectively.
For comparison, we also plot the theoretically predicted variance
$\gamma^\star_{i,t}$ from equation~\eqref{eq:variance-defn} for these
two cases. For Case 1 in panel (a), we see excellent agreement
between the empirical behavior and the theoretical prediction for the
full range of incoherence parameters. For Case 2 in panel (b), the
agreement is excellent for moderate and large values of the
incoherence parameters $\mu_i$ and $\nu_t$, and we observe some
discrepancies for very small values of incoherence (below $10^{-2}$).
This difference arises because in this extreme regime, certain
``higher-order'' terms---not part of the theoretical
prediction---become dominant. Overall, we see the entrywise
estimation error matches the theoretical scaling (hence is
statistically optimal) even when $\mu_i$ and/or $\nu_t$ are large, or
at least one of $\mu_i$ and $\nu_t$ is not vanishingly small. These
results support the practical validity of our algorithm and theory
over a wide range of problem instances.
\subsection{Scaling with rank}
In this section, we devised some experiments that expose the effect of
varying the matrix rank on our procedures. In order to do so, we set
$N = T = 500$, the noise level $\omega=1$, and consider the same
staggered adoption design with $k = 3$ groups with $N_{1} = T_{1} =
200$, $N_2 = T_2 = 200$ and $N_3 = T_3 = 100$.
\begin{figure}[h]
\begin{center}
\begin{tabular}{ccc}
\widgraph{0.45\textwidth}{figs/fig_rank_case1} &&
\widgraph{0.45\textwidth}{figs/fig_rank_case2} \\
(a) && (b)
\end{tabular}
\caption{ Red circles: Mean-squared error (MSE) of the estimate
$\widehat{M}_{i,t}$ of entry $(i,t) = (500, 500)$ returned by
Algorithm~\ref{alg:general} versus the rank $r$. Empirical MSE
was obtained as an average of $100$ Monte Carlo trials. Blue
lines: theoretically predicted scaling $\gamma_{i,t}^\star$ of
the MSE versus rank. (a) Plots for case $1$ where the theory
predicts scaling as $r$. (b) Plots for Case 2 where theory
predicts scaling as $r^{3/2}$.
\label{fig:rank}}
\end{center}
\end{figure}
As before, we consider estimating the missing entry
$(i,t)=(500,500)$. Under the sub-block
condition~\eqref{EqnSubmatrix-general}, the theoretical scaling
$\gamma_{i,t}^\star$ defined from equation~\eqref{eq:variance-defn} is
of the order
\begin{align*}
\ensuremath{\ensuremath{\gamma}^*}_{i, t} \asymp \frac{\ensuremath{\omega}^{2}}{N_1 T} \Vert \bm{U}_{i,
\cdot}^{\star} \Vert_2^2 + \frac{\ensuremath{\omega}^{2}}{N T_1} \Vert \bm{V}_{t,
\cdot}^{\star} \Vert_2^2.
\end{align*}
Although the rank $r$ does not appear explicitly, it arises implicitly
via the $\ell_2$-norm of two vectors $\bm{U}_{i,\cdot}^\star$ and
$\bm{V}_{t,\cdot}^\star \in \mathbb{R}^r$. Recall the definition of
$\bm{U}^\star (\cdot) \in \mathbb{R}^{N\times r}$ and $\bm{V}^\star
(\cdot) \in \mathbb{R}^{T \times r}$ from
equation~\eqref{eq:ensemble}, we consider the following two cases:
\begin{itemize}
\item {\bf Case 1:} We generate $\bm{M}^{\star}=\bm{U}^{\star} (1)
\bm{V}^{\star\top}(1)$ with $r$ varies between $1$ and $30$. This is
the ideal incoherent setting, and our theory predicts that the MSE
should scale as $\gamma_{i,t}^\star \asymp r$, so linearly in the
rank $r$.
\item {\bf Case 2:} We generate $\bm{M}^{\star}=\bm{U}^{\star}
(r^{1/4}) \bm{V}^{\star}(r^{1/4})^\top$ where $\bm{U}^\star (\cdot)
\in \mathbb{R}^{N\times r}$ and $\bm{V}^\star (\cdot) \in
\mathbb{R}^{T \times r}$ are defined in equations~\eqref{eq:U-mu}
and~\eqref{eq:V-nu}, with the rank $r$ varying between $1$ and
$30$. Under this setting, our theory predicts that the MSE should
scale as $\gamma_{i,t}^\star \asymp r^{3/2}$, so super-linearly in
the rank $r$.
\end{itemize}
\Cref{fig:rank} compares the mean-squared estimation
$\mathbb{E}[(\widehat{M}_{i,t}-M_{i,t}^\star)^2]$ and the theoretical
scaling $\gamma^\star_{i,t}$ for these two cases. Panel (a) shows
results for Case 1, where we expect the scaling with rank to be
roughly linear; here for the given dimensions $N = T = 500$, we see
good agreement between empirical and theory up to around rank $r
\approx 25$. On the other hand, panel (b) shows results for the more
challenging Case 2, where the theoretical effect of the rank is more
severe (growing as $r^{3/2}$). In this case, we see good agreement
between empirical and theoretical prediction up until around rank $r
\approx 20$. Given that the matrix dimension $N=T=500$ and $N_1 = T_1
= 200$, we can see that our algorithm and theory remain valid for over
a reasonably large range of ranks $r$.
\subsection{Computational costs}
Finally, we did some simple experiments to explore the computational
efficiency of our methods. Both of our procedures
(cf.~Algorithms~\ref{alg:4-blocks-Md}~and~\ref{alg:general}) involve
only the computation of partial SVDs, and solving least-squares
problems. Viewing the rank $r$ as an order one quantity, under the
four-block setting, the computational complexity of
Algorithm~\ref{alg:4-blocks-Md} is dominated by the cost of computing
the top-$r$ SVD of a $N$ by $T_1$ matrix and a $N_1$ by $T$ matrix;
the cost of doing so scales $O(N T_1+N_1 T)$, apart from log factors.
For the general staggered adoption setting with $k$ groups, the
computational complexity is
\begin{align*}
\sum_{i_0 = 1}^k \sum_{j_0 = k+1-i_0}^{k} \Bigg[ \bigg(
\sum_{i=1}^{k+1-j_0} N_i \bigg) \bigg( \sum_{j=1}^{j_0}T_j \bigg) +
\bigg( \sum_{i=1}^{i_0} N_i \bigg) \bigg( \sum_{j=1}^{k+1-i_0}T_j
\bigg) \Bigg].
\end{align*}
In rough terms, this upper bound scales as $O(k^2 NT)$, so that when
the number of groups $k$ is small, the computational cost is of the
same order as cost of reading the data (i.e., the $N T$ entries of the
matrix).
\begin{table}[h]
\centering
\begin{tabular}{c|c|c|c}
Dimension $N$ & Computation time (seconds) & Dimension $N$ &
Computation time (seconds) \\ \hline
512 & 0.10 & 4096 & 5.3 \\
1024 & 0.14 & 8192 & 22.3 \\
2048 & 0.94 & 16384 & 97.0 \\ \hline
\end{tabular}
\smallskip
\caption{The runtime of Algorithm~\ref{alg:general} vs.~the dimension
$N=T$. Reported runtimes are computed using a 2023 MacBook Pro with an Apple M2 Pro chip, and are averaged over 10
independent trials.
\label{table:runtime}}
\end{table}
To give some sense of the typical runtime of our algorithm, we
conducted some experiments on square matrices ($N = T$) with the side
length $N$ varying over the set $\{2^9,2^{10},\ldots,2^{15}\}$, and
with fixed rank $r = 3$. In all cases, we ran trials on a staggered
adoption design with $k = 3$ groups with $N_{1} = T_{1} = \lfloor 0.4N
\rfloor$ and $N_2 = T_2 = \lfloor 0.4N \rfloor$. Our code is
Python-based, making use of the standard \texttt{numpy} and
\texttt{scipy} routines for linear algebra.\footnote{One could likely
speed up our code substantially by incorporating various speed-ups,
e.g., from randomized numerical linear algebra; the current results
show runtimes without any such code optimization.} Based on our
discussion in the previous paragraph, we expect that the runtime of
Algorithm~\ref{alg:general} scales quadratically as $O(N^2)$; this is
confirmed by the simulation results shown in~\Cref{table:runtime}. We
note that this quadratic scaling is much faster than SDP-based
algorithms like convex relaxation considered in previous literature
(e.g.,~\cite{athey2021matrix}), which become computationally
burdensome even for moderate matrix dimensions (e.g., when $N >
3000$).
\section{Proofs}
\label{sec:Proof-outline}
In this section, we provide the proofs of certain results, including
our achievable guarantee in the four-block case
(\Cref{thm:distribution}) in~\Cref{subsec:proof-thm-distribution}; our
data-dependent procedure for computing confidence intervals
(\Cref{prop:CI}) in~\Cref{sec:proof-CI}; and our lower bounds
(\Cref{thm:CRLB}) in~\Cref{sec:proof-thm-crlb}.
In all cases, we defer the proofs of various intermediate results, of
a more technical nature, to the appendices. In addition, we defer the
proof of~\Cref{thm:distribution-general}
to~\Cref{sec:proof-general-distribution}.
\subsection{Proof of~\Cref{thm:distribution}}
\label{subsec:proof-thm-distribution}
Throughout this and other proofs, we use the shorthand notation
$\spicy \coloneqq \log(N + T)$.
\subsubsection{A master decomposition}
Our proof is based on a key decomposition, which we state as a
separate lemma in its own right. It involves the matrix
\begin{align}
\bm{Z} & \coloneqq \bm{U}_2^{\star} (\bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top} \bm{E}_{b} +
\bm{E}_{c} \bm{V}_{1}^{\star} (\bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1} \bm{V}_2^{\star\top} + \bm{E}_{c}
\bm{V}_{1}^{\star}(\bm{V}_{1}^{\star\top} \bm{V}_{1}^{\star})^{-1}
(\bm{\Sigma}^{\star})^{-1} (\bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top} \bm{E}_{b}, \label{eq:Z-defn}
\end{align}
as well as the error terms
\begin{align*}
\ensuremath{b}_1 (i,t) & \mydefn \frac{1}{\sqrt{N_1 T}}\bigg
(\frac{\ensuremath{\omega}^{3} }{\sigma_{r}^{\star2} N_1 } +
\frac{\ensuremath{\omega}^{2}}{\sigma_{r}^{\star} \sqrt{T_{1}}} \bigg ) \left \Vert
\bm{U}_{i,\cdot}^{\star} \right \Vert_{2} \sqrt{r + \spicy}, \\
\ensuremath{b}_2 (i,t) & \mydefn \frac{1}{\sqrt{N T_1} }\bigg
(\frac{\ensuremath{\omega}^{3} }{\sigma_{r}^{\star2} N_1 } +
\frac{\ensuremath{\omega}^{2}}{\sigma_{r}^{\star} \sqrt{N} } +
\frac{\ensuremath{\omega}^{2}}{\sigma_{r}^{\star} \sqrt{T_1} } \bigg ) \left \Vert
\bm{V}_{t,\cdot}^{\star} \right \Vert_{2} \sqrt{r + \spicy}, \\
\ensuremath{b}_3 (i,t) & \mydefn \bigg
(\ensuremath{\ensuremath{\rho}_{N,T}} + \ensuremath{\omega} \sqrt{\frac{r}{ N_{1} T } } + \ensuremath{\omega}
\sqrt{\frac{1}{\Nunit_1 \Time_1} \spicy } \bigg ) \left \Vert
\bm{U}_{i,\cdot}^{\star} \right \Vert _{2} \left \Vert
\bm{V}_{t,\cdot}^{\star} \right \Vert _{2}, \\
\ensuremath{b}_4 (i,t) & \mydefn \frac{1}{\sqrt{NT}} \frac{1}{\sqrt{N_1 T_1}} \, \bigg (\frac{\ensuremath{\omega}^{3} }{\sigma_{r}^{\star 2} \sqrt{N} }
+ \frac{\ensuremath{\omega}^{3}}{\sigma_{r}^{\star2} \sqrt{T_1} } + \frac{\ensuremath{\omega}^{4}}{\sigma_{r}^{\star3} N_1} \bigg )
\left(r + \spicy \right ).
\end{align*}
\begin{lemma}[\textsf{Master decomposition}]
\label{lem:master}
Suppose that the sub-block condition~\eqref{EqnSubmatrix} and the
noise condition~\eqref{eq:noise-condition-est} hold. Then we can write
\begin{subequations}
\begin{align}
\label{eq:Md-decomposition}
\widehat{\bm{M}}_{d} - \bm{M}_{d}^{\star} = \bm{Z} + \bm{\Delta},
\end{align}
for a perturbation matrix $\bm{\Delta}$ whose entries satisfy the
entrywise bound
\begin{align}
\big |\Delta_{i,t} \big | \leq C_4 \sum_{k=1}^4 \ensuremath{b}_k (i,t)
\qquad \mbox{for all entries $(i,t)$}
\end{align}
\end{subequations}
with probability at least $1 - O((N + T)^{-9})$.
\end{lemma}
\noindent See~\Cref{subsec:proof-lemma-master} for the proof. \\
\subsubsection{Main argument for~\Cref{thm:distribution}}
Armed with this result, we are now equipped to
prove~\Cref{thm:distribution}. We begin with the additive
decomposition $Z_{i,t} = \alpha_{i,t} + \beta_{i,t} + \lambda_{i,t}$,
where
\begin{subequations}
\label{eq:Zij-decom}
\begin{align}
\alpha_{i, t} & \mydefn (\bm{E}_{c})_{i,\cdot}
\bm{V}_{1}^{\star}( \bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1} \bm{V}_{t,\cdot}^{\star\top},
\\
\beta_{i, t} & \mydefn \bm{U}_{i,\cdot}^{\star} (
\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top}
( \bm{E}_{b})_{\cdot,t}, \quad \mbox{and} \\
\lambda_{i, t} & \mydefn ( \bm{E}_{c})_{i,\cdot} \bm{V}_{1}^{\star}(
\bm{V}_{1}^{\star\top} \bm{V}_{1}^{\star})^{-1}(
\bm{\Sigma}^{\star})^{-1}( \bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top}
(\bm{E}_{b})_{\cdot,t}.
\end{align}
\end{subequations}
Observe that
$\alpha_{i,t}$ and $\beta_{i,t}$ are independent zero-mean Gaussian
random variables with variances
\begin{subequations}
\label{eq:alpha-beta-var}
\begin{align}
\mathsf{var} \left (\alpha_{i,t} \right ) & = \frac{\ensuremath{\omega}^{2}}{NT} \big\Vert
\bm{U}_{i,\cdot}^{\star}( \bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star}){}^{-1} \bm{U}_{1}^{\star\top} \big
\Vert_{2}^{2} \overset{\text{(i)}}{\in} \bigg[ \frac{1}{\ensuremath{c_u}}
\frac{\ensuremath{\omega}^{2}}{N_{1} T} \Vert \bm{U}_{i,\cdot}^{\star}
\Vert_{2}^{2}, \frac{1}{\ensuremath{c_\ell}} \frac{\ensuremath{\omega}^{2}}{N_{1} T} \Vert
\bm{U}_{i,\cdot}^{\star} \Vert_{2}^{2} \bigg], \quad \mbox{and}
\\
\mathsf{var} \left (\beta_{i,t} \right ) & = \frac{\ensuremath{\omega}^{2}}{NT}
\big\Vert \bm{V}_{t,\cdot}^{\star}( \bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star}){}^{-1} \bm{V}_{1}^{\star\top}
\big\Vert_{2}^{2} \overset{\text{(ii)}}{\in} \bigg[
\frac{1}{\ensuremath{c_u}} \frac{\ensuremath{\omega}^{2}}{N \Time_1}\Vert
\bm{V}_{t,\cdot}^{\star}\Vert_{2}^{2}, \frac{1}{\ensuremath{c_\ell}}
\frac{\ensuremath{\omega}^{2}}{ N\Time_1}\Vert
\bm{V}_{t,\cdot}^{\star}\Vert_{2}^{2} \bigg],
\end{align}
\end{subequations}
where (i) and (ii) both follow from \eqref{EqnSubmatrix}. Then we can
check that
\begin{align}
\frac{b_1(i,t)}{\mathsf{var}^{1/2} (\alpha_{i,t} )} +
\frac{b_2(i,t)}{\mathsf{var}^{1/2} (\beta_{i,t} )} & \leq 3
\sqrt{\ensuremath{c_u}} ( \ensuremath{\ensuremath{\rho}_{N,T}}^2 + \ensuremath{\ensuremath{\rho}_{N,T}} ) \sqrt{r + \spicy}
\overset{\text{(i)}}{\leq} \frac{\delta}{4}, \nonumber
\\ \frac{b_3(i,t)}{\mathsf{var}^{1/2} (\alpha_{i,t} +
\beta_{i,t} )} & \leq \sqrt{\ensuremath{c_u}} \bigg ( \sqrt{r} (\ensuremath{\ensuremath{\rho}_N}
\nu_t + \ensuremath{\ensuremath{\rho}_T} \mu_i) + \min\left\{ \frac{\mu_{i}}{\sqrt{N_1}} ,
\frac{\nu_{t} }{\sqrt{T_1}} \right\} \sqrt{ r \spicy} \bigg )
\overset{\text{(ii)}}{\leq} \frac{\delta}{4}, \quad \text{and}
\nonumber \\ \frac{b_4(i,t)}{\mathsf{var}^{1/2} (\alpha_{i,t}
+ \beta_{i,t} )} & \leq 2 \sqrt{\ensuremath{c_u}} ( \ensuremath{\ensuremath{\rho}_{N,T}}^2 +
\ensuremath{\ensuremath{\rho}_{N,T}} ) \min \left\{ \frac{\ensuremath{\ensuremath{\rho}_T}}{\mu_i} , \frac{\ensuremath{\ensuremath{\rho}_N}}{\nu_t}
\right\} \left( \sqrt{r} + \frac{\spicy}{\sqrt{r}} \right )
\overset{\text{(iii)}}{\leq}
\frac{\delta}{4}. \label{eq:b-ratio}
\end{align}
Here step (i) follows from \eqref{eq:noise-condition-est}, step (ii)
utilizes the bound~\eqref{eq:incoherence-est}, whereas step (iii)
follows from the bounds~\eqref{eq:noise-condition-est}
and~\eqref{eq:signal-lb}. Conditioned on the random matrix
$\bm{E}_{b}$, we have
\begin{align*}
\lambda_{i,t} \sim \mathcal{N} \big(0, \frac{\ensuremath{\omega}^2}{NT} \big\Vert
\bm{V}_{1}^{\star}( \bm{V}_{1}^{\star\top} \bm{V}_{1}^{\star})^{-1}
(\bm{\Sigma}^{\star})^{-1} ( \bm{U}_{1}^{\star \top}
\bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star \top} (
\bm{E}_{b})_{\cdot,t} \big\Vert_{2}^{2} \big).
\end{align*}
Therefore, for any fixed $(i,t)$, with probability at least $1 - O((N
+ T)^{-10})$, we have
\begin{align*}
\left |\lambda_{i,t} \right | & \overset{\text{(i)}}{\leq} \frac{5 \ensuremath{\omega}}{\sqrt{NT}}
\big\Vert \bm{V}_{1}^{\star}( \bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1}( \bm{\Sigma}^{\star})^{-1}(
\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1}
\bm{U}_{1}^{\star\top}( \bm{E}_{b})_{\cdot,t} \big\Vert_{2}
\sqrt{\spicy} \nonumber \\
& \overset{\text{(ii)}}{\leq} \frac{5 C_g \ensuremath{\omega}^{2}}{NT} \big\Vert
\bm{V}_{1}^{\star}( \bm{V}_{1}^{\star\top} \bm{V}_{1}^{\star})^{-1}(
\bm{\Sigma}^{\star})^{-1}( \bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top} \big\Vert \sqrt{
\left (r + \spicy \right )\spicy}\nonumber \\
& \overset{\text{(iii)}}{\leq} \frac{5C_g}{\ensuremath{c_\ell}^2}
\frac{\ensuremath{\omega}^{2}}{\sigma_{r}^{\star}} \sqrt{\frac{ 1 }{ N T \Nunit_1 \Time_1} \left (r + \spicy \right )
\spicy}.
\end{align*}
Here the constant $5$ in step (i) comes from the fact that
$\mathbb{P}(X\geq 5\sqrt{\spicy} ) = O((N+T)^{-10})$ for
$X\sim\mathcal{N}(0,1)$, step (ii) uses the concentration bound
\[
\big\Vert \bm{\zeta} ( \bm{E}_{b})_{\cdot,t} \big\Vert_{2} \leq C_g \sigma \big\Vert
\bm{\zeta} \big\Vert \sqrt{
r + \spicy }, \quad \text{where} \quad \bm{\zeta} = \bm{V}_{1}^{\star}( \bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1}( \bm{\Sigma}^{\star})^{-1}(
\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1}
\bm{U}_{1}^{\star\top},
\]
which is a consequence of \Cref{lemma:gaussian-spectral} in \Cref{appendix:technical_lemmas}, whereas
step (iii) follows from \eqref{EqnSubmatrix}.
Taking the bound on $|\lambda_{i,t}|$ and \eqref{eq:signal-lb} together yields
\begin{align}
\frac{ | \lambda_{i,t} | }{ \mathsf{var}^{1/2} (\alpha_{i,t} + \beta_{i,t} ) } \leq \frac{5C_g \sqrt{\ensuremath{c_u}}}{\ensuremath{c_\ell}^2} \min \left\{
\frac{\ensuremath{\ensuremath{\rho}_T}}{\mu_i} , \frac{\ensuremath{\ensuremath{\rho}_N}}{\nu_t} \right\} \sqrt{\spicy + \spicy^2 / r} \leq \frac{\delta}{4}. \label{eq:lambda-ratio}
\end{align}
Let $G_{i,t} = \alpha_{i,t} + \beta_{i,t}$ and $\gamma_{i,t}^\star = \mathsf{var} (\alpha_{i,t}+\beta_{i,t})$, we can check that
\begin{align*}
\Biggr| \frac{1}{ (\ensuremath{\ensuremath{\gamma}^*}_{i,t})^{1/2} }
\big(\widehat{\bm{M}}_{d} - \bm{M}_{d}^{\star}\big)_{i, t} - G_{i, t}
\Biggr| \overset{\text{(i)}}{\leq} \sum_{j=1}^4
\frac{b_j(i,t)}{\mathsf{var}^{1/2} (\alpha_{i,t} + \beta_{i,t} )} +
\frac{ | \lambda_{i,t} | }{ \mathsf{var}^{1/2} (\alpha_{i,t} +
\beta_{i,t} ) } \overset{\text{(ii)}}{\leq} \delta
\end{align*}
where step (i) follows from~\Cref{lem:master} and step (ii) follows
from equations~\eqref{eq:b-ratio} and~\eqref{eq:lambda-ratio}.
\subsubsection{Proof of~\Cref{lem:master}}
\label{subsec:proof-lemma-master}
Our proof is based on decomposing the estimation error as
\begin{align}
\bm{M}_{d}- \bm{M}_{d}^{\star} & = \bm{U}_{2}( \bm{U}_{1}^{\top}
\bm{U}_{1})^{-1} \bm{U}_{1}^{\top} \bm{M}_{b}- \bm{U}_{2}^{\star}(
\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1}
\bm{U}_{1}^{\star\top} \bm{M}_{b}^{\star}\nonumber \\
& = \bm{U}_{2}( \bm{U}_{1}^{\top} \bm{U}_{1})^{-1} \bm{U}_{1}^{\top}(
\bm{M}_{b}- \bm{M}_{b}^{\star}) + \big[ \bm{U}_{2}( \bm{U}_{1}^{\top}
\bm{U}_{1})^{-1} \bm{U}_{1}^{\top}- \bm{U}_{2}^{\star}(
\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1}
\bm{U}_{1}^{\star\top} \big] \bm{M}_{b}^{\star} \nonumber \\
& = \bm{A} + \bm{B}, \label{eq:Md-error-decom}
\end{align}
where
\begin{align*}
\bm{A} \mydefn \bm{U}_{2}( \bm{U}_{1}^{\top} \bm{U}_{1})^{-1}
\bm{U}_{1}^{\top}( \bm{M}_{b}- \bm{M}_{b}^{\star}), \quad \mbox{and}
\quad \bm{B} \mydefn \big[ \bm{U}_{2}( \bm{U}_{1}^{\top}
\bm{U}_{1})^{-1} \bm{U}_{1}^{\top}- \bm{U}_{2}^{\star}(
\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1}
\bm{U}_{1}^{\star\top} \big] \bm{M}_{b}^{\star}.
\end{align*}
With these definitions in hand, we now state two auxiliary lemmas that
bound these two terms.
\begin{lemma}
\label{lemma:decom-A}
Under the conditions of \Cref{lem:master}, we have the decomposition
\begin{align}
\label{eq:decom-A}
\bm{A} = \bm{U}_{2}^{\star}( \bm{U}_{1}^{\star\top}
\bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top} \bm{E}_{b} +
\bm{E}_{c} \bm{V}_{1}^{\star}( \bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1}( \bm{\Sigma}^{\star})^{-1}(
\bm{U}_{1}^{\star\top} \bm{U}_{1}^{\star})^{-1} \bm{U}_{1}^{\star\top}
\bm{E}_{b} + \bm{\xi}_{ \bm{A}},
\end{align}
where $| ( \bm{\xi}_{ \bm{A}})_{i,t} | \leq C_A \sum_{k=1}^4 \ensuremath{b}_k (i,t) $ holds for each $(i,t)$ with probability at least $1-O((N + T)^{-9})$.
\end{lemma}
\noindent See~\Cref{subsec:proof-decom-A} for the proof.
\begin{lemma}
\label{lemma:decom-B}
Under the conditions of \Cref{lem:master}, there exists some universal constant $C_B>0$ such that we can decompose
\begin{align}
\label{eq:decom-B}
\bm{B} = \bm{E}_{c} \bm{V}_{1}^{\star}( \bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1} \bm{V}_{2}^{\star\top} + \bm{\xi}_{ \bm{B}},
\end{align}
where $| ( \bm{\xi}_{ \bm{B}})_{i,t} | \leq C_B \sum_{k=2}^3 \ensuremath{b}_k (i,t) $ holds for each $(i,t)$ with probability at least $1-O((N + T)^{-9})$.
\end{lemma}
\noindent See~\Cref{subsec:proof-decom-C} for the proof.
Combining equation~\eqref{eq:Md-error-decom} with these two lemmas
yields $ \bm{M}_{d} - \bm{M}_{d}^{\star} = \bm{Z} + \bm{\Delta}$,
where the random matrix $\bm{Z}$ was previously
defined~\eqref{eq:Z-defn}, and $\bm{\Delta} \mydefn \bm{\xi}_{\bm{A}}
+ \bm{\xi}_{\bm{B}}$. The residual matrix $\bm{\Delta}$ satisfies,
with probability exceeding $1-O((N+T)^{-9})$, that $ | \Delta_{ i , t}
| \leq C_4 \sum_{k=1}^4 b_k (i,t) $ holds for each $(i,t)$ with some
universal constant $C_4 = C_A + C_B$.
\subsection{Proof of~\Cref{prop:CI}}
\label{sec:proof-CI}
Naturally, a central step in the proof is to bound the
error incurred by using $\widehat{\gamma}_{i,t}$ as an estimate of
$\gamma_{i,t}^\star$. Let us summarize this auxiliary claim here:
with probability at least $1 - O((N+T)^{-10})$, there exists some constant $\ensuremath{C_{\gamma}}>0$ such that
\begin{align}
\label{eq:var-est-err}
|\widehat{\gamma}_{i,t} - \gamma_{i,t}^\star| & \leq \ensuremath{C_{\gamma}} \bigg( \tfrac{\delta}{\sqrt{\spicy}} + \sqrt{\tfrac{( N_{1} + T ) r
}{N_{1} T}} \bigg) \gamma_{i,t}^\star.
\end{align}
Here $\delta \in (0,1)$ is the tolerance parameter in the statement
of~\Cref{prop:CI}. We return to prove this result
in~\Cref{SecProofAuxClaim}.
\subsubsection{Main argument}
Taking the auxiliary claim~\eqref{eq:var-est-err} as given, let us
prove the coverage guarantee stated in~\Cref{prop:CI}. Define the
random variable
\begin{align*}
\ensuremath{V_{i,t}} \coloneqq (\widehat{M}_{i,t}-M_{i,t}^{\star}) /
\widehat{\gamma}_{i,t}^{1/2} - G_{i,t} ,
\end{align*}
where $G_{i,t}$ is the standard Gaussian variable defined
in~\Cref{thm:distribution}.
By applying~\Cref{thm:distribution}, we find that with probability at
least $1-O((N+T)^{-10})$,
\begin{align}
\label{EqnInterBound}
|\ensuremath{V_{i,t}}| & = \left | \tfrac{\widehat{M}_{i,t} - M_{i,t}^{\star}}
{(\gamma_{i,t}^\star)^{1/2}} - G_{i,t} \right | + \left |
\tfrac{\widehat{M}_{i,t} - M_{i,t}^{\star}}
{(\gamma_{i,t}^\star)^{1/2}} -
\tfrac{\widehat{M}_{i,t}-M_{i,t}^{\star}}
{(\widehat{\gamma}_{i,t})^{1/2}} \right | \; \leq \delta +
\left | \tfrac{\widehat{M}_{i,t} - M_{i,t}^{\star}}
{(\gamma_{i,t}^\star)^{1/2}} \right | \; \; \left |
1 -
\tfrac{(\gamma_{i,t}^\star)^{1/2}}{(\widehat{\gamma}_{i,t})^{1/2}}
\right |.
\end{align}
Next we observe that
\begin{align*}
\left | \tfrac{\widehat{M}_{i,t}-M_{i,t}^{\star}}
{(\gamma_{i,t}^\star)^{1/2}} \right | \leq \left |
\tfrac{\widehat{M}_{i,t}-M_{i,t}^{\star}}
{(\gamma_{i,t}^\star)^{1/2}} -G_{i,t} \right | + | G_{i,t}
| \overset{\text{(i)}}{\leq} \delta + 5\sqrt{\spicy}
\end{align*}
and
\begin{align*}
\left | 1 -
\tfrac{(\gamma_{i,t}^\star)^{1/2}}{(\widehat{\gamma}_{i,t})^{1/2}}
\right | = \tfrac{| \gamma_{i,t}^\star - \widehat{\gamma}_{i,t} |
}{(\widehat{\gamma}_{i,t})^{1/2} [ (\gamma_{i,t}^\star)^{1/2} +
(\widehat{\gamma}_{i,t})^{1/2} ] } \overset{\text{(ii)}}{\leq} 2 \ensuremath{C_{\gamma}} \tfrac{\delta}{\sqrt{\spicy}} + 2 \ensuremath{C_{\gamma}} \sqrt{\tfrac{(
N_{1} + T ) r }{N_{1} T}};
\end{align*}
where claim (i) follows from~\Cref{thm:distribution}, and step (ii)
follows from the auxiliary claim~\eqref{eq:var-est-err} and its consequence $\widehat{\gamma}_{i,t} \geq \gamma_{i,t}^\star / 2$,
provided that $\delta$ is sufficiently small and $\min\{N_1,T\}$ is sufficiently
large.
Combining these facts with our earlier bound~\eqref{EqnInterBound} and
defining $C_v = 12 \ensuremath{C_{\gamma}}$, we find that
\begin{align}
\label{eq:epsilon-bound}
|\ensuremath{V_{i,t}}| & \leq \delta + (\delta + 5\sqrt{\spicy} ) \left( 2 \ensuremath{C_{\gamma}} \tfrac{\delta}{\sqrt{\spicy}} + 2 \ensuremath{C_{\gamma}} \sqrt{\tfrac{( N_{1} + T
) r }{N_{1} T}} \right) \; \overset{\text{(iii)}}{\leq} C_v
\delta,
\end{align}
where step (iii) holds provided that $\min\{ N_1, T \} \geq
\delta^{-2} r \spicy$.
We now use this high probability bound on $\ensuremath{V_{i,t}}$ to establish our
coverage guarantee. For any $x \in \mathbb{R}$, from the definition of
$\ensuremath{V_{i,t}}$, we have
\begin{subequations}
\begin{align}
\label{EqnStar}
\mathbb{P} \bigg( \tfrac{\widehat{M}_{i,t} -
M_{i,t}^{\star}}{(\widehat{\gamma}_{i,t})^{1/2}} \leq x \bigg) &
\stackrel{(\star)}{=} \mathbb{P} \big( G_{i,t} + \ensuremath{V_{i,t}} \leq x, |
\ensuremath{V_{i,t}} | \leq C_v \delta \big) + \mathbb{P} \big( G_{i,t} + \ensuremath{V_{i,t}}
\leq x, | \ensuremath{V_{i,t}} | > C_v \delta \big) \\
& \leq \mathbb{P} \big( G_{i,t} \leq x + C_v \delta \big) + \mathbb{P}
\left( | \ensuremath{V_{i,t}} | > C_v \delta \right) \notag \\
& \overset{\text{(a)}}{=} \Phi \left ( x + C_v \delta \right )+ O
\left((N+T)^{-10}\right) \notag \\
\label{EqnBoundOne}
& \overset{\text{(b)}}{\leq} \Phi(x) + O\big( \delta + (N+T)^{-10} \big)
\end{align}
where step (a) follows from the bound~\eqref{eq:epsilon-bound}; and
step (b) follows since the normal CDF $\Phi$ is a
$1/\sqrt{2\pi}$-Lipschitz function.
Similarly, we have the lower bound
\begin{align}
\mathbb{P} \bigg( \tfrac{\widehat{M}_{i,t} -
M_{i,t}^{\star}}{(\widehat{\gamma}_{i,t})^{1/2}} \leq x \bigg) &
\overset{\text{(c)}}{\geq} \mathbb{P} \big( G_{i,t} + \ensuremath{V_{i,t}} \leq x,
| \ensuremath{V_{i,t}} | \leq C_v \delta \big) \notag
\\
& \geq \mathbb{P} \big( G_{i,t} \leq x - C_v \delta \big) \notag \\
& = \Phi
\left( x - C_v \delta \right) \notag \\
\label{EqnBoundTwo}
& \overset{\text{(d)}}{\geq} \Phi(x) -O (\delta),
\end{align}
\end{subequations}
where step (c) follows from equality $(\star)$ in
equation~\eqref{EqnStar}; and step (d) follows from the Lipschitz
property of $\Phi$.
Finally, by taking $x=\pm \Phi^{-1} (1-\alpha/2)$ and combining the
upper and lower bounds~\eqref{EqnBoundOne} and~\eqref{EqnBoundTwo},
we find that
\begin{align*}
\mathbb{P} \big( \mathsf{CI}_{i,t}^{1-\alpha} \ni M_{i,t}^\star
\big) = \mathbb{P} \bigg(
\tfrac{\widehat{M}_{i,t}-M_{i,t}^{\star}}{(\widehat{\gamma}_{i,t})^{1/2}}
\in \big[ \pm \Phi^{-1} (1-\alpha/2)\big] \bigg) = 1 - \alpha + O
\big( \delta + (N+T)^{-10} \big),
\end{align*}
as claimed in~\Cref{prop:CI}.
\subsubsection{Proof of the auxiliary claim~\eqref{eq:var-est-err}}
\label{SecProofAuxClaim}
To simplify presentation, we introduce the shorthand $\sigma =
\omega/\sqrt{NT}$, along with the associated estimate
$\widehat{\sigma}^2 = \widehat{\omega}^2 / NT$. Using this notation,
we write
\begin{align*}
& \widehat{\gamma}_{i,t} - \gamma_{i,t}^{\star} = \underbrace{(
\widehat{\sigma}^2 - \sigma^{2} )
\bm{U}_{i,\cdot}^{\star}(\bm{U}_{1}^{\star\top}\bm{U}_{1}^{\star})^{-1}
\bm{U}_{i,\cdot}^{\star\top} + (\widehat{\sigma}^{2} - \sigma^2 )
\bm{V}_{t,\cdot}^{\star}(\bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1}\bm{V}_{t,\cdot}^{\top}}_{\eqqcolon
\beta_1} \\
& \quad + \underbrace{\widehat{\sigma}^{2} \big( \bm{U}_{i,\cdot}
(\bm{U}_{1}^{\top}\bm{U}_{1})^{-1}\bm{U}_{i,\cdot}^{\top} -
\bm{U}_{i,\cdot}^{\star}(\bm{U}_{1}^{\star\top}\bm{U}_{1}^{\star})^{-1}
\bm{U}_{i,\cdot}^{\star\top} \big) + \widehat{\sigma}^{2} \big(
\bm{V}_{t,\cdot}^{\star}(\bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1}\bm{V}_{t,\cdot}^{\top} -
\bm{V}_{t,\cdot}^{\star}(\bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1}\bm{V}_{t,\cdot}^{\star\top} \big)
}_{\eqqcolon\beta_2}
\end{align*}
Next we bound each of the terms $\beta_1$ and $\beta_2$ in turn.
\paragraph{Bounding $\beta_1$:}
In order to bound $\beta_1$, we need to control the estimation error
associated with $\widehat{\sigma}^2$. The following result provides
such a bound:
\begin{lemma}
\label{LemSigBound}
Defining the constant $\ensuremath{C_{\sigma}} = 2 \sqrt{2\ensuremath{C_{\mathsf{b}}} + 1} + \ensuremath{C_{\mathsf{upper}}}^2 + 2\ensuremath{C_{\mathsf{b}}}$, we
have
\begin{align}
\label{eq:noise-est}
\vert \widehat{\sigma}^{2}-\sigma^{2}\vert \leq \alpha_1 + \alpha_2 +
\alpha_3 \leq \ensuremath{C_{\sigma}} \sigma^{2} \sqrt{\tfrac{( N_{1} + T ) r }{N_{1}
T}},
\end{align}
with probability at least $1 - O((N + T)^{-10})$.
\end{lemma}
\noindent See~\Cref{SecProofLemSigBound} for the proof. \\
We now apply the relation~\eqref{eq:noise-est} to bound $\beta_1$,
thereby obtaining
\begin{subequations}
\begin{align}
| \beta_1 | & \leq \ensuremath{C_{\sigma}} \sigma^{2} \sqrt{\tfrac{( N_{1} + T ) r
}{N_{1} T}} \big(
\bm{U}_{i,\cdot}^{\star}(\bm{U}_{1}^{\star\top}\bm{U}_{1}^{\star})^{-1}
\bm{U}_{i,\cdot}^{\star\top} +
\bm{V}_{t,\cdot}^{\star}(\bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1}\bm{V}_{t,\cdot}^{\top} \big) \nonumber \\
\label{EqnBetaOne}
& = \ensuremath{C_{\sigma}} \sqrt{\tfrac{( N_{1} + T ) r }{N_{1} T}} \gamma_{i,t}^\star.
\end{align}
\paragraph{Bounding $\beta_2$:}
In order to bound the term $\beta_2$, we need the following lemma
characterizing the error of the plug-in estimate.
\begin{lemma}
\label{lemma:CI-2}
Under the conditions of~\Cref{prop:CI}, with probability at least $1 -
O((N + T)^{-10})$,
\begin{align*}
\sigma^2 \big | \bm{U}_{i,\cdot}
(\bm{U}_{1}^{\top}\bm{U}_{1})^{-1}\bm{U}_{i,\cdot}^{\top} -
\bm{U}_{i,\cdot}^{\star}(\bm{U}_{1}^{\star\top}\bm{U}_{1}^{\star})^{-1}
\bm{U}_{i,\cdot}^{\star\top} \big | + \sigma^2 \big | \bm{V}_{t,\cdot}
(\bm{V}_{1}^{\top} \bm{V}_{1})^{-1}\bm{V}_{t,\cdot}^{\top} -
\bm{V}_{t,\cdot}^{\star}(\bm{V}_{1}^{\star\top}
\bm{V}_{1}^{\star})^{-1}\bm{V}_{t,\cdot}^{\star\top} \big | & \leq
\frac{\delta \, \gamma_{i,t}^\star }{\sqrt{\spicy}} .
\end{align*}
\end{lemma}
\noindent See~\Cref{sec:proof-lemma-CI-2} for the proof.
When $\min\{N_1,T\} \geq r$, a direct consequence of \eqref{eq:noise-est} is that $\widehat{\sigma}^2 \leq (\sqrt{2}\ensuremath{C_{\sigma}}+1) \sigma^2$. Combining this bound with \Cref{lemma:CI-2} yields
\begin{align}
\label{EqnBetaTwo}
| \beta_2 | & \leq ( \sqrt{2} \ensuremath{C_{\sigma}} + 1) \tfrac{\delta}{\sqrt{\spicy}}
\gamma_{i,t}^\star,
\end{align}
\end{subequations}
\vspace*{0.1in}
Collecting together the two bounds~\eqref{EqnBetaOne}
and~\eqref{EqnBetaTwo} on $\beta_1$ and $\beta_2$, respectively,
yields the upper bound
\begin{align*}
|\widehat{\gamma}_{i,t} - \gamma_{i,t}^\star| & \leq |\beta_1| +
|\beta_2| \leq \bigg( ( \sqrt{2} \ensuremath{C_{\sigma}} + 1 ) \tfrac{\delta}{\sqrt{\spicy}} +
\ensuremath{C_{\sigma}} \sqrt{\tfrac{( N_{1} + T ) r }{N_{1} T}} \bigg)
\gamma_{i,t}^\star.
\end{align*}
Let $\ensuremath{C_{\gamma}} = \sqrt{2} \ensuremath{C_{\sigma}} +1$ to achieve the claimed bound~\eqref{eq:var-est-err}.
\subsubsection{Proof of~\Cref{LemSigBound}}
\label{SecProofLemSigBound}
We first decompose $\widehat{\sigma}^{2}$ into a sum of three terms as
\begin{align*}
\widehat{\sigma}^{2} & = \frac{1}{N_{1} T} \big\Vert
\bm{M}_{\mathsf{upper}}^{\star} + \bm{E}_{\mathsf{upper}} -
\widehat{\bm{M}}_{\mathsf{upper}} \big\Vert_{\mathrm{F}}^{2} \\ &=
\frac{1}{N_{1}T} \Vert \bm{E}_{\mathsf{upper}} \Vert_{\mathrm{F}}^{2}
+ \frac{1}{N_{1} T} \big\Vert \bm{M}_{\mathsf{upper}}^{\star} -
\widehat{\bm{M}}_{\mathsf{upper}} \big\Vert_{\mathrm{F}}^{2} +
\frac{2}{N_{1} T} \big\langle \bm{M}_{\mathsf{upper}}^{\star} -
\widehat{\bm{M}}_{\mathsf{upper}} , \bm{E}_{\mathsf{upper}}
\big\rangle.
\end{align*}
As our analysis will show, the first quantity concentrates around
$\sigma^2$. Thus, it is natural to bound the estimation error
associated with $\widehat{\sigma}^2$ in the following way:
\begin{align*}
\vert \widehat{\sigma}^{2} - \sigma^{2} \vert & \leq
\underbrace{\Big\vert \frac{1}{N_{1}T} \Vert\bm{E}_{\mathsf{upper}}
\Vert_{\mathrm{F}}^{2} - \sigma^{2} \Big\vert }_{\eqqcolon
\alpha_1} + \underbrace{\frac{1}{N_{1}T} \big\Vert
\bm{M}_{\mathsf{upper}}^{\star} - \widehat{\bm{M}}_{\mathsf{upper}}
\big\Vert_{\mathrm{F}}^{2} }_{\eqqcolon \alpha_2} +
\underbrace{\frac{2}{N_1 T} \big\Vert \bm{M}_{\mathsf{upper}}^{\star}
- \bm{M}_{\mathsf{upper}} \big\Vert_{\mathrm{F}} \Vert
\bm{E}_{\mathsf{upper}} \Vert_{\mathrm{F}}}_{\eqqcolon \alpha_3}.
\end{align*}
By applying Bernstein's inequality \citep[Theorem
2.8.1]{vershynin2016high}, we can conclude that for some constant
$\ensuremath{C_{\mathsf{b}}}>0$, as long as $N_1 T \geq \spicy$, we find that
\begin{subequations}
\label{eq:E-upper-fro}
\begin{align}
\alpha_1 = \frac{1}{N_1 T} \Big | \sum_{i=1}^{N_1} \sum_{j=1}^T (
E_{i,j}^2 - \sigma^2 ) \Big | \leq \ensuremath{C_{\mathsf{b}}} \sigma^{2}\sqrt{\frac{ \spicy
}{N_1 T}} + \ensuremath{C_{\mathsf{b}}} \sigma^2 \frac{\spicy}{N_1 T} \leq 2 \ensuremath{C_{\mathsf{b}}}
\sigma^{2}\sqrt{\frac{ \spicy }{N_1 T}}
\end{align}
with probability at least $1-O((N+T)^{-10})$. As a consequence, we
have
\begin{align}
\Vert \bm{E}_{\mathsf{upper}} \Vert_{\mathrm{F}}^2 \leq N_1 T \sigma^2
+ N_1 T \alpha_1 \leq (2 \ensuremath{C_{\mathsf{b}}} + 1 ) N_1 T \sigma^2,
\end{align}
\end{subequations}
provided that $N_1 T \geq \spicy$.
\noindent In order to bound the terms $\alpha_2$ and $\alpha_3$, we
make use of the following lemma.
\begin{lemma}
\label{lemma:CI-1}
Under the conditions of~\Cref{prop:CI}, we have
\begin{align*}
\big \Vert \bm{M}_{\mathsf{upper}}^{\star} -
\widehat{\bm{M}}_{\mathsf{upper}}\big\Vert_{\mathrm{F}} & \leq \ensuremath{C_{\mathsf{upper}}}
\sigma \sqrt{ ( T + N_{1} ) r}
\end{align*}
with probability at least $1 - O((N + T)^{-10})$.
\end{lemma}
\noindent See~\Cref{sec:proof-lemma-CI-1} for the proof.\\
Combining~\Cref{lemma:CI-1} with the bounds~\eqref{eq:E-upper-fro}
yields
\begin{align}
\alpha_2 \leq \ensuremath{C_{\mathsf{upper}}}^2 \frac{\sigma^{2} \left( T + N_{1} \right)
r}{N_{1} T} \qquad \text{and} \qquad \alpha_3 \leq 2 \sqrt{2\ensuremath{C_{\mathsf{b}}} +
1} \ensuremath{C_{\mathsf{upper}}} \sigma^{2} \sqrt{\frac{ ( N_{1} + T ) r }{N_{1} T}}.
\end{align}
Putting together the upper bounds for $\alpha_1$, $\alpha_2$ and
$\alpha_3$ yields the claimed bound~\eqref{eq:noise-est}.
\subsection{Proof of lower bounds}
\label{sec:proof-thm-crlb}
In this section, we prove the Cram\'{e}r--Rao and local minimax bounds
stated in the main text. While our formal statement (\Cref{thm:CRLB})
only covers the four-block case, here we actually provide a proof for
the multi-block case, so that we match our general set of achievable
results (\Cref{thm:distribution-general}).
Let us recall the set-up for our genie-aided problem associated with
estimating entry $\bm{M}_{i,t}^{\star}$ of the matrix. The only
unknown quantities are the vectors $\xstar \coloneqq
\bm{X}_{i,\cdot}^{\star}$ and $\ystar \coloneqq
\bm{Y}_{t,\cdot}^{\star}$, and our goal is to estimate the inner
product $\bm{M}_{i,t}^{\star} = \inprod{\xstar}{\ystar}$ based on
linear measurements of the form
\begin{align*}
\big\{ \bm{M}_{i, s} = \inprod{\xstar}{\bm{Y}_{s, \cdot}^{\star}} +
E_{i, s} \big\}_{s=1}^{\bar{T}_1} \qquad \text{and} \quad \big\{
\bm{M}_{k, t} = \inprod{\bm{X}_{k, \cdot}^{\star}}{\ystar} + E_{k, t}
\big\}_{k=1}^{\bar{\Nunit}_1}.
\end{align*}
\paragraph{Cram\'{e}r--Rao lower bound:}
We first compute the Cram\'{e}r--Rao lower bound (CRLB) for the
problem. Using $\bftheta = (\xstar, \ystar)$ to denote the unknown
parameters, our goal is to estimate the functional $\Psi(\bftheta)
\coloneqq \inprod{\xstar}{\ystar}$. The CRLB is given by $\nabla
\Psi(\bftheta)^T \FishMat^{-1}(\bftheta) \nabla \Psi(\bftheta)$, where
$\FishMat(\bftheta)$ is the Fisher information matrix, and the
gradient takes the form $\nabla \Psi(\bftheta) = \begin{bmatrix}
\ystar & \xstar
\end{bmatrix}^T$.
Given that our measurements are linear with Gaussian observation noise
with variance $\ensuremath{\omega}^2/(N T)$, the Fisher information matrix takes
the form
\begin{align*}
\FishMat(\bftheta) = \frac{NT}{\ensuremath{\omega}^{2}} \mbox{blkdiag} \Big(
\sum_{s = 1}^{\bar{T}_{1}} \bm{Y}_{s,\cdot}^{\star \top}
\bm{Y}_{s,\cdot}^{\star}, \quad \sum_{k = 1}^{\bar{N}_{1}}
\bm{X}_{k,\cdot}^{\star \top} \bm{X}_{k, \cdot}^\star \Big).
\end{align*}
Putting together the pieces, we find that the CRLB is given by
\begin{align*}
\mathsf{CRLB} \left (M_{i,t}^{\star} \right) & =
\bm{Y}_{t,\cdot}^{\star} \big ( \frac{NT}{\ensuremath{\omega}^{2}} \sum_{s =
1}^{\bar{T}_{1}} \bm{Y}_{s,\cdot}^{\star \top}
\bm{Y}_{s,\cdot}^{\star} \big )^{-1} \bm{Y}_{t,\cdot}^{\star \top} +
\bm{X}_{i,\cdot}^{\star} \big ( \frac{NT}{\ensuremath{\omega}^{2}} \sum_{k =
1}^{\bar{N}_{1}} \bm{X}_{k,\cdot}^{\star \top}
\bm{X}_{k,\cdot}^{\star} \big )^{-1} \bm{X}_{i,\cdot}^{\star \top} \\
& \overset{\text{(i)}}{=} \frac{\ensuremath{\omega}^{2}}{NT} \Big \{
\bm{V}_{t,\cdot}^{\star} \big (\sum_{s = 1}^{\bar{T}_{1}}
\bm{V}_{s,\cdot}^{\star \top} \bm{V}_{s,\cdot}^{\star} \big )^{-1}
\bm{V}_{t,\cdot}^{\star \top} + \bm{U}_{k,\cdot}^{\star} \big (\sum_{k
= 1}^{\bar{N}_{1}} \bm{U}_{k,\cdot}^{\star \top}
\bm{U}_{k,\cdot}^{\star} \big )^{-1} \bm{U}_{i,\cdot}^{\star \top}
\Big \} \\
& = \underbrace{\frac{\ensuremath{\omega}^2}{NT} \Big \{ \bm{V}_{t,\cdot}^{\star}(
\bm{V}_{1}^{\star \top} \bm{V}_{1}^{\star})^{-1}
\bm{V}_{t,\cdot}^{\star \top} + \bm{U}_{i,\cdot}^{\star}(
\bm{U}_{1}^{\star \top} \bm{U}_{1}^{\star})^{-1}
\bm{U}_{i,\cdot}^{\star \top} \Big \}}_{ \equiv \gamma_{i,t}^\star}
\end{align*}
Here step (i) uses the fact that $\bm{X}^\star = \bm{U}^\star
(\bm{\Sigma}^\star)^{1/2}$ and $\bm{Y}^\star = \bm{V}^\star
(\bm{\Sigma}^\star)^{1/2}$, along with some algebra.
\paragraph{Computation of the local minimax bound:}
We now turn to computation of the local minimax lower bound, where we
make use of a Bayesian form of the CRLB. We begin by observing that
for any $\varepsilon > 0$, the local minimax risk for estimating
$\bm{M}^\star_{i,t}$ can be lower bounded as
\begin{align*}
R^\star(\varepsilon) & \geq \inf_{\widehat{M}_{i,t}}
\mathop{\mathbb{E}}_{ (\bm{x} , \bm{y}) \sim \pi } \mathbb{E} \left[
\big (\widehat{M}_{i,t} - \inprod{\bm{x}}{\bm{y}} \big)^2 \right],
\end{align*}
where $\pi$ is any prior supported on the Cartesian product space.
$\ensuremath{\mathcal{B}}_\infty (\bm{x}^\star , \varepsilon ) \times \ensuremath{\mathcal{B}}_\infty
(\bm{y}^\star , \varepsilon )$.
Following a standard avenue for computing Bayesian
CRLBs~\cite{polyanskiy2023information}, we first define the squared
cosine density function $g(u) = \cos^2 (\pi u / 2)$ for $u \in [-1 ,
1]$. Using this building block, we then define the product prior
distribution $\pi(\theta_1, \ldots, \theta_{2r}) = \prod_{j=1}^{2 r}
\pi_j(\theta_j)$ with marginals
\begin{align*}
\pi_j(\theta) & = \begin{cases} \frac{1}{\varepsilon} g \left(
\frac{\theta - x_{j}^\star }{\varepsilon} \right) & \mbox{for
$\theta_j \in [ x_j^\star - \varepsilon, x_j^\star +
\varepsilon]$ and $j = 1, \ldots, r$, and} \\
\frac{1}{\varepsilon} g \left( \frac{\theta -
y_{j-r}^\star}{\varepsilon} \right) & \mbox{for $\theta_{j} \in
[y_{j-r}^\star - \varepsilon, y_{j-r}^\star + \varepsilon ]$ and
$j = r+1, \ldots, 2r$.}
\end{cases}
\end{align*}
With this choice of prior, we can compute the lower bound
\begin{align*}
R^\star(\varepsilon) & \stackrel{\text{(i)}}{\geq}
\bm{Y}_{t,\cdot}^{\star} \big ( \frac{NT}{\ensuremath{\omega}^{2}} \sum_{s =
1}^{\bar{T}_{1}} \bm{Y}_{s,\cdot}^{\star \top}
\bm{Y}_{s,\cdot}^{\star} + \frac{\pi^2 }{\varepsilon^2} \bm{I}_r
\big )^{-1} \bm{Y}_{t,\cdot}^{\star \top} + \bm{X}_{i,\cdot}^{\star}
\big ( \frac{NT}{\ensuremath{\omega}^{2}} \sum_{k = 1}^{\bar{N}_{1}}
\bm{X}_{k,\cdot}^{\star \top} \bm{X}_{k,\cdot}^{\star} + \frac{\pi^2
}{\varepsilon^2} \bm{I}_r \big )^{-1} \bm{X}_{i,\cdot}^{\star \top}
\\
& \stackrel{\text{(ii)}}{=} \underbrace{\frac{\ensuremath{\omega}^{2}}{NT}
\bm{V}_{t,\cdot}^{\star} \big ( \bm{V}_{1}^{\star \top}
\bm{V}_{1}^{\star} + \frac{\pi^2 \ensuremath{\omega}^{2} }{\varepsilon^2 NT}
(\bm{\Sigma}^\star)^{-1} \big )^{-1} \bm{V}_{t,\cdot}^{\star \top}
+ \frac{\ensuremath{\omega}^{2}}{NT} \bm{U}_{i,\cdot}^{\star} \big (
\bm{U}_{1}^{\star \top} \bm{U}_{1}^{\star} + \frac{\pi^2
\ensuremath{\omega}^{2} }{\varepsilon^2 NT} (\bm{\Sigma}^\star)^{-1} \big
)^{-1} \bm{U}_{i,\cdot}^{\star \top}}_{ \eqqcolon
\widetilde{\gamma}_{i,t}(\epsilon)}
\end{align*}
where step (i) makes use of some standard computations with squared
cosine priors (e.g., \S 29.2 in the
book~\cite{polyanskiy2023information}); and step (ii) follows from the
spectral decomposition of $\bm{M}^\star$, along with some algebra.
Thus far, we proven that $R^\star(\varepsilon) \geq
\widetilde{\gamma}_{i,t} (\varepsilon) $. It remains to show that
with the choice $\varepsilon = \varepsilon_{N,T} \equiv
\sqrt{\sigma_r^\star} \max \{1/\sqrt{N}, 1/\sqrt{T} \}$ we have
\begin{align}
\label{EqnAux}
\Delta_{i,t}^\star \coloneqq \gamma_{i,t}^\star -
\widetilde{\gamma}_{i,t} (\varepsilon) \leq \gamma_{i,t}^\star \;
\frac{\ensuremath{c_u} \pi^2}{c_\ell} \ensuremath{\ensuremath{\rho}_{N,T}}^2,
\end{align}
from which it will follow that $R^\star(\varepsilon) \geq
\widetilde{\gamma}_{i,t} (\varepsilon) = \gamma^\star_{i,t} -
\Delta^\star_{i,t} \geq \big(1 - \frac{\ensuremath{c_u} \pi^2}{c_\ell}
\ensuremath{\ensuremath{\rho}_{N,T}}^2 \big) \gamma^\star_{i,t}$, as claimed.
From the definitions of $\widetilde{\gamma}_{i,t}(\varepsilon)$ and
$\gamma_{i,t}^\star$, we have
\begin{align*}
\Delta^*_{i,t} = \gamma_{i,t}^\star - \widetilde{\gamma}_{i,t} (\varepsilon) & =
\frac{\pi^2 \ensuremath{\omega}^{4} }{\varepsilon^2 (NT)^2}
\bm{V}_{t,\cdot}^{\star}( \bm{V}_{1}^{\star \top}
\bm{V}_{1}^{\star})^{-1} (\bm{\Sigma}^\star)^{-1} \big (
\bm{V}_{1}^{\star \top} \bm{V}_{1}^{\star} + \frac{\pi^2 \ensuremath{\omega}^{2}
}{\varepsilon^2 NT} (\bm{\Sigma}^\star)^{-1} \big )^{-1}
\bm{V}_{t,\cdot}^{\star \top} \\
& \qquad + \frac{\pi^2 \ensuremath{\omega}^{4} }{\varepsilon^2 (NT)^2}
\bm{U}_{i,\cdot}^{\star}( \bm{U}_{1}^{\star \top}
\bm{U}_{1}^{\star})^{-1} (\bm{\Sigma}^\star)^{-1} \big (
\bm{U}_{1}^{\star \top} \bm{U}_{1}^{\star} + \frac{\pi^2 \ensuremath{\omega}^{2}
}{\varepsilon^2 NT} (\bm{\Sigma}^\star)^{-1} \big)^{-1}
\bm{U}_{i,\cdot}^{\star \top} \\
& \leq \frac{\pi^2 \ensuremath{\omega}^{4} }{\varepsilon^2 (NT)^2} \Big(
\Vert (
\bm{V}_{1}^{\star \top} \bm{V}_{1}^{\star})^{-1}
\Vert^2 \Vert
(\bm{\Sigma}^\star)^{-1} \Vert \Vert \bm{V}_{t,\cdot}^{\star}
\Vert_2^2 + \Vert ( \bm{U}_{1}^{\star \top} \bm{U}_{1}^{\star})^{-1}
\Vert \Vert^2 (\bm{\Sigma}^\star)^{-1} \Vert \Vert
\bm{U}_{i,\cdot}^{\star} \Vert_2^2 \Big).
\end{align*}
Now observe that $\Vert (\bm{\Sigma}^\star)^{-1} \Vert \leq
\frac{1}{\sigma_r^\star}$, and that the sub-block condition
condition~\eqref{EqnSubmatrix-general} ensures that
\begin{align*}
\Vert (\bm{V}_{1}^{\star \top} \bm{V}_{1}^{\star})^{-1} \Vert^2 \leq
\tfrac{1}{\ensuremath{c_\ell}^2} \tfrac{T^2}{T_1^2} \quad \mbox{and} \quad \Vert
(\bm{U}_{1}^{\star \top} \bm{U}_{1}^{\star})^{-1} \Vert^2 \leq
\tfrac{1}{\ensuremath{c_\ell}^2} \tfrac{N^2}{N_1^2}.
\end{align*}
Combining these facts, we find that
\begin{subequations}
\begin{align}
\label{EqnSierraOne}
\Delta^*_{i,t} & \leq \frac{\pi^2 \ensuremath{\omega}^{4} }{ \ensuremath{c_\ell}^2 \varepsilon^2
\sigma_r^\star (NT)^2} \Big( \frac{T^2}{T_1^2} \Vert
\bm{V}_{t,\cdot}^{\star} \Vert_2^2 + \frac{N^2}{N_1^2} \Vert
\bm{U}_{i,\cdot}^{\star} \Vert_2^2 \Big).
\end{align}
On the other hand, using the $\ensuremath{c_u}$-bound in the sub-block
condition~\eqref{EqnSubmatrix-general}, we have
\begin{align}
\label{EqnSierraTwo}
\frac{1}{\ensuremath{c_u}} \ensuremath{\omega}^{2} \frac{1}{N_{1} T} \Vert
\bm{U}_{i,\cdot}^{\star} \Vert_{2}^{2} + \frac{1}{\ensuremath{c_u}}
\ensuremath{\omega}^{2} \frac{1}{N T_{1}} \Vert \bm{V}_{t,\cdot}^{\star}
\Vert_{2}^{2} \leq \gamma_{i,t}^\star.
\end{align}
\end{subequations}
Combining inequalities~\eqref{EqnSierraOne} and~\eqref{EqnSierraTwo}
yields the bound
\begin{align*}
\Delta^*_{i,t} & \leq \frac{\ensuremath{c_u} \pi^2 \ensuremath{\omega}^{2}}{ \ensuremath{c_\ell}^2
\varepsilon^2 \sigma_r^\star} \max \Big\{ \frac{1}{N T_1},
\frac{1}{N_1 T} \Big\} \gamma_{i,t}^\star \leq \frac{\ensuremath{c_u} \pi^2 \ensuremath{\omega}^{2} }{ \ensuremath{c_\ell}^2 \sigma_r^{\star, 2}}
\frac{1}{\min \{T_1, N_1 \}}
\underbrace{\frac{\sigma_r^\star}{\epsilon^2} \max \{ \frac{1}{N},
\frac{1}{T} \}}_{=1} \, \gamma_{i,t}^\star \overset{\text{(i)}}{\leq} \frac{\ensuremath{c_u} \pi^2 }{ \ensuremath{c_\ell}^2 }
\ensuremath{\ensuremath{\rho}_{N,T}}^2 \gamma_{i,t}^\star,
\end{align*}
where step (i) follows from our choice of $\varepsilon$, and the
definition of $\ensuremath{\ensuremath{\rho}_{N,T}}$. This establishes the auxiliary
claim~\eqref{EqnAux}, and thus the claimed lower bound.
\section{Discussion}
In this paper, we proposed and analyzed a simple algorithm for
estimation and inference for panel data with missingness induced by
staggered adoption design. It is appealing from a computational point
of view, since it is non-iterative, requiring only elementary matrix
operations and singular value decomposition. At the same time, we
demonstrated that its statistical properties are also attractive, in
that it achieves estimation accuracy, in an elementwise sense, that
matches non-asymptotic lower bounds applicable to any estimator.
Moreover, we showed how this theory enables data-driven construction
of confidence intervals for the missing potential outcomes. We note
that our development relies on a more general inferential toolbox for
the SVD algorithm in the matrix denoising model. It introduces the
idea of ``leave-one-block-out'', an analysis technique that we suspect
will be useful for other problems.
Moving forward, let us outline some open questions, along with broader
directions for future study. Our current theory relies on a noise
condition~\eqref{eq:noise-condition-est} and an incoherence
condition~\eqref{eq:incoherence}; it remains unclear whether or not
these conditions are improvable, for instance in terms of dependence
on the rank $r$. Moreover, it would be interesting to extend our
guarantees to accommodate matrices that are only approximately
low-rank, along with more general noise models. Finally, while the
the staggered adoption design arises frequently, it would be
interesting to understand optimal procedures for more general
mechanisms of missing data.
\subsection*{Acknowledgements}
Y.~Yan was supported by the Norbert Wiener Postdoctoral Fellowship
from MIT. M.~J.~Wainwright was partially supported by ONR grant
N00014-21-1-2842 and NSF grant DMS-2311072. We thank Eric Xia for
helpful discussions.
\printbibliography