Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
125,522 characters · 31 sections · 89 citation commands
Causal Matrix Completion
Matrix completion is the study of recovering an underlying matrix from its noisy and partial observations. Given its widespread applicability, the field of matrix completion has grown tremendously in recent years. To establish statistical guarantees for the various algorithms that exist for matrix completion, it is typically assumed that: (i) the underlying noiseless matrix has latent structure, e.g., it is low-rank, and (ii) the entries of this matrix are missing completely at random (MCAR), i.e., an entry is missing independent of everything else and with uniform probability. However, numerous modern applications of interest violate the latter assumption. Below, we consider two motivating examples.
First, arguably the most well-known application of matrix completion is recommender systems, which are ubiquitous in modern online platforms. Typically, data is collected in the form of a matrix, where the rows index users and columns index items; the $(i,j)$-th entry, therefore, corresponds to the rating supplied by user $i$ for item $j$. In such scenarios, observations are often subject to {\em selection-biases}. For instance, in movie recommendations, a fan of fantasy fiction will almost certainly watch and highly rate the {\em Harry Potter} series. Similarly, in restaurant recommendations, a vegetarian is unlikely to enjoy nor rate a steakhouse restaurant. While these examples demonstrate self-selection biases from the end of the users, systems also exhibit targeted suggestions. For example, when a user searches for trails at the Grand Canyon, an ad placement system is more likely to display an ad for hiking boots than wedding shoes; in turn, this can increase the user's likelihood to purchase and rate hiking boots. In all of these cases, the user's preferences and/or the system's beliefs in its users' preferences, influence the sparsity pattern of the observation matrix.
A second example is panel data settings in econometrics. Here, observations of units (e.g., individuals, geographic locations) are collected over time as they undergo different interventions (e.g., promotions, socio-economic policies). The induced matrix has rows index units and columns index time-intervention pairs; the $(i, (a, t))$-th entry then corresponds to the potential outcome of unit $i$ under the $a$-th intervention at time step $t$; here $(a, t)$ represents the $j$-th column, i.e., columns are double indexed by both intervention and time. As with recommender systems, observations in panel data settings are unlikely to occur completely at random. For instance, policy-makers strategically recommend programs that are designed to achieve certain desirable outcomes based on numerous socio-economic factors surrounding the geographic region under their purview. Further, competing programs with disagreeing agendas cannot be simultaneously adopted for a specific region during the same time period, i.e., if the $(i, (a, t))$-th entry is observed, then the $(i, (a', t))$-th entry must be missing. Notably, similar matrices and observation patterns can arise in sequential decision-making paradigms within machine learning such as online learning, contextual bandits, and reinforcement learning with time-intervention pairs being replaced by state-action pairs.
In both examples, the missingness pattern of the matrix is dependent on the underlying values in that matrix, and observing the outcome of one entry can alter the probability of observing another. That is, the entries are missing {\em not} at random (MNAR). To address the above challenges, there has been exciting recent progress on matrix completion with MNAR data, including schnabelfwang16, ma2019missing, zhu2019high, sportisse2020imputation_low_rank, sportisse2020estimation_PCA, causal_recommender_systems, yang2021tenips, bhattacharya2021matrix. Through numerous empirical studies, these works have shown that algorithms that account for MNAR data outperform conventional algorithms that are designed for MCAR data. With respect to theoretical analysis, however, critical aspects of matrix completion with MNAR data remain to be explored. In particular, as highlighted in ma2019missing, there are two common limiting assumptions in the literature: (i) the revelation of each entry in the matrix is independent of all other entries, and (ii) each entry has a nonzero probability of being observed.
Another recent exciting line of work that we build upon is that of panel data and matrix completion, see RSC, mRSC, arkhangelsky2019synthetic, bai2019matrix, fernandez2020low, athey2021matrix, agarwal2019robustness,agarwal2020principal, agarwal2021synthetic, agarwal2021causal. Some of these works allow for MNAR data and entries of a matrix to be deterministically missing. However, they consider very restricted sparsity patterns that are not particularly suitable for important applications of matrix completion. For example, the most common sparsity pattern considered in the panel data literature is where for a given row $i$, if a column $j$ is missing, then entries for all columns $j' > j$ in row $i$ are also missing; such a pattern is unlikely to arise in recommendation systems or sequential decision-making. Further, to the best of our knowledge, none of these works within the panel data literature provide meaningful results for matrix completion with MCAR data. The statistical parameters these works aim to estimate are also less meaningful for these other applications of matrix completion. The most common statistical parameter these works consider is the average outcome for all missing entries in a given row $i$; in say recommendation systems, this would correspond to the average rating a user $i$ would have given for all movies they did not rate. This is not particularly meaningful for an online platform---ideally, a platform would like to do accurate inference for each $(i, j)$ pair.
The focus of this work is to propose a formal causal framework and an algorithm with provable guarantees to analyze matrix completion with MNAR data where the probability that an entry of the matrix is missing can: (i) depend on the underlying values in the matrix itself; (ii) depend on which other entries are missing; (iii) potentially be deterministically zero. Further, we want to allow for more general missingness patterns and estimate more refined statistical parameters than considered in the panel data literature thus far.
Indeed, it is both commonly said that
We hope this work further bridges the rich and growing fields of causal inference and matrix completion.
As further motivation for why it is important to carefully think about the underlying mechanism for why data is missing, we now provide illustrative empirical simulations. In particular we run three experiments, each with a different mechanism for how data is missing. In Experiment 1, data is missing via a MCAR mechanism i.e., each entry is missing independently at random with probability $0.35$; the induced sparsity pattern is depicted in Figure (ref). In Experiment 2, data is missing in a MNAR fashion, i.e., each entry has a different probability of being missing; the induced sparsity pattern is depicted in Figure (ref). However, we ensure key assumptions made thus far in the matrix completion literature with MNAR data are maintained; in particular, (i) the revelation of entries are entry-wise independent and (ii) each entry has a nonzero probability of being observed. In Experiment 3, data is missing in a MNAR fashion, but we violate conditions (i) and (ii) above; the induced sparsity pattern is depicted in Figure (ref). For exact details on the missingness mechanism in Experiment 2 and 3, refer to Section (ref) and (ref), respectively.
In all experiments, we first create a sample of true “ratings”, which are invariant across all three experiments. We enforce these ratings to go from 1 to 5, as is standard in many online platforms. The distribution of true/revealed ratings are plotted in light/dark blue in Figures (ref), (ref), and (ref), respectively. As expected, the distribution of the revealed ratings in the MCAR setup matches that of the true ratings. However, the set of ratings that are revealed in both MNAR settings are severely biased, i.e., their distribution does not match that of the true underlying ratings.
We use three matrix completion algorithms and see whether they can recover the distribution of true ratings given the revelead entries in all three experiments. The algorithms are: Universal singular value thresholding (USVT) Chatterjee15, which is a popular spectral based method;\footnote{Surprisingly, we find that the original USVT algorithm performs better in all three experiments. See Appendix (ref).} Softimpute (softImpute) softimpute14, which is a popular optimization based method; Synthetic nearest neighbours (SNN), which is our proposed method for matrix completion with MNAR data, and is a combination of the approach taken in nearest neighbour style and panel data methods in econometrics. USVT and softImpute are not designed for MNAR data, as is, but we de-bias them for MNAR data as is done in bhattacharya2021matrix and ma2019missing, respectively. See details in Section (ref).
We see that in Figure (ref), under the MCAR setting, softImpute and SNN both recover the distribution of true ratings very well, while USVT cannot. Once we go to the limited MNAR setting, depicted in Figure (ref), where conditions (i) and (ii) are upheld, SNN is still able to recover the underlying distribution of true ratings, but now both softImpute and USVT have non-negligible bias. In the general MNAR setting, depicted in Figure (ref), \texttt{SNN} continues to accurately recover the distribution, but the bias of \texttt{softImpute} is significantly worsened.
This empirical illustration highlights the sensitivity of these traditional matrix completion methods to the missingness mechanism and strongly motivates the need for a rigorous framework for tackling the general MNAR setting where conditions (i) and (ii) above are violated. Providing such a framework is what we set out to do in this work.
We now formally introduce our setup. Consider a signal matrix $\bA = [A_{ij}] \in \Rb^{m \times n}$, a noise matrix $\bE = [\varepsilon_{ij}] \in \Rb^{m \times n}$, and a propensity score matrix $\bP = [p_{ij}] \in [0,1]^{m \times n}$. All three matrices are entirely latent, i.e., unobserved. Let $\boldsymbol{Y} = [Y_{ij}] \in \Rb^{m \times n}$ denote the “noisy” version of $\bA$, with $\mathbb{E}[\boldsymbol{Y}] = \bA$; we denote $\varepsilon_{ij} = Y_{ij} - A_{ij}$. We assume $\boldsymbol{Y}$ itself is partially observed. In particular, we denote $\bD = [D_{ij}] \in \{0,1\}^{m \times n}$ with $\mathbb{E}[\bD] = \bP$ as the missingness mask matrix that indicates which entries of $\boldsymbol{Y}$ are observed. For convenience, we encode our observations into $\widetilde{\boldsymbol{Y}} = [\widetilde{Y}_{ij}] \in \{\Rb \cup \{\star\}\}^{m \times n}$ such that for $(i,j) \in [m] \times [n]$,
In words, if $D_{ij}=1$ then $A_{ij}$ is noisily observed, and if $D_{ij}=0$ then $A_{ij}$ remains unknown. For concreteness, let us return to the recommender system example. Here, $\bA$ represents the expected rating for every user-item pair and $\bP$ dictates the probability these expected ratings are revealed, both of which are unknown. $\boldsymbol{Y}$ in relation to $\bA$ then models the inherent randomness in how users rate items; that is, $\boldsymbol{Y}$ can be interpreted as a “noisy” instance of $\bA$. Another interpretation of what $\varepsilon_{ij}$ represents is that many online platform only allow users to input integer valued ratings (e.g. integer between 1 to 5 or a binary 0/1). Hence, $Y_{ij}$ can be interpreted as a “noisy” {\em discretized} observation of $A_{ij}$, which may actually be {\em continuous} (i.e., lie within the continuous interval $[1,5]$ or $[0, 1]$). Observationally, we have access to $\bD$ and $\widetilde{\boldsymbol{Y}}$; the former refers to the collection of ratings users have supplied to the system while the latter refers to the corresponding realized “noisy” ratings. Finally, we remark that (ref) also agrees with standard panel data setups in econometrics, where each observation is assumed to be corrupted by an idiosyncratic shock, which is represented by $\varepsilon_{ij}$.
In terms of the type of MNAR data this work considers, we allow for $\bD$ and $\boldsymbol{Y}$ to be dependent, provided $\bD \protect\mathpalette{\protect\independenT}{\perp} \boldsymbol{Y} | \bA$, where $\bA$ is latent. In fact, we allow $\bD$ to be any arbitrary function of $\bA$, random or deterministic, subject to suitable observation patterns which we discuss in the forthcoming sections. Notably, our framework also allows the entries in $\bD$ to be dependent with each other across both rows and columns, and the minimum value of $\bP$ to be $0$, which are important departures from the current matrix completion literature. Under these conditions, we propose an algorithm that provably recovers $\bA$ from $\widetilde{\boldsymbol{Y}}$ with entry-wise (i.e., max-norm) guarantees.
{\bf Section (ref): Related works.} We provide an overview of the current literature on matrix completion under the different models of missingness proposed by rubin1976, little2019statistical: (i) missing completely at random (MCAR); (ii) missing at random (MAR); (iii) missing not at random (MNAR). We note the usage of the terms MCAR, MAR, and MNAR is inconsistent across the previous works on matrix completion, and so we hope that our literature survey helps give a more comprehensive and unified overview of the different regimes of missingness considered in these works.
{\bf Section (ref): Causal framework for matrix completion.} We propose a formal causal framework for matrix completion using the language of potential outcomes, see neyman, rubin. We interpret $\boldsymbol{Y}$ as the matrix of potential outcomes and $\bP$ as the matrix of intervention assignments. Building upon the recent work of agarwal2021synthetic, we propose a framework that allows (i) correlation between $\bD$ and $\boldsymbol{Y}$, i.e., hidden confounding; (ii) correlation between the entries of $\bD$; (iii) the minimum value of $\bP$ to be $0$, i.e., entries of $\widetilde{\boldsymbol{Y}}$ can be deterministically missing; (iv) $\bP$ to not exhibit low-dimensional structure as is required in the panel data literature, i.e., we consider significantly more general missingness patterns. To the best of our knowledge, our framework, and associated algorithm, is the first within the MNAR matrix completion literature that allows for conditions (i)-(iv) to simultaneously hold. Additionally, we do {\em not} make any parametric or distributional assumptions on $\bP$, as is common in previous works on matrix completion. Nevertheless, we establish an identification result in Theorem (ref), which effectively states that $\bA$ can be learned from $\widetilde{\boldsymbol{Y}}$ in an entry-wise sense. We believe our proposed framework provides a unified causal view for a variety of applications that can be posed as matrix completion problems with MNAR data.
{\bf Section (ref): An algorithmic solution.} We combine the nearest neighbours approach for matrix completion ---popularly known as collaborative filtering---with the synthetic controls approach for panel data, to design a novel two-step algorithm, which we call “synthetic nearest neighbors” (SNN), to estimate $\bA$ from $\widetilde{\boldsymbol{Y}}$. Pleasingly, each step of SNN enjoys a simple closed-form solution. In order to efficiently execute SNN in practice, we provide an algorithm to automatically find the “neighbors” for any $(i,j)$ pair in a data-driven manner. To do so, we relate this task to the well-known problem of finding the “maximum” biclique in a bipartite graph. Since SNN is a generalization of the recently proposed synthetic interventions (SI) estimator of agarwal2021synthetic, which itself is a generalization of the popular synthetic controls algorithm of abadie1, abadie2, this subroutine may be of independent interest to the synthetic controls and panel data literatures.
{\bf Section (ref): Theoretical results.} We establish entry-wise finite-sample consistency and asymptotic normality of SNN, i.e., we provide theoretical guarantees for $A_{ij}$ for each $(i,j)$ pair. Hence, our analysis implies new theoretical results, in a max-norm sense, for the literature on matrix completion with MNAR data. As a special case, this also provides novel entry-wise finite-sample consistency and asymptotic normality results for the traditional matrix completion with MCAR data literature. Collectively, our identification, consistency, and asymptotic normality results, coupled with SNN, can be seen as a generalization of the SI framework proposed in agarwal2021synthetic.
{\bf Section (ref): Experimental validation.} We run comprehensive experiments, both with simulated and real-world data, to test the empirical efficacy of SNN against a collection of state-of-the-art matrix completion algorithms for MNAR data. Some key takeaways are as follows: (i) SNN is robust to the various forms of missingness across all experiments, while the previous methods are relatively sensitive to it. (ii) we find the approaches to de-bias estimators for MNAR data are not particularly effective, i.e., their performance is similar to their MCAR analogues; this is in line with the empirical findings of ma2019missing.
For a matrix $\bX \in \Rb^{m \times n}$, we denote its operator (spectral), nuclear, Frobenius, and max element-wise norms as $\|\bX\|_{2}$, $\|\bX\|_*$, $\|\bX\|_F$, and $\|\bX\|_{\max}$, respectively. For a matrix $\bX$ with orthonormal columns, let $\Pc_X = \bX \bX^T$ denote the projection matrix onto the subspace spanned by the columns of $\bX$. For a vector $v \in \Rb^m$, let $\|v\|_p$ denote its $\ell_p$-norm. For a random variable $v$, we define its sub-gaussian (Orlicz) norm as $\|v\|_{\psi_2}$. Let $\circ$ denote component-wise multiplication and let $\otimes$ denote the outer product. For a positive integer $a$, let $[a] = \{1, \dots, a\}$. For index sets $\Ic_1 \subseteq [m] $ and $\Ic_2 \subseteq [n]$, let $\bX_{\Ic_1, \Ic_2}$ denote the $| \Ic_1| \times | \Ic_2|$ sub-matrix of $\bX$ whose rows and columns are indexed by $\Ic_1$ and $\Ic_2$, respectively. As a shorthand, let $\bX_{\Ic_1, \cdot}$ denote the $|\Ic_1| \times n$ sub-matrix of $\bX$ that retains the columns of $\bX$ but only considers those rows indexed by $\Ic_1$; we define $\bX_{\cdot, \Ic_2}$ analogously. Unless stated otherwise, we index rows with $i \in [m]$ and columns with $j \in [n]$.
Let $f$ and $g$ be two functions defined on the same space. We say $f(n)$ = $O(g(n))$ if and only if there exists a positive real number $M$ and a real number $n_0$ such that for all $n \ge n_0, |f (n)| \le M|g(n)|$. Analogously we say: $f (n) = \Theta(g(n))$ if and only if there exists positive real numbers $m, M$ such that for all $n \ge n_0, \ m|g(n)| \le |f(n)| \le M|g(n)|$; $f (n) = o(g(n))$ if for any $m > 0$, there exists $n_0$ such that for all $n \ge n_0, |f(n)| \le m|g(n)|$. We adopt the standard notations and definitions for stochastic convergences. As such, we denote $\xrightarrow{d}$ and $\xrightarrow{p}$ as convergences in distribution and probability, respectively. We will also make use of $O_p$ and $o_p$, which are probabilistic versions of the commonly used deterministic $O$ and $o$ notations. More formally, for any sequence of random vectors $X_n$, we say $X_n = O_p(a_n)$ if for every $\varepsilon>0$, there exists constants $C_\varepsilon$ and $n_\varepsilon$ such that $\Pb( \| X_n \|_2 > C_\varepsilon a_n) < \varepsilon$ for every $n \ge n_\varepsilon$; equivalently, we say $(1/a_n) X_n$ is “uniformly tight” or “bounded in probability”. Similarly, $X_n = o_p(a_n)$ if for all $\varepsilon, \varepsilon'> 0$, there exists $n_\varepsilon$ such that $\Pb(\| X_n \|_2 > \varepsilon' a_n) < \varepsilon$ for every $n \ge n_\varepsilon$. Therefore, $X_n = o_p(1) \iff X_n \xrightarrow{p} 0$. Additionally, we denote: $\text{plim} \ X_n = a \iff X_n \xrightarrow{p} a$. We say a sequence of events $\Ec_n$, indexed by $n$, holds “with high probability” (w.h.p.) if $\Pb(\Ec_n) \rightarrow 1$ as $n \rightarrow \infty$, i.e., for any $\varepsilon > 0$, there exists a $n_\varepsilon$ such that for all $n > n_\varepsilon$, $\Pb(\Ec_n) > 1 - \varepsilon$. More generally, a multi-indexed sequence of events $\Ec_{n_1,\dots, n_d}$, with indices $n_1,\dots, n_d$ with $d \geq 1$, is said to hold w.h.p. if $\Pb(\Ec_{n_1,\dots, n_d}) \rightarrow 1$ as $\min\{n_1,\dots, n_d\} \rightarrow \infty$. We also use $\mathcal{N}(\mu, \sigma^2)$ to denote a normal or Gaussian distribution with mean $\mu$ and variance $\sigma^2$---we call it {\em standard} normal if $\mu = 0$ and $\sigma^2 = 1$.
Given the vastness of the matrix completion literature, we do not strive to do an exhaustive review of it. Instead, we focus on a few representative works that propose and analyze algorithms designed for the three different models of missingness: MCAR, MAR, and MNAR. In Section (ref), we give an overview of the type of algorithms for matrix completion studied thus far in existing works. In Section (ref), we discuss the different models of missingness considered in the matrix completion literature, and representative algorithms for these various models. Finally, in Section (ref), we discuss the growing literature exploring the intersection of matrix completion and causal inference; in particular, the panel data literature in econometrics.
Algorithms for matrix completion broadly fall into two classes: empirical risk minimization (ERM) methods and matching (i.e., collaborative filtering) methods, with ERM methods being relatively more popular. We give an overview of both class of methods below.
{\bf Empirical Risk Minimization (ERM) Methods.} Empirical risk minimization (ERM) is arguably the de facto approach to recover the underlying signal matrix $\bA$ given $\widetilde{\boldsymbol{Y}}$. Specifically, ERM approaches aim to solve the following program:
Here, $\Omega \subseteq [m] \times [n]$, $d(\cdot, \cdot)$ is an appropriate distance measure (e.g., squared loss), $T_{ij}$ is a “simple” transformation of $\widetilde{Y}_{ij}$ (e.g. $\mathds{1}(D_{ij}=1) \cdot \widetilde{Y}_{ij}$), $\textsf{regularize}(\cdot)$ is a regularization term and $\lambda > 0$ is the regularization hyper-parameter. For certain algorithms, they replace the regularizer (i.e., set $\lambda = 0$) with a constraint, $\textsf{constraint}(\cdot)$.
In order to prove statistical guarantees about these various estimators, structure is placed on $\bA$. The assumptions made guide the specific choices of the above parameters, which then define the algorithm. For instance, if the singular values of $\bA$ are assumed to be moderately sparse (i.e., only few are non-zero), then a natural convex regularizer would penalize solutions with large nuclear norm, i.e., $\textsf{regularize}(\bQ) = \| \bQ \|_*$ CandesTao10, Recht11. Indeed, choosing $\Omega = \{(i,j): D_{ij} = 1\}$ as the collection of observed entries, $T_{ij} = \widetilde{Y}_{ij}$, and $d(\cdot, \cdot)$ as the squared loss yields the popular softImpute algorithm of softimpute, softimpute14. As another example, if $\bA$ is assumed to be exactly low-rank, then a natural constraint would be the rank of the output matrix. More specifically, $\textsf{contraint}(\bQ)$ can be defined as $\text{rank}(\bQ) \le \mu$ for some pre-specified integer $\mu > 0$. Then, choosing $\Omega = [m] \times [n]$, $T_{ij} = \mathds{1}(D_{ij}=1) \cdot \widetilde{Y}_{ij} $, and $d(\cdot, \cdot)$ as the squared loss yields a suite of spectral based methods KeshavanMontanariOh10a, KeshavanMontanariOh10b, donoho14, Chatterjee15. Other notable algorithms within the broader ERM class include maximum-margin matrix factorization (MMMF) mmmf, probabilistic matrix factorization (PMF) pmf, and SVD++ koren08 to name a few.
Broadly speaking, it is commonly assumed that $\bA$ follows some form of a latent variable model; in particular, $A_{ij} = f(u_i, v_j)$, where $f$ is a sufficiently “smooth” latent function (e.g., H\"older continuous), and $u_i, v_j$ are low-dimensional latent variables associated with row $i$ and column $j$, respectively. Such latent variable models imply that $\bA$ is (approximately) low-rank, i.e., $A_{ij} \approx \langle u_i, v_j \rangle$, where $u_i, v_j \in \Rb^r$ and $r \ll \min\{m, n\}$, e.g., xu2017rates, udell_low_rank, agarwal2019robustness. For an excellent overview on standard assumptions made on $\bA$ and the subsequent guarantees proven for the estimation error, please refer to ieee_matrix_completion_overview.
When every entry is revealed with uniform probability (i.e., $p_{ij} = p$), (ref) is an unbiased estimate of the full loss function with all entries revealed (i.e., $\bD$ is an all ones matrix). When $p_{ij}$ are nonuniform, however, recent works have provably and empirically shown that (ref) is biased schnabelfwang16, ma2019missing. As such, these works advocate to de-bias the standard ERM objective by re-weighting each observation inversely by its propensity score $p_{ij}$. This technique is often known in the causal inference literature as inverse propensity scoring (IPS) or weighting (IPW), see imbens_rubin_2015, little2019statistical. This yields the following adapted program:
where $\widehat{p}_{ij}$ is an estimate of $p_{ij}$. In words, (ref) requires learning $\bP$ prior to carrying out the standard ERM of (ref). Faithful matrix recovery under more general missingness patterns thus requires structure on not only $\bA$, but also $\bP$ and $\bD$. We overview standard assumptions on these quantities in Section (ref).
{\bf Matching methods.} For traditional applications of matrix completion, such as recommendation systems, K nearest neighbour (KNN) methods have been popular (e.g., goldberg1992using, linden2003amazon, kleinberg2008using, koren2015advances, LeeLiShahSong16, chen2018explaining). In KNN, to impute a missing entry $(i, j)$, the first step is to select $K$ rows for which the entry in the $j$-th column is not missing. Of all the rows for which the $j$-th column is not missing, the $K$ rows are selected such that they are the “closest” to row $i$. In particular, a hyper-parameter of KNN is the metric that is chosen to define “closeness” between any two given rows; the most commonly used metric is the mean squared distance between the commonly revealed entries for a given two rows. Once these $K$ “neighbour rows” are chosen, the estimate for the missing entry $(i, j)$ is the average $\frac{1}{K}\sum_{k \in \text{neighbour rows}} \widetilde{Y}_{k j}$. An attractive quality of these KNN methods is that they do not require imputing missing values by $0$. A related literature that shares similarities with KNN is that of synthetic controls abadie1, abadie2. A key difference is that to impute $(i, j)$, uniform weights (i.e., $1/K$) are not used for the neighbouring rows; classically in synthetic controls, these weights are constrained to lie within the simplex, i.e., the weights are non-negative and sum to $1$ (if the weights are restricted to be $1/K$, this is known in the panel data literature as “difference-in-differences”). However, as discussed earlier, synthetic controls methods have been designed to handle restricted sparsity patterns naturally arising in the panel data setting. Given the growing literature on synthetic controls, we do a detailed literature review of it in Section (ref).
Below, we utilize the useful taxonomy set in rubin1976, little2019statistical to discuss the three primary mechanisms that lead to missing data and how previous works fit within these regimes.
{\bf Missing completely at random (MCAR).} MCAR is the most standard model of missingness assumed in the matrix completion literature and is characterized by the following properties: (i) $\bD \protect\mathpalette{\protect\independenT}{\perp} \boldsymbol{Y}$; (ii) $D_{ij} \protect\mathpalette{\protect\independenT}{\perp} D_{ab}$ for all $(i,j) \neq (a,b)$; (iii) $p_{ij} = p > 0$ for all $(i,j)$. In words, MCAR assumes each element of $\bD$ is an independent and identically distributed (i.i.d.) Bernoulli random variable (r.v.) with parameter $p \in (0, 1]$. This implies that the missingness pattern is independent of the values in $\boldsymbol{Y}$. We note that this condition $p_{ij} > 0$ is known in the causal inference literature as “positivity”, see imbens_rubin_2015. It follows that the maximum likelihood estimator $\widehat{p}_{ij} = \widehat{p}$ for all $(i,j)$, where $\widehat{p}$ is the fraction of observed entries in $\widetilde{\boldsymbol{Y}}$. As previously mentioned, given MCAR data, (ref) is an unbiased estimator of the ideal loss function where all entries observed. Though MCAR is likely unrealistic outside experimental settings, the MCAR regime remains a popular abstraction in machine learning and statistics to study the inherent trade-offs between the observation probability $p$, properties of the noise $\bE$, and the structure imposed on the signal $\bA$, in terms of the estimation error between $\bhA$ and $\bA$. Methods such as singular value thresholding explicitly impute missing values in $\boldsymbol{Y}$ (denoted as $\star$) by $0$ and re-weight all non-missing values in $\boldsymbol{Y}$ by $1 / \widehat{p}$, where $\widehat{p}$ is the fraction of observed entries. This can be interpreted as a form of uniform IPW. Other methods such as nuclear norm minimization, alternating least squares, and nearest neighbour methods do not require imputing missing values by $0$. However, existing theoretical analysis of these algorithms do still require that $\mathbb{E}[D_{ij}] = p$, and that $D_{ij}$ is independent of the all other randomness in the model.
{\bf Missing at random (MAR).} \footnote{Many works in the matrix completion literature do not differentiate between MAR and MNAR, and call both regimes MNAR. We differentiate between them to be more in line with models of missingness proposed by rubin1976, little2019statistical.} MAR is a more challenging setting than MCAR. The three key assumptions of MAR are as follows. (i) $\bD \protect\mathpalette{\protect\independenT}{\perp} \boldsymbol{Y} ~|~ \Oc$, where $\Oc$ represents observed covariates about the rows and columns of the matrix (e.g., covariates about users and movies in the context of recommender systems)---concretely, these observed variables, $\Oc$, often include features or covariates $(X_i, \tilde{X}_j)$, which are associated with row $i$ and column $j$, respectively, and observed outcomes $\widetilde{Y}_{ij}$. (ii) $D_{ij} \protect\mathpalette{\protect\independenT}{\perp} D_{ab}$ for all $(i,j) \neq (a,b)$. (iii) $p_{ij} > 0$ for all $(i,j)$. Here, the entries of $\bD$ continue to obey positivity and remain independent Bernoulli r.v.'s.
Below, we overview two popular propensity estimation techniques of schnabelfwang16. To aid the following discussion, let $\bX = \{(X_i, \tilde{X}_j): (i, j) \in [m] \times [n]\}$ denote the set of observed features, and $\bH$ denote the set of hidden features. The first approach is via Naive Bayes, which assumes that $p_{ij} = \mathbb{E}[D_{ij} | \bX, \bH, \boldsymbol{Y}] = \mathbb{E}[D_{ij} | \widetilde{Y}_{ij}]$. Under this assumption, the maximum likelihood estimator $\widehat{p}_{ij}$ can be solved using Bayes formula; however, such an approach requires a small sample of MCAR data, see schnabelfwang16. The second estimation strategy is based on logistic regression. Here, it is assumed that there exists model parameters $\phi$ such that $p_{ij} = \mathbb{E}[D_{ij} | \bX, \bH, \boldsymbol{Y}] = \mathbb{E}[D_{ij} | X_i, \tilde{X}_j, \phi]$; within the causal inference literature, this is often known as “selection on observables”, see imbens_rubin_2015. Typically, it is posited that $\phi = (\omega_1, \omega_2, \alpha, \gamma)$ and $\mathbb{E}[D_{ij} | X_i, \tilde{X}_j, \phi] = \sigma(\langle \omega_1, X_i \rangle + \langle \omega_2, \tilde{X}_j \rangle + \alpha_i + \gamma_j)$, where $\sigma(\cdot)$ takes a simple parametric form such as the sigmoid function. Some notable works in the MAR literature include Liang2016, WangNeurips2018, WangAAAI2018, wang2019.
{\bf Missing not at random (MNAR).} MNAR is the most challenging missingness model in matrix completion with a comparatively sparser literature. In its fullest generality, in MNAR the following conditions are allowed: (i) $\bD$ can depend on $\boldsymbol{Y}$ and other unobserved variables; (ii) $D_{ij}$ can be correlated with $D_{ab}$ for all $(i,j) \neq (a,b)$; (iii) $\min p_{ij} = 0$. The first condition implies that $\bD$ and $\boldsymbol{Y}$ remain dependent even conditional on observed covariates. The second condition allows the revelation of one outcome to alter the probability of another outcome being revealed. Finally, the third condition can restrict certain outcomes from ever being revealed. Hence, the literature has thus far only considered a {\em limited} version of MNAR with conditions cf. ma2019missing, bhattacharya2021matrix, yang2021tenips. In particular, they continue to make the following assumptions: $p_{ij}$ is a (nice) function solely of latent factors associated with entry $(i, j)$; each entry of $\bD$ is an independent (not necessarily identically distributed) Bernoulli r.v. with a strictly positive probability of being revealed, which are the assumptions as in MAR. These assumption are what allow the weighted ERM framework of (ref) to continue being valid. The methods proposed in ma2019missing, bhattacharya2021matrix, yang2021tenips work for this {\em limited} MNAR setting by positing that $\bP$ is (approximately) low-rank, and recovers $\bP$ from $\bD$ via matrix completion algorithms. This is a generalization of the MAR setting as such an approach circumvents the requirement of meaningful auxiliary features $\bX$ to conduct propensity score estimation. Additional works within the MNAR literature include zhu2019high, sportisse2020imputation_low_rank, sportisse2020estimation_PCA, causal_recommender_systems.
As previously mentioned, our work operates under greater generality than the {\em limited} MNAR regime thus far considered in the literature. More specifically, our framework allows $\bD$ and $\boldsymbol{Y}$ to be dependent, provided $\bD \protect\mathpalette{\protect\independenT}{\perp} \boldsymbol{Y} | \bA$, and for $\bD$ to be any arbitrary function of $\bA$, subject to suitable observation patterns. We also allow for conditions (ii) and (iii) described above to hold, i.e., the entries in $\bD$ can be highly correlated and the minimum probability of observation can be deterministically set to $0$. In Section (ref), we will formally introduce our causal framework to rigorously discuss these properties.
{\bf Summary of matrix completion results.} Across the various models of missingness, the key theoretical results for low-rank matrix completion typically have error bounds that scale in the following form (see ieee_matrix_completion_overview):
for $\delta \geq 0$ and where $\text{poly}(\cdot)$ denotes polynomial dependence. Here, $p_\min = \min p_{ij}$ and $r$ refers to the (approximate) rank of $\bA$. The most studied metric in the literature is the average error across all entries, $(1/mn) \| \bhA - \bA \|^2_F$, though recent works have begun to analyze stronger metrics such as the maximum average error across all columns, $(1/m) \| \bhA - \bA \|^2_{2, \infty}$ (e.g., agarwal2019robustness, agarwal2020principal, agarwal2021causal), and the maximum entry-wise error, $\| \bhA - \bA \|_\max$ (e.g., LeeLiShahSong16). Crucially, all of these error bounds scale with the inverse of $\text{poly}(p_\min)$. As discussed above, this immediately rules out settings where $p_\min = 0$, i.e., condition (iii) of MNAR above. Finally, we remark that the literature studying the asymptotic properties of $\| \bhA - \bA \|$ (e.g., proving asymptotic normality) is relatively small. Some notable works on the asymptotic analyses of matrix completion estimators under MCAR include chen2019inference, cai2020uncertainty, bhattacharya2021matrix.
In Section (ref), we propose a causal framework for matrix completion that draws inspiration from the rich and growing literature in econometrics on panel data and matrix completion; some relevant works include RSC, mRSC, arkhangelsky2019synthetic, bai2019matrix, fernandez2020low, athey2021matrix, agarwal2019robustness, agarwal2021synthetic, agarwal2021causal. As is common in matrix completion, these works impose a (approximate) low-rank factor model on the signal matrix (i.e., $\bA$), also known as an interactive fixed effects model, to capture structure across units and time (i.e., the rows and columns of the matrix, respectively).
{\bf Panel data & matrix completion: an overview.} As described in Section (ref), the sparsity structure considered in these works is one where for each row $i$, there is a column $j_i \in [n]$ such that $D_{i j} = 1$ for all $j < j_i$ and $D_{i j} = 0$ for $j \geq j_i$. That is, all entries for a given row $i$ are observed till some column $j_i$, after which they are all missing. The motivation for such a sparsity pattern comes from socio-economic policy making where $Y_{ij}$ represents unit $i$'s potential outcome at time step $j$ under “control”, i.e., if no socio-economic intervention has yet been applied on unit $i$. The time steps $[1, j_i -1]$ represent the period when unit $i$ is under control, and time steps $[j_i, n]$ represent the period when unit $i$ has undergone an intervention. Hence, $Y_{ij}$ for $j > j_i$ is missing and the goal is to estimate the counterfactual of what would have happened to unit $i$ had it remained under control during $[j_i, n]$. This particular setting is also known in the econometrics literature as “synthetic contorls” abadie1, abadie2. The statistical/causal parameter that is most commonly studied is for a “treated” unit $i$, to estimate $ \frac{1}{n - j_i} \sum^n_{j = j_i} Y_{ij}. $ That is, the average potential outcome of unit $i$ under control during the “post-intervention” period. Most of these works make the additional assumption that each unit either remains under control for the entire time period under consideration, or undergoes an intervention at a time step that is common across all units. athey2021matrix is one notable work that allows for different post-intervention periods for each unit.
{\bf Connections to matrix completion with MNAR data.} An attractive quality of this literature is that in some ways it allows for more relaxed conditions on $\bD$ and $\bP$ than those considered in the matrix completion with MNAR data literature discussed earlier, see ma2019missing, yang2021tenips, sportisse2020estimation_PCA, sportisse2020imputation_low_rank, wang2019. In particular, the panel data literature allows the entries of $\bD$ to be correlated, e.g., if $D_{i j} = 0$, then $D_{i j'} = 0$ for $j' > j$. Further, $\min p_{ij}$ is allowed to be $0$ and $j_i$ is allowed to depend on $\bA$. On the other hand, the sparsity pattern considered in the panel data literature is far more restrictive compared to the works on matrix completion with MNAR data---as discussed above, in panel data settings, all columns for a given row are observed till a specific point, after which they are all missing (i.e., $D_{i j} = 1$ for all $j < j_i$ and $D_{i j} = 0$ for $j \geq j_i$). Note that this also implies that $\bP$ is low-rank. Such a sparsity pattern is unrealistic for many important applications for matrix completion, including recommendation systems and sequential decision-making. Further, it is not straightforward to see how the target statistical/causal parameter $ \frac{1}{n - j_i} \sum^n_{j = j_i} Y_{ij} $ is particularly meaningful outside the synthetic controls literature. Hence, our aim with this work is to combine the best of both worlds, where we: (i) allow entries of $\bD$ to be correlated; (ii) allow $\min p_{ij} = 0$; (iii) make no parametric assumptions about $\bP$; (iv) allow $\bP$ to not be low-rank; (v) allow for general missingness patterns in the matrix that includes MCAR data as a special case. Further the target parameter we aim to estimate (in expectation) is {\em each entry} $Y_{ij}$ for every $(i, j)$ pair. Also, by formally bridging the panel data literature to more classical applications of matrix completion such as recommendation systems, we hope this spurs further investigation into the unexplored connections between these two fields.
{\bf Comparison with synthetic interventions.} Our proposed framework framework builds upon the recent work of agarwal2021synthetic, called synthetic interventions (SI). SI is a causal inference method to do tensor completion with MNAR data, where the dimensions of the order-3 tensor of interest are units, measurements, and interventions. That is, an entry $Y_{ijd}$ of the tensor considered in SI refers to the potential outcome of the $i$-th unit, its $j$-th measurement, under the $d$-th intervention. Their setup can be made a special case of ours by effectively {\em flattening} the tensor into a matrix, where the rows of the induced matrix still correspond to units, but a column is a double index for a measurement and an intervention, i.e., the $(i, j, d)$-th entry of the tensor corresponds to the $(i, (j, d))$-th entry of the induced matrix. Given this simple reduction, we generalize the framework, algorithm, and theoretical results in agarwal2021synthetic in the following ways. First, we formally extend the SI framework, to recover matrices under more general missingness patterns than that considered in agarwal2021synthetic. Doing so allows us to apply our framework to a wider variety of applications such as recommender systems, while the SI framework was introduced in the context of personalized policy evaluation and synthetic A/B testing. Third, this work establishes point-wise finite-sample consistency and asymptotic normality of our proposed SNN algorithm, which was absent in agarwal2021synthetic with respect to the \textsf{SI} algorithm. Indeed, in the context of the panel data literature, establishing point-wise asymptotic normality for each unit, (intervention, time)-tuple is of independent interest.
In this section, we develop a formal causal framework for matrix completion with MNAR data. In Section (ref), we show how to causally interpret matrix completion with MNAR data using the language of potential outcomes in Section (ref). We then state and justify our assumptions in Section (ref), define our causal estimand in Section (ref), and present our identification result in Section (ref).
We follow the potential outcomes framework of neyman, rubin. In particular, we let the r.v. $Y_{ij} \in \Rb$, as defined in Section (ref), denote the potential outcome associated with each pair $(i,j)$ if it is revealed. For instance, in the case of recommender systems, $Y_{ij}$ can be interpreted as the rating user $i$ {\em would have} given to item $j$ had they rated it. In the context of healthcare for example, $Y_{ij}$ could represent patient $i$'s health metric of interest (e.g. heart rate) {\em had they} been given treatment $j$. Finally, in the case of panel data setting, as discussed in Section (ref), $Y_{i, (a, t)}$ can denote the metric of interest for unit $i$ (e.g. revenue generated, socio-economic indicator), if they {\em would have} received the $a$-th socio-economic policy at time step $t$; here, $(a, t)$ represents the $j$-th column.
If $D_{ij} = 1$, then by (ref) we see that we actually do observe the $(i,j)$-th potential outcome, i.e. $\widetilde{Y}_{ij} = Y_{ij}$. That is, in the language of potential outcomes, we can interpret $\bD$ as the matrix of intervention assignments. Through this perspective, we remark that (ref) is an implicit assumption that is known in the causal inference literature as “consistency” or “stable-unit-treatment-value assumption” (SUTVA). As discussed earlier, the fact that $\boldsymbol{Y} \not\perp\!\!\!\perp \bD$ (e.g. a user's preference for a movie can determine whether they rate it) means that the potential outcomes are not independent of the intervention assignments. This dependence is known in the causal inference literature as “confounding”. Lastly, as alluded to earlier, we generalize the standard potential outcomes framework in that a given unit can receive multiple interventions. Traditionally, it is assumed that a unit receives exactly one intervention. However, in applications like movie recommendation systems, a user can “intervene” and rate multiple movies. Lastly, this framework also generalizes panel data settings, as we allow each unit to receive different interventions at different time steps; as discussed earlier, it is typically assumed that units are in control for a period of time, and then some subset of units receive one intervention for the remaining time steps.
Below, we state our causal assumptions and then provide their corresponding interpretations.
{\em Neighbourhood rows and columns.} For the remainder of this work, for a given column $j$, we refer to $\textsf{NR}(j) = \{a \in [m]: D_{a j} = 1\}$ as “neighborhood rows”, i.e., rows where entries in column $j$ are not missing. Similarly, for a given row $i$, we refer to $\textsf{NC}(i) = \{b \in [n]: D_{ib} = 1\}$ as “neighborhood columns”, i.e., columns where entries in row $i$ are not missing. See Figure (ref) for a visual depiction of $\textsf{NR}(j)$ and $\textsf{NC}(i)$.
{\bf Interpretation of Assumptions (ref) to (ref)} By the tower law, Assumption (ref) implies that $\mathbb{E}[ \bE | \bU, \bV] = 0$. This together with Assumption (ref) posits that $\mathbb{E}[\boldsymbol{Y} | \bU, \bV]$ is a low-rank matrix with rank $r$. As discussed in Section (ref), this is a standard assumption within the matrix completion literature. Next, we remark that Assumption (ref), coupled with Assumption (ref), implies that
That is, the potential outcomes are mean independent of the intervention assignments, conditioned on the latent row and column factors. This has been termed as “selection on latent factors”, see agarwal2021synthetic. Similar conditional independence conditions have been explored in athey2021matrix, kallus2018causal. Lastly, given Assumption (ref), it follows that Assumption (ref) is rather mild. To see this, suppose $\text{span}(\{u_\ell: \ell \in \Ic\}) = \Rb^r$, i.e., $\text{rank}(\bU_{\Ic, \cdot}) = r$. Then, Assumption (ref) immediately holds as $u_i \in \Rb^r$. More generally, if the rows of $\bU$ are randomly sampled sub-gaussian vectors, then $\text{span}(\{u_\ell: \ell \in \Ic\}) = \Rb^r$ for any set $\Ic$ holds w.h.p., provided $\mu \ge r$ is chosen to be sufficiently large; see vershynin2018high for details.
Define
Note that given Assumptions (ref) and (ref), the definition of $\bA$ in (ref) is consistent with the definition of $\bA$ used in Section (ref). We are now equipped to define our target causal estimand, which is $A_{ij}$; for the remainder, of the paper we focus on a particular pair $(i, j)$, without loss of generality. Note given Assumptions (ref) and (ref), we can write
In words, $A_{ij}$ translates as the expected potential outcome for the $(i,j)$-th pair, conditioned on its row and column latent vectors $(u_i, v_j)$. For instance, returning to recommender systems, $A_{ij}$ represents the expected rating user $i$ would supply for item $j$, conditioned on the latent features that characterize user $i$ and item $j$. In panel data settings, letting $j = (a, t)$, $A_{ij}$ represents the potential outcome of unit $i$ had it received the $a$-th intervention at time step $t$.
The following identification results establishes that each entry of $\bA$ can be learned from observable quantities, i.e., from $\widetilde{\boldsymbol{Y}}$. Practically speaking, this means that matrix completion with MNAR data for any pair $(i, j)$ is possible.
{\bf Interpretation.} Theorem (ref) states that despite the missingness pattern being MNAR, if Assumptions (ref) to (ref) hold, and given knowledge of the linear model parameter $\beta$, the causal estimand $A_{ij}$ can be expressed in terms of quantities that can be estimated from observed data, namely $\mathbb{E}[\widetilde{\boldsymbol{Y}}_{\Ic, j}]$; this is known in the causal inference literature as “identification”. Note, $\Ic$ is deterministic given $\bD$. The key requirement of the missingness pattern $\bD$ is that $\Ic \subseteq \textsf{NR}(j)$ is sufficiently large, which is parameterized by $\mu$, i.e., we require $\mu \gg r$ where $r$ is the rank of $\bA$. That is, the number of rows for which column $j$ is observed is sufficiently large. Thus, Theorem (ref) suggests that the key quantity that enables the recovery of $A_{ij}$ is $\beta$. In Section (ref), we provide an algorithm to estimate $\beta$, which in turn, allows us to estimate $A_{ij}$.
In this section, we introduce an algorithm, synthetic nearest neighbors (SNN), for matrix completion with MNAR data. Towards this, we introduce helpful notation that will be used for the remainder of this work. Again, without loss of generality, we consider imputing the $(i,j)$-th entry of the matrix .
{\bf Notation.} Let $\textsf{AR} \subseteq \textsf{NR}(j)$ and $\textsf{AC} \subseteq \textsf{NC}(i)$ denote a subset of rows and columns, respectively, of $\widetilde{\boldsymbol{Y}}$ that satisfy $D_{ab} = 1$ for all $(a,b) \in \textsf{AR} \times \textsf{AC}$. We refer to $\textsf{AR}$ and $\textsf{AC}$ as the “anchor rows” and “anchor columns” of pair $(i,j)$, respectively. Collectively, $\textsf{AR}$ and $\textsf{AC}$ form a fully observed sub-matrix of $\widetilde{\boldsymbol{Y}}$; for ease of notation. We refer to this $|\textsf{AR}| \times |\textsf{AC}|$ sub-matrix as $\bS \coloneqq [\widetilde{Y}_{ab}: (a,b) \in \textsf{AR} \times \textsf{AC}]$. See Figure (ref) for a visual depiction of $\textsf{AR}$, $\textsf{AC}$, and $\bS$. Note, by construction $\bS$ is such that if entries from row $a$ are present in $\bS$, then $D_{a j} = 1$; similarly, if entries from column $b$ are present in $\bS$, then $D_{i b} = 1$ Additionally, let $q \coloneqq [\widetilde{Y}_{ib}: b \in \textsf{AC}]$ and $x \coloneqq [\widetilde{Y}_{aj}: a \in \textsf{AR}]$. $q \in \Rb^{|\textsf{AC}|}$ refers to the columns in row $i$ which correspond to $\textsf{AC}$; similarly, $x \in \Rb^{|\textsf{AR}|}$ refers to the rows in column $j$ which correspond to $\textsf{AR}$. By construction, all the elements in $q$ and $x$ are not missing. See Figure (ref) for a visual depiction of $q$ and $x$.
We now present SNN in Algorithm (ref) to impute the $(i,j)$-th entry. It has $K \in \Nb$ and $\lambda^{(k)} \in \Rb$ for $k \in [K]$ as hyper-parameters.
Note, for ease of notation, in Algorithm (ref) we suppress the dependence on $i$ and $j$ in the definitions of $\{(\textsf{AC}^{(k)}, \textsf{AR}^{(k)}): k \in [K]\}$, $\bS^{(k)}$, $\widehat{\beta}^{(k)}$, $q^{(k)}$, and $x^{(k)}$. That is, these quantities will change depending on which $(i,j)$-th entry of the matrix we aim to impute. We continue to suppress this dependence for the remainder of the paper. For a visual depiction of the SNN algorithm for $K = 1$, refer to Figure (ref). For $K > 1$, we simply re-run the SNN algorithm seperately for the $K$ disjoint subsets $\{\textsf{AR}^{(k)}: k \in [K]\}$, and take the average of the estimates $\hA_{ij}^{(k)}$ for $k \in [K]$, produced by each iteration.
{\bf Interpretation.} SNN draws inspiration from the popular $K$ Nearest Neighbour (KNN) algorithm, described in Section (ref). However, the key assumption underlying KNN is that there {\em do exist} $K$ rows that are close to identical to the $i$-th row, with respect to some pre-defined metric. However, it is not necessary that these $K$ rows exist {\em even} for a rank $1$ matrix. As a simple example, consider a matrix $\bM \in \Rb^{m \times n}$ where $M_{i \cdot} = [i, 2i, \dots, ni]$. By construction $\bM$ is rank $1$, but for any row, there does not exist any other row that is close to it in a mean squared sense; hence, it has no nearest neighbours.
The SNN algorithm overcomes this hurdle by first constructing $K$ “synthetic” neighbors of row $i$ from $\textsf{NR}(j)$, where the $k$-th synthetic neighboring row is formed by a linear combination, defined by $\widehat{\beta}^{(k)}$, of the rows in $\textsf{AR}^{(k)}$. Then, similar to KNN, SNN estimates $A_{ij}$ by taking an average of the observed outcomes for column $j$ that are associated with the $K$ synthetic neighbors of row $i$. In words, $\widehat{\beta}^{(k)}$ is precisely the set of estimated linear weights that best recreates the observed outcomes of row $i$ from the rows in $\textsf{AR}^{(k)}$, using observations from the columns in $\textsf{AC}^{(k)}$. This idea of matching rows via a linear re-weighting takes inspiration from the synthetic controls literature---see Section (ref) for details. To ensure the linear fit is appropriately regularized, a spectral sparsity constraint is imposed on $\bS^{(k)}$, which is parameterized by $\lambda^{(k)}$. This constrained regression is known in the literature as principal component regression (PCR) (see agarwal2019robustness, agarwal2020principal, agarwal2021synthetic, agarwal2021causal). We note that in lieu of requiring that there exist $K$ close neighbouring rows as in KNN, SNN requires that the $i$-th row lies in the linear span of the rows in $\textsf{AR}^{(k)}$; that is, given Assumption (ref) holds, we require $|\textsf{AR}^{(k)}| \ge \mu$. Note that for the matrix $\bM$ described above, for any particular row, all other rows satisfy this linear span inclusion condition, i.e., $\mu = 1$ since $\bM$ is rank $1$.
{\bf Choosing $\lambda^{(k)}$} There exist a number of principled heuristics to select the hyper-parameter $\lambda^{(k)}$, and we name a few here. As is standard within the statistics and ML literatures, the most popular data-driven approach is to use cross-validation. Another common approach is to use a universal thresholding scheme that preserves the singular values above a precomputed threshold (see donoho14, Chatterjee15). Finally, a human-in-the-loop approach is to inspect the spectral characteristics of $\bS^{(k)}$ and choose $\lambda^{(k)}$ to be the natural “elbow” point that partitions the singular values into those of large and small magnitudes; in such a setting, the large magnitude singular values, which typically correspond to signal, are retained while the small magnitude singular values, which are often induced by noise, are filtered out. See the exposition on choosing the hyper-parameter for PCR in agarwal2019robustness, agarwal2020principal, agarwal2021synthetic.
{\bf Another Perspective on SNN.} SNN imputes $A_{ij}$ by building synthetic neighbors of row $i$ from $\textsf{NR}(j)$. In Proposition (ref), we demonstrate that $A_{ij}$ can be equivalently estimated by building synthetic neighbors of column $j$ from $\textsf{NC}(i)$ through a simple “transposition” of Algorithm (ref).
Note the SNN algorithm takes as input $\{(\textsf{AC}^{(k)}, \textsf{AR}^{(k)}): k \in [K]\}$. However, the question remains that given the matrix $\bD$, how to find these anchor rows and columns, with the additional constraint that the $K$ set of anchor rows $\{\textsf{AR}^{(k)}: k \in [K]\}$ are mutually joint. In Section (ref), we provide a practical algorithm $\texttt{AnchorSubMatrix}$ in Algorithm (ref) to find $\{(\textsf{AC}^{(k)}, \textsf{AR}^{(k)}): k \in [K]\}$ for a given pair $(i,j)$. In Section (ref), we discuss some motivating applications where anchor rows and columns are naturally induced.
In particular, we reduce our task of finding anchor rows and columns to a well-known problem in the graph theory literature known as finding “maximum bicliques”. We briefly explain how to do this simple reduction. We first introduce some standard notation from graph theory. Let $\mathcal{G} = (\mathcal{V}_1, \mathcal{V}_2, \Ec)$ denote a bipartite graph, where $(\mathcal{V}_1, \mathcal{V}_2)$ are the disjoint vertex sets and $\Ec \in \mathcal{V}_1 \times \mathcal{V}_2$ is the edge set, i.e., $(v_1, v_2) \in \Ec$ if there an edge between $v_1$ and $v_2$. Another way of representing $\mathcal{G}$ is via a bipartite incidence matrix $\bB \in \{0, 1\}^{|\mathcal{V}_1| \times |\mathcal{V}_2|}$ (or adjacency matrix). In particular, $B_{ij} = 1$ if $(v_i, v_j) \in \Ec$. If a sub-graph of $\mathcal{G}$ is complete, also called a biclique, then we denote it as $\BCc \subset \mathcal{G}$, i.e., there is an edge between any pair of nodes $(v_1, v_2) \in \BCc$. Now to see how to do the reduction between finding anchor rows and columns to the maximum biclique problem, recall $\bD \in \{0, 1\}^{m \times n}$ is our matrix of intervention assignments. Note, $\bD$ immediately induces a bipartite graph with $|\mathcal{V}_1| = m$ and $|\mathcal{V}_2| = n$, i.e., the vertex sets $\mathcal{V}_1$ and $\mathcal{V}_2$ correspond to the rows and columns of $\bD$, respectively. We define $\Ec$ as follows, $(v_i, v_j) \in \Ec$ if $D_{ij} = 1$. In other words, the incidence matrix $\bB \in \{0, 1\}^{m \times n}$ induced by this graph is exactly equal to $\bD$, i.e., $B_{ij} = 1$ if and only if $D_{ij} = 1$.
Given this reduction, we now describe how to practically implement the AnchorSubMatrix algorithm. We assume access to two algorithms: createGraph and maxBiclique. The former, $\texttt{createGraph}: \bB \rightarrow \mathcal{G}$, takes as input a bipartite incidence matrix $\bB$ (or adjacency matrix) and returns a bipartite graph $\mathcal{G}$; we note that the Python package NetworkX is an excellent resource to generate such graphs. The latter, $\texttt{maxBiclique}: \mathcal{G} \rightarrow \{\BCc^{(\ell)}\}_{\ell \in [L]}$, takes as input a bipartite graph $\mathcal{G}$ and returns a set of $L$ maximal bicliques $\{\BCc^{(\ell)}\}_{\ell \in [L]}$; we refer the interested reader to biclique, biclique1, biclique2, biclique3 and references therein for example algorithms.
Given $(\textsf{AC}, \textsf{AR})$ from Algorithm (ref), we can construct $\{(\textsf{AC}^{(k)}, \textsf{AR}^{(k)}): k \in [K]\}$ as follows: First, we assign $\textsf{AC}^{(k)} \leftarrow \textsf{AC}$ for every $k$, i.e., the anchor columns for each subgroup $k$ are all identically equal to $\textsf{AC}$. Second, we (randomly) partition $\textsf{AR}$ into $K$ subgroups of equal size and then assign $\textsf{AR}^{(k)}$ as the $k$-th subgroup of $\textsf{AR}$ such that $|\textsf{AR}^{(k)}| \sim |\textsf{AR}|/K$; in doing so, we ensure that $\{\textsf{AR}^{(k)}: k \in [K]\}$ are mutually disjoint sets. Note, for the purposes of theoretical analysis, we do not necessarily need to have $\textsf{AC}^{(k)}$ be identical across all $K$. In Section (ref), we show how the estimation error of SNN scales with $|\textsf{AR}^{(k)}|$ and $|\textsf{AC}^{(k)}|$. In short, our theoretical results suggest that we want $\{(\textsf{AC}^{(k)}, \textsf{AR}^{(k)}): k \in [K]\}$ to be large on average; It is sufficient that we choose $|\textsf{AC}^{(k)}|, |\textsf{AR}^{(k)}|$ such that $\min_{k \in [K]}\{|\textsf{AC}^{(k)}|, |\textsf{AR}^{(k)}| \}$ is as large as possible; this is essence what Step 5 of Algorithm (ref) is doing.
In this section, we discuss the typical sparsity pattern in recommender systems and sequential decision-making paradigms, which include panel data settings, reinforcement learning, and sequential A/B testing. We argue why these applications have a sparsity pattern where anchor rows and columns are naturally induced.
{\bf Recommender systems.} As stated earlier, one of the key motivating applications for matrix completion is recommender systems. It has been noted in ma2019missing that real-world recommender systems exhibit block-sparse structure; further the sparsity pattern is such that there is dependent missingness (i.e., $D_{ij} \not\perp\!\!\!\perp D_{ab}$ and zero probability of observing certain entries (i.e., $p_\min = 0$). An extreme version of this selection-bias would induce a sparsity structure as shown in Figure (ref). Within the context of movie recommender systems, a narrative for this missingness pattern is one where users only watch films that belong to genre(s) that they like and nothing else. However, in many recommender system applications, there exists a dense sub-matrix which corresponds to items that all users commonly rate---this corresponds to the rightmost columns of Figure (ref). This could occur if say a platform ask new users to indicate a subset of films that they enjoy. Indeed, this is a common practice for online platforms such as {\em Hulu, Netflix, StitchFix} to quickly learn a new user's preferences in order to provide a “warm-start” to their recommendation engine. Alternatively, many a time there are a small subset of iconic films (e.g., {\em Titanic} or {\em Star Wars}) that a large majority of users have watched. In this example in Figure (ref), all users can be used as anchor rows, and the set of items that are commonly rated across all users can be used as anchor columns. Further, we remark that in this example $\bP$ is not low-rank, thus violating the key assumption required to learn $p_{ij}$ in ma2019missing, cai2020uncertainty, bhattacharya2021matrix.
{\bf Sequential decision-making.} As described earlier, in sequential decision-making, data is collected across units (e.g., individuals, customer types, geographic locations) over time in a sequential manner, where each unit is likely to be observed under a single or small set of interventions out of many at any time period. Many sequential decision-making problems can be phrased this way, including (i) panel data settings in econometrics; (ii) reinforcement learning and its variants (e.g., online learning, contextual bandits); (iii) sequential A/B testing. In (ii), an intervention denotes both the action picked and the observed state for that given time period; meanwhile in (iii), platforms run experiments on different customer types in a sequential and/or adaptive manner over time. The induced matrix in these settings has rows index units and columns index time-intervention pairs. It is common in many of these sequential decision-making settings that there is a time period when all units are under the same intervention. This is usually done to collect “control” data about each unit to establish its baseline. For example, in an e-commerce setting, companies commonly estimate the baseline engagement level of a customer to understand the treatment effect of a discount policy; similarly, in clinical trials, pharmaceutical companies collect health metrics of patients to establish the treatment effect a particular therapy has. Further, the assumption that such a control period exists is standard in the synthetic controls literature. For an illustration of the sparsity pattern in the induced matrix with a control period, see Figure (ref). Hence, this “control” period in sequential decision-making can serve as our anchor columns, and all units can serve as anchor rows.
Below, we establish the statistical properties of the SNN algorithm. Without loss of generality, we consider a specific pair $(i,j)$. Recall from our discussion earlier, we suppress dependencies on $(i,j)$, e.g., all anchor rows and columns $\textsf{AR}^{(k)}, \textsf{AC}^{(k)}$ are defined with respect to $(i,j)$. In Section (ref), we state additional assumptions required to establish the theoretical results. In Sections (ref) and (ref), we establish finite-sample consistency and asymptotic normality of the SNN algorithm for a given entry $(i, j)$. In Sections (ref) and (ref), we discuss our assumptions and theoretical results, respectively.
{\bf Notation.} For every vector $v \in \Rb^a$, let $\| v \|_p$ denotes its $\ell_p$-norm. For the remainder of this work, let $\Ec = \{\bU, \bV, \bD\}$, i.e., the collection of latent factors and the observed missingness pattern. Recall the definition of $\bS^{(k)}$, $\textsf{AR}^{(k)}$ and $\textsf{AC}^{(k)}$ from Section (ref). Moreover, for every $k \in [K]$, we denote the SVD of $\mathbb{E}[\bS^{(k)} ~|~ \Ec]$ as
here, $r^{(k)} = \text{rank}(\mathbb{E}[\bS^{(k)} ~|~ \Ec])$. We denote $\bU^{(k)} \in \Rb^{|\textsf{AR}^{(k)}| \times r^{(k)}}$ and $\bV^{(k)} \in \Rb^{|\textsf{AC}^{(k)}| \times r^{(k)}}$ as the matrices of left and right singular vectors, respectively, i.e., $u^{(k)}_\ell \in \Rb^{|\textsf{AR}^{(k)}|}$ and $v^{(k)}_\ell \in \Rb^{|\textsf{AC}^{(k)}|}$ form the $\ell$-th columns of $\bU^{(k)}$ and $\bV^{(k)}$, respectively.
We state additional assumptions required to establish guarantees for the SNN algorithm. In Section (ref) we provide interpretations for Assumptions (ref) and (ref); Assumptions (ref) and (ref) are relatively standard and self-explanatory. Below, $k$ is indexed over $[K]$, where recall $K$ is a hyper-parameter of the SNN algorithm.
The following result establishes that the SNN algorithm outputs entry-wise consistent estimates of $\bA$, i.e., we establish consistency in $\| \cdot \|_\max$-norm. To simplify notation, we will henceforth absorb dependencies on $\sigma$ into the constant within $O_p(\cdot)$. That is, we assume there exists an absolute constant $C \ge 0$ such that $\sigma \le C$.
Note, Theorem (ref) does not require $N \rightarrow \infty$ to establish consistency of the SNN estimator. Rather, that $|\textsf{AR}^{(k)}|, |\textsf{AC}^{(k)}|$ is growing {\em on average} (ignoring logarithmic factors and dependence on $\beta^{(k)}, r^{(k)}, \sigma$). However, we state Corollary (ref) to help further interpret our results in Section (ref).
{\bf Implication for matrix completion with MCAR data.} Proposition (ref) below shows that SNN provides uniform entry-wise consistency for matrix completion with MCAR data as a special case if $p$, the probability of observing an entry, is sufficiently large.
The following establishes that the entry-wise estimate $\hA_{ij}$ of the SNN algorithm is asymptotically normal around the target causal parameter $A_{ij}$.
{\bf Interpretation of Assumption (ref).} Assumption (ref) requires that the nonzero singular values of $\mathbb{E}[\bS^{(k)} ~|~ \Ec]$ are well-balanced. Such an assumption is quite standard with the econometrics factor model and matrix completion literature. For example, it is analogous to incoherence-style conditions; see Assumption A of bai2019matrix and the discussion of theoretical results in agarwal2019robustness. It is also closely related to the notion of pervasiveness, see Proposition 3.2 of fan2018eigenvector. Indeed, the assumption that there is a gap between the top few singular values of a matrix of interest, and the remaining singular values has been widely adopted in the econometrics literature of large dimensional factor analysis dating back to chamberlainfactor. Crucially though, these works within econometrics (e.g. bai2019matrix, fan2018eigenvector, chamberlainfactor) aim to accurately estimate the factors themselves, which require making additional assumptions about the spectra of the matrix of interest to ensure these factors are uniquely identifiable. Instead we simply require that these low-rank factors exist, but do not explicitly require accurately estimating them. Assumption (ref) has also been shown to hold with high-probability for the canonical probabilistic generating process used to analyze probabilistic principal component analysis in bayesianpca and probpca; here, the observations are assumed to be a high-dimensional embedding of a low-rank matrix with independent sub-Gaussian entries (see Proposition 4.2 of agarwal2019robustness). Within the matrix/tensor completion literature, for an overview of where the well-balanced spectra assumption is utilized, see cai2021nonconvex and references therein. Practically speaking, Assumption (ref) can be empirically validated by plotting the spectrum of $\bS^{(k)}$, defined in Algorithm (ref); if there is a natural “elbow” point in the singular spectrum of $\bS^{(k)}$, i.e., there are a relatively small number of singular values that have a large and approximately equal magnitude, and the remaining singular values are significantly smaller, then Assumption (ref) is likely to hold. For further discussion of this empirical robustness check, please refer to the related discussion in agarwal2021synthetic.
{\bf Interpretation of Assumption (ref).} Recall from Algorithm (ref) that we learn the model $\widehat{\beta}^{(k)}$ by regressing $q^{(k)}$ on $\bS^{(k)}$. $\hA_{ij}^{(k)}$ is then estimated by applying the model $\widehat{\beta}^{(k)}$ on the outcomes in $x^{(k)}$ (i.e., the entries in the $j$-th column of the rows $\textsf{AR}^{(k)}$). The key question that remains is why would a model learned between $q^{(k)}$ and $\bS^{(k)}$, {\em generalize} well to accurately estimate $A_{ij}^{(k)}$ using $\langle x^{(k)}, \widehat{\beta}^{(k)} \rangle$. Normally, in statistical learning, such generalization requires making distributional assumptions about the training data (i.e., $\bS^{(k)}$) and the testing data (i.e., $x^{(k)}$). For example, each column of $\bS^{(k)}$ and $x^{(k)}$ are sampled i.i.d. However, we do not want to make such an assumption as it is unrealistic in setting such as recommendation systems, e.g., the ratings users give different movies is likely to be neither identically nor independently distributed. Indeed, by conditioning on $\Ec$, we are implicitly conditioning on $\bU$ and $\bV$, which requires our analysis to be {\em instance dependent}, i.e., has to hold for the specific sampling of the latent factors $\bU$ and $\bV$. To circumvent making any distribution assumptions, we make the natural assumption that in expectation, $x^{(k)}$ lies within the linear span of $\bS^{(k)}$. Such a condition is necessary as well for generalization, e.g., if every entry of $\mathbb{E}[\bS^{(k)} ~|~ \Ec]$ is equal to $0$, then no meaningful model $\widehat{\beta}^{(k)}$ can learned. Such an assumption has also been explored in agarwal2020principal, agarwal2021synthetic, agarwal2021causal. In particular, in agarwal2021synthetic the authors provide a data-driven hypothesis test to verify when such a condition holds.
To ease the discussion of the interpretation of the results, we will ignore dependence on logarithmic factors, and $\beta^{(k)}, r^{(k)}, \sigma$.
{\bf Sample complexity.} Note that even if $D_{ij} = 1$, estimating $A_{ij}$ is not straightforward; we never get to observe $A_{ij}$, rather we only observe $Y_{ij}$, where $Y_{ij} = A_{ij} + \varepsilon_{ij}$. That is, even if $D_{ij} = 1$, our observation of $A_{ij}$ is corrupted by noise and we only get a single sample of it. Remarkably, despite having access to (at most) a single noisy sample of $A_{ij}$, the estimate $\hA_{ij}$ produced by the SNN algorithm is consistent and asymptotically normal around $A_{ij}$. Of course, this is assuming a low-rank factor model and a suitable observation pattern. Hypothetically, if we get $K$ independent noisy samples of $A_{ij}$, denoted by $Y^{(1)}_{ij}, \dots, Y^{(K)}_{ij}$, the maximum likelihood estimator would be the empirical mean, $\frac{1}{K} \sum^K_{k = 1} Y^{(k)}_{ij}$. In this hypothetical scenario, this empirical mean would concentrate around $A_{ij}$ with error scaling as $O(K^{-1/2})$, i.e., with this estimation procedure to obtain an additive error of $O(\delta)$, we would need $K = \Omega(\delta^{-2})$ independent copies.
Now in comparison to the hypothetical scenario above where we have access to $K$ independent samples, Corollary (ref) effectively establishes that with access to at most $N^2 \times K$ observations, the error of the SNN estimator scales as $O(\max(N^{-1/4}, K^{-1/2}))$; this is assuming $|\textsf{AC}^{(k)}|, |\textsf{AR}^{(k)}| = N$ for all $k \in [K]$, as in Corollary (ref). This implies for any $\delta > 0$, $A_{ij}$ can be estimated to within an additive error of $O(\delta)$, if $N = \Omega(\delta^{-1/4})$ and $K = \Omega(\delta^{-2})$. Hence, compared to if we had $K$ independent noisy samples of each $A_{ij}$, we pay an additional cost of $N^2$ in terms of the number of samples needed, and $N^{-1/4}$ in terms of the estimation error rate even though we either do not observe a sample of $A_{ij}$ (i.e., it is missing), or only observe a single, noisy instantiation of it in $\widetilde{Y}_{ij}$. That is, to obtain estimation error of $O(\delta)$, it requires $O(\delta^{-2} \times \delta^{-4})$ observations (across different entries).
Further, in the hypothetical scenario where we get $K$ independent noisy copies for each $(i, j)$, if we wanted to estimate $A_{ij}$ to within error $O(\delta)$ for all $(i, j)$, this would require $m \times n \times K$ observations, with $K = \Omega(\delta^{-2})$. In contrast, for SNN, if we assume that for all $(i, j)$, we can use the same set of anchor rows and columns, i.e., $\{|\textsf{AC}^{(k)}|, |\textsf{AR}^{(k)}|\}_{k \in [K]}$ can be chosen to be the same for all $(i, j)$, then one can easily verify that the number of observations we need to recover each $A_{ij}$ to within error $O(\delta)$ is at most $N^2 \times K + m \times N + n \times (N \times K)$,\footnote{Technically, we only need $N^2 \times K + (m - N - K) \times N + (n - N) \times (N \times K)$.} with $N = \Omega(\delta^{-1/4})$ and $K = \Omega(\delta^{-2})$. See Figure (ref) for a visual depiction of the observation pattern for which this holds. Thus, for any fixed $\delta > 0$, we can recover every entry $A_{ij}$ to within additive error $O(\delta)$, with access to only $O(m + n)$ observations, rather than $O(m \times n)$ observations as would be naively required.
{\bf Connections to causal transportability, transfer learning, learning with distribution shift.} We note that this problem of generalizing well without making an i.i.d assumption is known by a variety of terms across many fields of study; these include “causal transportability”, “transfer learning”, “learning with distribution shift”. Given that subspace inclusion, i.e., Assumption (ref) holds, we show that generalization is possible without making any distributional assumptions about the underlying signal matrix $\bA$. Indeed, our theoretical results in Theorem (ref) and (ref) can be interpreted as point-wise out-of-sample generalization error bounds, which are distribution free (i.e., instance dependent). This might be of independent interest.
The objective of this section is to compare the imputation accuracy of SNN against the state-of-the-art matrix completion algorithms for MNAR data. We describe these algorithms in Section (ref). We do two case studies. In Section (ref), we apply these various algorithms in the setting of recommender systems with different missingness patterns. In Section (ref), we do the same but using data from a classic panel data case study in the econometrics literature called “California Prop 99” abadie2.
In particular, we compare two types of algorithms for matrix completion against SNN; we choose these benchmarks to be in line with those considered in ma2019missing. The first group of algorithms does not account for entries being MNAR; these include PMF (pmf), SVD (funk2006netflix), SVD++ (koren08), softImpute (softimpute14), and KNN (LeeLiShahSong16); we remark that the algorithm proposed in athey2021matrix is similar to \texttt{softImpute} with the addition of separate fixed effects terms. In particular, both the algorithm design and associated analysis of these algorithms is for MCAR data. In contrast, the second group does account for the limited MNAR setting as described in Section (ref); these include \texttt{MaxNorm} (cai2016matrix), \texttt{ExpoMF} (Liang2016), and \texttt{WTN} (WTN).
With the exception of ExpoMF, we further consider IPW-variants of the other benchmark algorithms, i.e., for each algorithm, we first de-bias the loss function given in (ref) via propensity scores. We do not do so with ExpoMF as their algorithm does not lend itself to be de-based via propensity scores in a straightforward manner (also see ma2019missing). The propensity scores are estimated in two ways, which are in line with the MAR and limited MNAR setting described in Section (ref). (i) MAR setting: We provide meaningful additional covariates $(X_i, \tilde{X}_j)$ for row $i$ and column $j$ and use logistic regression to learn $\widehat{p}_{ij}$; if a matrix completion algorithm is de-biased in this way, we add LR in front of it, e.g., LR-PMF means the PMF algorithm is used to estimate $\bhA$ and the loss function is de-biased using logistic regression. (ii) Limited MNAR setting: We do not provide additional covariates and directly estimate $\widehat{p}_{ij}$ using the observed mask matrix $\bD$; this is done using the 1bitMC algorithm in davenport20141bit and algorithms de-biased in this manner has a prefix of \texttt{1bitMC} added to them.; this approach to de-bias MNAR data is in line with what is proposed in ma2019missing, yang2021tenips, bhattacharya2021matrix.
We consider two error metrics, root mean-squared-error (RMSE) and mean-absolute-error (MAE). For all benchmark algorithms, we use 5-fold cross validation to tune their hyper-parameters through grid search for every error metric, i.e., for each benchmark algorithm, we find its best performing hyper-parameters with respect to RMSE and MAE on the validation set and report the error metric-specific hyper-parameters on the test set for each error metric. For SNN, we choose $K=1$ and $\lambda^{(1)}$ as per donoho14, i.e., we do not tune the hyper-parameters of SNN nor do we optimize it for each error metric. We emphasize that only the algorithms with the LR prefix use the additional row and column covariates $X_i, \tilde{X}_j$.
We begin with recommendation systems, which is arguably the canonical matrix completion application. Through the recommendation systems setting, we present two MNAR missingness patterns---one obeys the standard assumptions on MNAR in the literature (which we often refer to as limited MNAR) while the other considers a more general MNAR setting. To better understand the effect of the underlying mechanism which leads to missingness on each algorithm's ability to perform imputation, we consider the “noiseless” case, i.e., $\widetilde{Y}_{ij} = A_{ij}$ if $D_{ij} =1$ and $\widetilde{Y}_{ij} = \star$ otherwise. We study the effect of additional noise $\varepsilon_{ij}$ in the panel data setting in Section (ref).
In our first illustration, our observation pattern reflects the self-selection bias phenomena where most users tend to provide ratings if they particularly liked or disliked an item. However, they are much less inclined to provide a rating for an item that they are lukewarm about. Our simulated setup also consists of “core users” and “core movies”. We use core users to represent movie fanatics or critics, for instance, who provide explicit feedback for a significant number of films. In the setting of movie recommendations, we use core items to represent iconic movies such as {\em Star Wars} or {\em Titanic} that have influenced future films and popular culture, and are largely viewed by the general audience. These can also represent the subset of items that online platforms such as {\em Hulu, Netflix, Stichfix} display to new users when prompting for their preferences.
{\bf Experimental setup.} We consider $m = 80$ users and $n=80$ movies. We choose the dimension of the latent space as $r=5$. We generate the latent user matrix $\bU \in \Rb^{m \times r}$ as follows: (i) we first choose $m_\text{core} =20$ core users and construct $\bU_0 \in \Rb^{m_\text{core} \times r}$ by sampling entries i.i.d. from a standard normal distribution; (ii) next, we construct $\bU_1 = \bB \bU_0 \in \Rb^{(m-m_\text{core}) \times r}$, where the entries in $\bB \in \Rb^{(m-m_\text{core}) \times m_\text{core}}$ are sampled i.i.d. from a Dirichlet distribution, which ensures that the new factors lie in the same intervals as the factors in $\bU_0$. In doing so, every row of $\bU_1$, representing the latent factors corresponding to the “standard” users, is a linear combination of that of core users $\bU_0$, i.e., every standard user can be expressed as a weighted combination of core users. We then define $\bU = [\bU_0, \bU_1]$ such that the first $m_\text{core}$ rows of $\bU$ correspond to the core users. We construct $\bV = [\bV_0, \bV_1] \in \Rb^{n \times r}$ similarly, where $\bV_0 \in \Rb^{n_\text{core} \times r}$ and $\bV_1 \in \Rb^{(n - n_\text{core}) \times r}$ represent the matrix of latent factors associated with core movies and standard movies, respectively; here, we choose $n_\text{core} = 20$. We form $\bA = \bU \bV^T \in \Rb^{m \times n}$ and scale the values to lie within the interval $[1, 5]$; by construction, $\bA$ is a low-rank matrix.
Finally, we generate user and movie covariates matrices $\bX = \bU \bQ_1 \in \Rb^{m \times 3}$ and $\tilde{\bX} = \bV \bQ_2 \in \Rb^{n \times 3}$, where the entries in $\bQ_1 \in \Rb^{r \times 3}$ and $\bQ_2 \in \Rb^{r \times 3}$ sampled i.i.d. from a standard normal $\mathcal{N}(0,1)$; additionally, we normalize the columns in $\bQ_1$ and $\bQ_2$ to have unit $\ell_2$-norm.
Next, we describe our generative model for the propensity matrix $\bP \in \Rb^{m \times n}$. Without loss of generality, we denote $\mathcal{C}_\text{core} := \{(i,j): i \le m_\text{core}, j \le n_\text{core}\}$ as the subset of core users and core movies, $\mathcal{C}_\text{user} := \{(i,j): i \le m_\text{core}, j > n_\text{core}\}$ as the subset of core users and standard movies, $\mathcal{C}_\text{item} := \{(i,j): i > m_\text{core}, j \le n_\text{core}\}$ as the subset of standard users and core movies, and $\mathcal{C}_\text{standard} := \{(i,j): i > m_\text{core}, j > n_\text{core}\}$ as the subset of standard users and standard movies. These will represent our four cohorts of interest. Next, for some $t \in (1, 5)$, $\kappa_{ij} > 0$, and $\alpha_{ij} \in (0,1]$,
In our setting, we choose our threshold $t = 2.3$. Here, $\alpha_{ij}$ is a parameter that controls the MNAR effect: $\alpha_{ij} = 1$ is MCAR while $\alpha_{ij} \rightarrow 0$ only reveals $1$ and $5$ rated movies. We choose $\alpha_{ij} = 0.7$ for $(i,j) \in \mathcal{C}_\text{core}$, $\alpha_{ij} = 0.35$ for $(i,j) \in \mathcal{C}_\text{user}, \mathcal{C}_\text{item}$, and $\alpha_{ij}=0.1$ for $(i,j) \in \mathcal{C}_\text{standard}$. For every $(i,j)$ pair, $\kappa_{ij}$ is set so that the expected number of revealed ratings within the cohort is equal to some value. We choose the expected number of observations within $\mathcal{C}_\text{core}$ as 90%, within $\mathcal{C}_\text{user}$ as 70%, within $\mathcal{C}_\text{item}$ as 70%, and within $\mathcal{C}_\text{standard}$ as 5%. This sampling process ensures the two key assumptions in the limited MNAR setting of the entries of $\bD$ being independent and $p_\min > 0$ are satisfied. See Figure (ref) for a visual depiction of empirical sparsity pattern under this missingness mechanism.
{\bf Results.} In the following simulations, we obey the generative process above. In particular, we sample $\bA$ and $\bP$ once, as well as $\bX$ and $\tilde{\bX}$, and perform 10 experimental repeats where the only randomization lies in the sparsity pattern, i.e., we observe 10 independent realizations of $\bD$. We report the average RMSEs and MAEs, as well as their respective standard deviations, over the 10 experimental runs in Table (ref). We find that with respect to MAE, SNN achieves the best result along with MaxNorm (and its variants); with respect to RMSE, SNN is a close second with SVD++ (and its variants), after MaxNorm (and its variants). Although positivity and independence between entries in $\bD$ are upheld, we remark that debiasing via 1bitMC and \texttt{LR-} do not always yield stronger results, e.g., see \texttt{PMF} and \texttt{softImpute}.
In this simulation, we violate two key assumptions in the current literature on MNAR data: (i) positivity and (ii) independence between the entries in $\bD$. Towards this, we continue the notion of core movies, for which all users provide ratings. For the remaining movies, users only provide ratings if a movie belong to their favorite genre. This deterministically sets every entry in $\bP$ (and thus $\bD$) to either $0$ or $1$, and correlates the entries in $\bD$, which yields a sparsity pattern similar to that shown in Figure (ref). See Figure (ref) for a visual depiction of empirical sparsity pattern under this missingness mechanism.
{\bf Experimental setup.} In particular, we consider $m=80$ users and $n=80$ items. We choose the dimension of the latent space as $r=5$. We generate the latent user matrix $\bU \in \Rb^{m \times r}$ by sampling entries i.i.d. from a standard normal distribution. To generate the latent item matrix $\bV \in \Rb^{n \times r}$, we first choose $n_{\text{core}}=30$ core items (to be defined in greater detail below) and construct $\bV_0 \in \Rb^{n_{\text{core}} \times r}$ by sampling entries i.i.d. from a standard normal. Next, we construct $\bV_1 = \bB \bV_0 \in \Rb^{(n-n_{\text{core}}) \times r}$, where the entries in $\bB \in \Rb^{(n-n_{\text{core}}) \times n_{\text{core}}}$ are sampled i.i.d. from a Dirichlet distribution. In doing so, every row of $\bV_1$ is a linear combination of rows in $\bV_0$, i.e., every item can be expressed as a weighted combination of core items. We then define $\bV = [\bV_0, \bV_1]$ such that the first $n_{\text{core}}$ rows of $\bV$ correspond to the core items. We form $\bA = \bU \bV^T \in \Rb^{m \times n}$ and scale the values to lie within the interval $[1,5]$. Finally, we generate user and item feature matrices $\bX = \bU \bQ_1 \in \Rb^{m \times 10}$ and $\tilde{\bX} = \bV \bQ_2 \in \Rb^{n \times 10}$, where the entries in $\bQ_1 \in \Rb^{r \times 10}$ and $\bQ_2 \in \Rb^{r \times 10}$ sampled i.i.d. from a standard normal $\mathcal{N}(0,1)$; additionally, we normalize the columns in $\bQ_1$ and $\bQ_2$ to have unit $\ell_2$-norm. We generate higher dimensional covariates $\bX$ and $\tilde{\bX}$ to see if improves the relative performance of the MAR algorithms, denoted by the prefix LR, which use this additional information to estimate the propensities.
To describe the generating process for the observation pattern $\bD$, we begin by providing an interpretation of the above quantities. First, we interpret $r$ as the number of latent genres. In turn, the $(i,k)$-th entry in $\bU$ can be interpreted as user $i$'s preference for genre $k$; similarly, the $(j, k)$-th entry in $\bV$ can be interpreted as the level to which item $j$ is composed of genre $k$. We consider the setting where all users provide ratings for all core items, i.e., $D_{ij} = 1$ for every user $i \in [m]$ and core item $j \in [n_{\text{core}}]$. For the remaining entries in $\bD$, we posit that every user will only rate items from their favorite genre. More specifically, given the above interpretation, we define user $i$'s favorite genre $k^*(i)$ as $k^*(i) = \operatorname*{\arg\!\max}_{k \in [r]} U_{ik}$; similarly, we classify an item $j$ as belonging to genre $k^\sharp(j)$ if $k^\sharp(j) = \operatorname*{\arg\!\max}_{k \in [r]} V_{jk}$. Hence, for every user $i \in [m]$ and non-core item $j > n_{\text{core}}$, we have $D_{ij} = 1$ if $k^*(i) = k^\sharp(j)$ and $0$ otherwise. We underscore that this model violates the standard operating assumptions within the current MNAR literature as entries in $\bD$ are deterministically set to $0$ (i.e., the minimum element in $\bP$ is $0$), and are dependent on one another.
{\bf Results.} In the following simulations, we obey the generative process above. In particular, we sample $\bV$ and $\bX$ once, and perform 10 experimental repeats where the only randomization lies in the re-sampling of $\bU$; this is done to model new users coming into the system with the movies fixed. We report the average RMSEs and MAEs, as well as their respective standard deviations, over the 10 experimental runs in Table (ref). We find that SNN achieves the best RMSE and MAE, with 1bitMC-MaxNorm as a close second with respect to RMSE and MAE. As with the limited MNAR setting experiment, we find that de-biasing does not always improve results. This is reasonable given that our generative process violates the typical assumptions underlying propensity estimation methods. The relative improvement of SNN shows its robustness to the general MNAR setting, where entry-wise positivity and independence of $\bD$ are violated. The fact that $\texttt{KNN}$ performs relatively poorly indicates that matching via linear weights is indeed more expressive than matching with uniforms weights as done in $\texttt{KNN}$. WE also note that the various state-of-the-art algorithms are still relatively robust to the general MNAR setting. This may warrant further investigation into the potential gap between theory and practice on the robustness of these methods to different missingness patterns.
We now compare SNN against the same benchmark matrix completion algorithms using a classic case study of California smoking data of abadie2, which has been widely utilized within the econometrics literature. We do so as this setting has a MNAR sparsity pattern which is quite distinct from what is seen in recommendation systems. We now give a brief overview of the case study. In 1988, California introduced the first modern-time large-scale anti-tobacco legislation in the United States (Proposition 99). There was interest in estimating the effect of this legislation on tobacco sales in California. Towards this, per-capita cigarette sales data was collected across 39 U.S. states from 1970 to 2000. Among the 39 states, there was one “treated” state, California, which implemented the legislation; the remaining 38 states were chosen as “control” states as they neither instituted a tobacco control program nor raised cigarette sales taxes by 50 cents or more. These other 38 control states were then used to build a “synthetic California”, i.e., a synthetic trajectory of cigarette sales in California if it had not introduced any tobacco legislation.
{\bf Experimental setup.} We consider the time horizon of $n=31$ years and restrict our focus to the $m = 38$ control units in the original dataset. This data is encoded into a $38 \times 31$ matrix, $\boldsymbol{Y}$, where the entry $Y_{ij}$ represents the potential outcome of per-capita cigarette sales (in packs) for state $i$ in year $j$ under control, i.e., without any intervention in place. To generate MNAR data, we artificially introduce interventions to a subset of states in 1989, where the probability a state adopts an intervention (e.g., tobacco control program) depends on their change in cigarette sales pre- and post-1989. More specifically, we consider the following treatment adoption protocol: First, we cluster states into three categories---mild, moderate, or severe---based on their change in average cigarette sales during 1989-2000 compared to that during 1970-1988; we note that in this context, a negative change means that the cigarette sales in the post-intervention period are lower than that in the pre-intervention period. As such, we define (i) mild states as those whose change is at least one standard deviation above the average change across all states; (ii) severe states as those whose change is at least one standard deviation below the average change across all states; (iii) and moderate states as the remaining states whose change is within one standard deviation.
We then designate the probability of intervention for mild, moderate, and severe states as $10\%$, $30\%$, and $50\%$, respectively. In words, this setup reflects the scenario in which a state is more likely to adopt an intervention if their average sales in the post-intervention period is relatively closer to their pre-intervention sales compared to that of their peer states. In the language of causal inference, this is exactly confounding, i.e., there is a correlation between the treatment assignment and the eventual outcome.
For an example of a mild, moderate, and severe state, please see Figure (ref). Additionally, we remark that once an intervention is adopted, all sales under control during the post-intervention period are, by definition, unobserved, i.e., $\widetilde{Y}_{ij} = \star$ for any intervened on state $i$ and for all $j \ge 19$ (after 1988); this yields the observation pattern shown in Figure (ref). Finally, to employ logistic regression, i.e., LR to de-bias the estimates, we use state covariate data from abadie2, which include average retail price of cigarettes, per capita state personal income (logged), the percentage of the population age 15-24, and per capita beer consumption. We note that SNN does not use this auxiliary data.
{\bf Results.} Using the above setup, we apply the various matrix completion methods to impute the missing counterfactual cigarette sales associated with the artificial intervention states during the post-intervention period. We report the average root mean-squared-errors (RMSEs) and mean absolute errors (MAEs), as well as their respective standard deviations, over 10 experimental runs in Table (ref). As the table shows, SNN significantly outperforms all baseline algorithms under both error metrics. The only exception is KNN, which performs similarly to SNN; this is interesting as KNN is in essence, the difference-in-differences estimator, a standard method within the panel data econometrics literature. Further, SVD++ and MaxNorm (and its variants), which performed strongly in the recommendation systems example, now incur a significant error. We display a few representative results in Figure (ref). Collectively across all three studies, we find that \texttt{SNN} is robust under varying missingness mechanisms.