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.
69,737 characters · 10 sections · 29 citation commands
Robust Matrix Estimation with Side Information
\doparttoc \faketableofcontents
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 \fi
\spacingset{1.12}
{\it Keywords: Matrix completion, Nuclear norm penalization, Causal inference, Non-linear estimation, Covariate information}
\part
\spacingset{1.7}
Recent technological progress has made it possible to gather and process high-volume data that are conveniently organized in matrix form, often with both dimensions scaling up rapidly. Accordingly, high-dimensional matrix estimation problems such as matrix denoising and matrix completion have attracted considerable attention, and many impressive results have been obtained from both statistical and computational perspectives. However, although side information is often available in addition to the target outcome data, traditional approaches typically use only the outcome data for matrix estimation. Incorporating side information can enrich the underlying model and improve estimation and prediction accuracy. As our ability to access auxiliary covariate data continues to grow, developing matrix estimation methods that effectively leverage side information has become an important and timely research direction.
A number of computational algorithms, along with their statistical properties, have been proposed recently. Arguably, the most popular model that incorporates additional information in matrix estimation is the Inductive Matrix Completion (IMC) model xu2013speedup,jain2013provable,zhang2018fast. The standard IMC model takes the form: $$ Y = M = X L Z^\top $$ where $Y = [y_{it}]_{i \leq N,, t \leq T}$ is the outcome matrix, $X = [x_1, \ldots, x_N]^\top$ is the $N \times d_1$ row-feature matrix, $Z = [z_1, \ldots, z_T]^\top$ is the $T \times d_2$ column-feature matrix, and $M$ is the target matrix. Here, $L$ is assumed to be a low-rank $d_1 \times d_2$ matrix. A typical estimation approach is to solve $$ \min_{L \in \mathbb{R}^{d_1 \times d_2}} \left\Vert\mathcal{P}_{\Omega} (XLZ^\top - Y)\right\Vert_F^2 + \lambda \left\VertL\right\Vert_{*}, $$ where $||\cdot||_*$ denotes the nuclear norm, $\mathcal{P}_{\Omega}(A) = \Omega \circ A$, and $\Omega$ is the $N \times T$ indicator matrix for observability in matrix completion.
Although IMC is a popular and useful approach to matrix estimation with side information, it has several important limitations. First, it requires that the features be present on both sides and also interact linearly. Moreover, it predicates upon the informativeness of both row and column features and can break down if features are weak or irrelevant. Several extensions of IMC have been proposed in recent years to address these shortcomings. ledent2023generalization incorporate a noise component $E$ into the IMC model and derive bounds on the expected $\ell$-risk. zhong2019provable allow a nonlinear relationship between $(X, Z)$ and $M$. wang2018high consider settings in which the rank of $L$ can be large. Notably, chiang2015matrix propose the so-called “dirty” IMC model, which augments the standard IMC formulation with an additional low-rank term $R$ and estimates $(L, R)$ by solving $$ \min_{L ,R } \left\Vert\mathcal{P}_{\Omega} (XLZ^\top + R - Y)\right\Vert_F^2 + \lambda_1 \left\VertL\right\Vert_{*} + \lambda_2 \left\VertR\right\Vert_{*}. $$ However, this model does not include a noise term $E$ and components explained only by one-sided characteristics. In addition, it does not allow a nonlinear relationship between $(X, Z)$ and $M$, and it still requires $L$ to be low-rank. Overall, each extension addresses only part of the limitations and still leaves other issues unresolved.
Another notable line of research on matrix completion with covariates includes mao2019matrix and ma2025statistical. These papers consider the model $$ Y = X B^\top + R + E $$ where $B$ is an unknown coefficient matrix and $R$ is a low-rank matrix. Because this approach does not incorporate column characteristics $Z$, it cannot capture interaction terms involving both $X$ and $Z$ (such as $XLZ^\top$) or components explained solely by $Z$. In addition, it does not allow for nonlinear effects. As a result, the model still has some limitations.
Lastly, a related strand of work studies PCA or factor analysis in settings without missing data fan2016projected,chiang2016robust,niranjan2017provable,zhu2016personalized,xue2017side. For example, fan2016projected consider the model, $Y = (G(X) + \Gamma) F^\top + E = G(X) F^\top + \Gamma F^\top + E$, where $G(X)$ is a part of loading defined by an unknown function of $X$, and $F$ denotes unobserved factors. In contrast, zhu2016personalized study the model, $Y = XB^\top +AZ^\top + R + E$, where $R$ is low-rank. The former framework cannot capture interaction terms involving both $X$ and $Z$ (such as $XLZ^\top$) or components explained solely by $Z$, whereas the latter does not include an interaction term between $X$ and $Z$. Moreover, this model primarily emphasizes linear relationships.
To overcome the limitations of existing approaches, we consider the following model:
where $G_1, G_2, Q_1,$ and $Q_2$ are unknown matrix-valued functions, and $W_1, W_2, V_1,$ and $V_2$ are unobserved matrices. This model is more general and nests the above models. For example, the models in xu2013speedup,jain2013provable,wang2018high correspond to the special case $M = M_1$, and the model in chiang2015matrix corresponds to $M = M_1 + M_4$. In addition, the models in fan2016projected,mao2019matrix,ma2025statistical can be viewed as special cases of $M = M_2 + M_4$. For instance, the model in fan2016projected can be represented as $M = M_2 + M_4$ with $V_1 = V_2$. Hence, our estimation approach under this model is less likely to suffer from model misspecification. Moreover, if the data contain components that existing models do not account for, our estimator is expected to perform better than estimators based on those restricted models. As discussed in Section (ref), the convergence rates of our estimator demonstrate the robustness of our method across models, and the simulation results in Section (ref) are consistent with these theoretical findings.
Our estimation is based on a sieve projection method. Using projection matrices constructed from sieve bases for $X$ and $Z$, we estimate each component of $M$ separately and then obtain an estimator of $M$ by summing these estimates. Thanks to the sieve projection, our method can accommodate potentially nonlinear effects of $X$ and $Z$ on $M$. In addition, estimating each component separately allows us to fully exploit the model structure in (ref). Together, these features make our estimator more likely to outperform methods based on more restrictive models when the data contain components that those restrictive models do not account for.
Another important feature of our approach is the use of nuclear-norm penalization, which corresponds to a soft-thresholding procedure. Hence, if some of $M_2$, $M_3$, and $M_4$ are exactly zero, then our estimators for those components are also exactly zero with high probability. This property enhances the robustness of our estimator.
In contrast, if we use a spectral method to estimate each component, we must estimate the rank of each part, and existing rank estimators may produce incorrect (nonzero) estimates when the corresponding component is weak due to a low signal-to-noise ratio. As a result, spectral methods may perform poorly when some of $M_2$, $M_3$, and $M_4$ are zero or close to zero. By comparison, our nuclear-norm–penalized estimator does not require estimating the rank of each component or the signal strength of each component; therefore, small values of $M_2$, $M_3$, and $M_4$ do not pose a problem.
Another important contribution of this paper is that we also consider a setting with missing entries, where the missingness is not at random. While many papers use side information for imputation under MAR (missing at random), to the best of our knowledge, no existing matrix completion work incorporates side information under MNAR (missing not at random). Since the seminal work of athey2021matrix, which demonstrated that matrix completion techniques can be very useful for causal panel data models, matrix completion has become a popular tool for estimating unobserved potential outcomes under the untreated (control) condition. However, the potential-outcome matrix under the untreated condition usually exhibits a missingness pattern that does not follow random missingness. Consequently, matrix completion under MNAR—and its applications to causal inference—has been actively studied recently (see, e.g., athey2021matrix,bai2021matrix,agarwal2023causal,choi2024matrix,yan2024entrywise). We propose a novel matrix completion method that leverages side information under MNAR. As shown in our real-data experiment in Section (ref), our method outperforms standard matrix-completion approaches in imputing unobserved potential outcomes and demonstrates its usefulness for treatment-effect estimation.
The remainder of this paper is organized as follows. Section (ref) introduces our model and our estimation method, which uses sieve projection with nuclear norm penalization. Section (ref) presents asymptotic error bounds for the estimator and discusses the robustness of our method across different models. Section (ref) extends our estimation strategy to the case in which the outcome matrix is partially observed. Importantly, we consider the MNAR setting as well as the MAR setting. Section (ref) presents numerical studies using simulated and real data to demonstrate the advantages of our method. All proofs are relegated to the supplement due to space limitations.
In this paper, we consider the following panel model: $$ Y = M + E, \qquad M = M_1 + M_2 + M_3 + M_4, $$ where \(Y = (y_{it})_{i\le N,\, t\le T}\) is the outcome matrix, \(E = (\epsilon_{it})_{i\le N,\, t\le T}\) is the noise matrix, and \(M = (m_{it})_{i\le N,\, t\le T}\) is the matrix of interest. We decompose \(M\) into four parts: (i) \(M_1\), a component well explained by both \(X\) and \(Z\); (ii) \(M_2\), a component explained by \(X\) but irrelevant to \(Z\); (iii) \(M_3\), a component explained by \(Z\) but irrelevant to \(X\); and (iv) \(M_4\), a component irrelevant to both \(X\) and \(Z\), where \(X = (x_i)_{i\le N}\) and \(Z = (z_t)_{t\le T}\) are observable characteristics corresponding to the row and column indices, respectively.
More specifically, each part can be represented as
where \(G_1(X) = (g_{1,k}(x_i))_{i\le N,\, k\le K_1}\), \(G_2(X) = (g_{2,k}(x_i))_{i\le N,\, k\le K_2}\), \(Q_1(Z) = (q_{1,k}(z_t))_{t\le T,\, k\le K_1}\), and \(Q_2(Z) = (q_{2,k}(z_t))_{t\le T,\, k\le K_3}\) for some unknown functions \(g_{1,k}(\cdot)\), \(g_{2,k}(\cdot)\), \(q_{1,k}(\cdot)\), and \(q_{2,k}(\cdot)\). Here, \(W_1 = (w_{1,ik})_{i\le N,\, k\le K_3}\) and \(W_2 = (w_{2,ik})_{i\le N,\, k\le K_4}\) capture the components not explained by \(X\), while \(V_1 = (v_{1,tk})_{t\le T,\, k\le K_2}\) and \(V_2 = (v_{2,tk})_{t\le T,\, k\le K_4}\) capture the components not explained by \(Z\). This model is general and encompasses many existing models.
\paragraph{Estimation.}
To properly accommodate and exploit the structure of our model in (ref), we propose estimating \(M\) using a sieve projection method. For two sets of basis functions \(\{\phi_1(x), \ldots, \phi_{J}(x)\}\) and \(\{\psi_1(z), \ldots, \psi_{J}(z)\}\) (e.g., B-splines, Fourier series, wavelets, or polynomial series), define
where \(d_1\) and \(d_2\) are the dimensions of \(x_i\) and \(z_t\), respectively. The corresponding projection matrices are \[ P_X = \Phi(X)\big(\Phi(X)^{\top}\Phi(X)\big)^{-1}\Phi(X)^{\top}, \qquad P_Z = \Psi(Z)\big(\Psi(Z)^{\top}\Psi(Z)\big)^{-1}\Psi(Z)^{\top}, \] where \(\Phi(X) = \big[ \boldsymbol{\phi}(x_1), \ldots, \boldsymbol{\phi}(x_N) \big]^\top\) and \(\Psi(Z) = \big[ \boldsymbol{\psi}(z_1), \ldots, \boldsymbol{\psi}(z_T) \big]^\top\). Note that as long as \(g_{1,k}(\cdot)\), \(g_{2,k}(\cdot)\), \(q_{1,k}(\cdot)\), and \(q_{2,k}(\cdot)\) are sufficiently smooth, for any \(\iota \in (1,2)\) we have
Moreover, \(\|P_X E P_Z\|_F\) can be much smaller than \(\|E\|_F\) due to the orthogonality between \((X,Z)\) and \(E\). Leveraging this property, we propose estimating \(M\) as follows.
To understand how this estimator works, note that by (ref) and basic properties of nuclear-norm penalization, we have
under suitable conditions on the noise and on the smoothness of \(g_{1,k}(\cdot)\), \(g_{2,k}(\cdot)\), \(q_{1,k}(\cdot)\), and \(q_{2,k}(\cdot)\). Importantly, the terms involving \(P_X W_\iota\) or \(P_Z V_\iota\) cancel out when we sum the four estimators. Hence, without imposing any orthogonality conditions between \(X\) and \(W\) (or between \(Z\) and \(V\)), our final estimator \(\widehat{M}\) can estimate \(M\) well.
In addition, because nuclear-norm penalization acts as a thresholding estimator, when some (or all) of \(M_2\), \(M_3\), and \(M_4\) are zero or sufficiently small, it helps us obtain a tighter bound.
In this section, we present the convergence rate of our estimator. We begin by imposing the following conditions.
The independence and sub-Gaussianity assumptions are used to derive tight bounds for \(\|P_X E P_Z\|\), \(\|P_X E\|\), and \(\|E P_Z\|\). We can generalize this condition to weakly dependent noise with suitable moment conditions, at the cost of additional \(J\)-dependent terms in the bound in Theorem (ref).
This condition is standard in the sieve-estimation literature (e.g., fan2016projected,chen2023semiparametric). Because we focus on the case where \(J d_1 \ll N\) and \(J d_2 \ll T\), it follows from the law of large numbers and therefore is not overly restrictive.
Condition (i) is a standard assumption in sieve estimation. For example, if \(g_{\iota,k}(\cdot)\) has an additive form $g_{\iota,k}(x_i) = \sum_{l=1}^{d_1} g_{\iota,kl}(x_{il})$ and each \(g_{\iota,kl}(\cdot)\) belongs to the H\"older class \(\mathcal{H}(\rho_\iota^G, \tau_\iota^G)\), where \[ \mathcal{H}(\rho,\tau) = \left\{ h : \big| h^{(\rho)}(s) - h^{(\rho)}(t) \big| \leq C |s-t|^\tau \right\} \] for some \(C>0\), then \(\gamma_\iota^G = \rho_\iota^G + \tau_\iota^G\) for typical choices of basis functions (see, e.g., chen2007large). On the other hand, condition (ii) requires sufficient smoothness of the functions \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\). Note that \(\gamma_\iota^G\) and \(\gamma_\iota^Q\) can be viewed as smoothness parameters for \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\), respectively. Therefore, if \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth, then \(\gamma_\iota^G\) and \(\gamma_\iota^Q\) will be large, and condition (ii) can be satisfied even when \(J\) increases slowly.
Lastly, we impose the following moment conditions.
We are now in a position to state the statistical properties of our estimators. The following theorem provides the convergence rate of the proposed estimator.
Some immediate remarks are in order. First, note that the dominating part of the error bound for our estimator does not depend on \(K_1\). Thus, we can allow \(K_1\) to be large as long as \(g_{1,k}(\cdot)\) and \(q_{1,k}(\cdot)\) are sufficiently smooth. The last term, \[ \sqrt{NT}\left[ \frac{K_1}{J^{\gamma_1^G}} + \frac{K_1}{J^{\gamma_1^Q}} + \frac{K_2}{J^{\gamma_2^G}} + \frac{K_3}{J^{\gamma_2^Q}} \right], \] arises from the sieve approximation (smoothing) error. When the functions \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth (i.e., \(\gamma_\iota^G\) and \(\gamma_\iota^Q\) are large), this term is small and dominated by the other terms.
Theorem (ref) illustrates the robustness of our estimator. First, consider the most favorable case, \(M=M_1\), where \(M\) is well explained by \(X\) and \(Z\). In this case, the convergence rate of our estimator is \(O_p(J)\) provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth. This matches the rate of the “double projection” estimator \(P_X Y P_Z\), which is the most suitable estimator when we know a priori that \(M=M_1\). By contrast, if we estimate \(M\) using a standard low-rank method such as nuclear-norm penalization in this setting, the convergence rate would be \(O_p(\sqrt{K_1(N+T)})\), which is much larger than ours. Moreover, even when \(M\) contains an additional component \(M_2+M_3+M_4\) beyond \(M_1\), as long as this component is small (in the sense that \(\|M_2+M_3+M_4\|_F \ll \sqrt{N+T}\)), our rate \(O_p\big(J + \sqrt{K^*}\,\|M_2+M_3+M_4\|_F\big)\) is smaller than that of the usual low-rank estimator \(O_p\big(\sqrt{K(N+T)}\big)\), where \(K^*=\max\{K_2,K_3,K_4\}\) and \(K\) is the rank of \(M\).
Next, consider the least favorable case, \(M=M_4\), where \(M\) is unrelated to \(X\) and \(Z\) and the side information is uninformative. In this case, the convergence rate of our estimator is \(O_p(\sqrt{K_4(N+T)})\), provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth and \(J \ll \sqrt{K_4(N+T)}\). Note that this rate coincides with that of standard low-rank estimators. Hence, even in the least favorable case, the error bound for our estimator is comparable to that of a typical low-rank method. In contrast, the “double projection” estimator \(P_X Y P_Z\) is inconsistent in this case. Moreover, if \(M\) contains an additional small component \(M_1+M_2+M_3\) beyond \(M_4\), the convergence rate of our estimator becomes \(O_p(\sqrt{K^*(N+T)})\), whereas that of a typical low-rank estimator remains \(O_p(\sqrt{K(N+T)})\). Thus, when \(K_1\) is large, our method can yield a tighter bound.
Lastly, consider the case \(M=M_2\). In this case, the convergence rate of our estimator is \(O_p(\sqrt{K_2 T})\), provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth and \(J \ll \sqrt{K_2 T}\). By comparison, the convergence rate of a typical low-rank estimator is \(O_p(\sqrt{K_2(N+T)})\). Hence, when \(N \gg T\), our method yields a tighter bound. In addition, if \(M\) contains an additional small component \(M_1+M_3+M_4\) beyond \(M_2\), the convergence rate of our estimator becomes \(O_p(\sqrt{(K_2+K_4)T})\), whereas that of a typical low-rank estimator is \(O_p(\sqrt{K(N+T)})\). Thus, when \(K_1\) is large or \(N \gg T\), our method can yield a better bound. A similar discussion applies to the case \(M=M_3\).
Table (ref) summarizes the convergence rates of the estimators in the cases discussed above. We can see that, in every case, the convergence rate of our estimator is at least as good as that of the other estimators, provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth.
Lastly, as a corollary, we present convergence rates for the estimated singular vectors, since the factors and loadings are often of interest to researchers (see, e.g., bai2008large). Let \(\widehat{U} \in \mathbb{R}^{N \times K}\) and \(\widehat{V} \in \mathbb{R}^{T \times K}\) be the left and right singular vectors of \(\widehat{M}\), respectively. Similarly, let \(U\) and \(V\) denote the left and right singular vectors of \(M\), respectively. For notational convenience, denote the upper bound on \(\|\widehat{M}-M\|_F\) in Theorem (ref) by
Then, we obtain the following convergence rates.
Next, we extend our estimation strategy to the case where the outcome matrix is only partially observed. The base model is the same as in Section (ref), and we additionally assume that researchers observe \(\Omega \circ Y\) instead of \(Y\), where \(\Omega = (\omega_{it})_{i \leq N,\, t \leq T} \in \{0,1\}^{N \times T}\).
In this section, we consider the case where outcome entries are missing at random. Specifically, we assume that \((\omega_{it})_{i \leq N,\, t \leq T}\) are i.i.d.\ Bernoulli random variables with mean \(p\), as is common in the matrix completion literature.
\paragraph{Estimation.}
Similarly to the fully observed case above, we use the projection method to exploit the structure of the model in (ref). However, when entries are missing, a key difficulty is that we cannot directly observe \(\Omega \circ \big(P_X Y (I_T - P_Z)\big)\) and \(\Omega \circ \big((I_N - P_X) Y P_Z\big)\) when we aim to estimate \(P_X M (I_T - P_Z)\) or \((I_N - P_X) M P_Z\) via nuclear-norm penalization. On the other hand, we can still estimate \(P_X M P_Z\) accurately using the projection estimator \(p^{-1} P_X (\Omega \circ Y) P_Z\).
Hence, in the presence of missing entries, we proceed as follows.
Note that, by the projection relation (ref) and the usual properties of nuclear-norm penalization in matrix completion, we have
under conditions on the noise and on the smoothness of \(g_{1,k}(\cdot)\), \(g_{2,k}(\cdot)\), \(q_{1,k}(\cdot)\), and \(q_{2,k}(\cdot)\) similar to those in the previous section. As above, the terms involving \(P_X W_\iota\) or \(P_Z V_\iota\) cancel out when we add the two estimators. Hence, our final estimator \(\widehat{M}\) can estimate \(M\) well. In particular, because we estimate the \(M_1\) component using the projection method rather than a low-rank estimator, our approach can have advantages when \(K_1\) is large or when \(M_1\) is large relative to \(M_2\), \(M_3\), and \(M_4\).
\paragraph{Asymptotic result.}
We now present the convergence rate of our estimator. We begin by introducing several additional assumptions.
In addition to Assumption (ref), we require a slightly different condition on the sieve approximation error than in the fully observed case.
The following theorem provides the convergence rate of our estimator.
Similar discussions to the fully observed case apply here. The error bound for our estimator does not depend on \(K_1\) as long as \(g_{1,k}(\cdot)\) and \(q_{1,k}(\cdot)\) are sufficiently smooth, because we estimate \(M_1\) using the projection method rather than a low-rank estimator. Hence, we can allow \(K_1\) to be large.
In addition, the estimator enjoys a robustness property. For simplicity, assume that \(M_{\max}\) and \(\|M\|_{\infty}\) are bounded. First, consider the most favorable case, \(M=M_1\). In this case, the convergence rate of our estimator is \(O_p(J/\sqrt{p})\) provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth. By contrast, if we estimate \(M\) using standard low-rank completion methods (e.g., nuclear-norm penalization), the convergence rate would be \(O_p(\sqrt{K_1(N+T)/p})\), which is much larger than ours. Moreover, even when \(M\) contains an additional component \(M_2+M_3+M_4\) beyond \(M_1\), as long as this component is small (i.e., \(\|M_2+M_3+M_4\|_F \ll \sqrt{(N+T)/p}\)), our rate \[ O_p\!\Big( J + \sqrt{K^*}\,\|M_2+M_3+M_4\|_F \Big) \] is smaller than that of the usual low-rank completion methods \(O_p\!\big(\sqrt{K(N+T)/p}\big)\), where \(K\) is the rank of \(M\) and \(K^*=\max\{K_2,K_3,K_4\}\).
Next, consider the least favorable case, \(M=M_4\), where \(M\) is unrelated to \(X\) and \(Z\). In this case, the convergence rate of our estimator is \(O_p(\sqrt{K_4(N+T)/p})\), provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth. Note that this rate coincides with that of standard low-rank completion methods. Hence, even in the least favorable case, our estimator is comparable to typical low-rank completion methods. Moreover, if \(M\) contains an additional small component \(M_1+M_2+M_3\) beyond \(M_4\), the convergence rate of our estimator becomes \(O_p(\sqrt{K^*(N+T)/p})\), whereas that of a typical low-rank completion method is \(O_p(\sqrt{K(N+T)/p})\). Thus, when \(K_1\) is large, our method can yield a tighter bound. Similar discussions apply to the cases \(M=M_2\) and \(M=M_3\): our estimator attains the same rate as standard low-rank completion methods, and it can yield a better bound when \(M\) also contains an \(M_1\) component with large \(K_1\).
Although the missing-at-random assumption is common in the matrix completion literature, it can be inappropriate for some important applications, such as imputing control potential outcomes in causal panel models, where treatment is assigned to a subset of units starting at a certain time (or in a staggered fashion). In such settings, it may be more appropriate to treat the missingness pattern as fixed (i.e., nonrandom).
Following the literature on matrix completion under missing not at random (MNAR), e.g., bai2021matrix,choi2024matrix,yan2024entrywise, we assume that the missingness pattern takes the form shown in Figure (ref). In this setting, all (or some) entries in the “miss” submatrix are unobserved, while all entries in the “tall” and “wide” submatrices are observed. This pattern is prevalent in causal panel data: the “wide” submatrix corresponds to observations for the control group over all time periods, and the “tall” submatrix corresponds to observations for all units in the pre-treatment period, where the outcome is the potential outcome under control.
\paragraph{Estimation.}
Note that the entries in the “tall” and “wide” submatrices are fully observed. Hence, we can apply Algorithm (ref) to the “tall” and “wide” submatrices to estimate \(M_{\mathrm{tall}} = (m_{it})_{i\leq N,\, t \leq T_0}\) and \(M_{\mathrm{wide}} = (m_{it})_{i\leq N_0,\, t \leq T}\). Then, as noted in Corollary (ref), we can estimate the left and right singular vectors of \(M_{\mathrm{tall}}\) and \(M_{\mathrm{wide}}\), respectively. Importantly, the left singular vectors of \(M_{\mathrm{tall}}\) and \(M\) span the same space. Similarly, the right singular vectors of \(M_{\mathrm{wide}}\) and \(M\) span the same space. Hence, by combining the estimator of the left singular vectors of \(M_{\mathrm{tall}}\) with the estimator of the right singular vectors of \(M_{\mathrm{wide}}\), with an appropriate rotation adjustment, we can estimate \(M\).
Specifically, we estimate \(M\) as follows.
Because this estimator is built on Algorithm (ref) for the fully observed case, we expect it to share similar advantages to those discussed in Section (ref).
\paragraph{Asymptotic result.}
To make this point more precise, we present the convergence rate of our estimator. We begin with an additional assumption.
This assumption can be viewed as an incoherence condition ensuring that the left singular vectors of \(M\) are not dominated by either treated or untreated units, and that the right singular vectors are not dominated by either pre-treatment or post-treatment periods. It is common in the MNAR matrix completion literature (e.g., Assumption D in bai2021matrix and Theorem 3.1(v) in choi2024matrix) and allows us to relate the properties of the submatrices \(M_{\mathrm{tall}}\) and \(M_{\mathrm{wide}}\) to those of \(M\). For example, if \(\{u_i\}_{i \leq N}\) is stationary, then \[ \frac{N}{N_0} \sum_{i \leq N_0} u_i u_i^\top \approx \frac{1}{N}\sum_{i \le N} u_i u_i^\top = I_K, \] and the condition is satisfied.
We are now in a position to state the convergence rate of the estimator. Let \(\mathcal{R}_{\mathrm{tall}}\) and \(\mathcal{R}_{\mathrm{wide}}\) denote the upper bounds on \(\|\widehat{M}_{\mathrm{tall}} - M_{\mathrm{tall}}\|_F\) and \(\|\widehat{M}_{\mathrm{wide}} - M_{\mathrm{wide}}\|_F\), respectively, as given by Theorem (ref):
where \[ P_{X,\mathrm{sub}} = \Phi_{\mathrm{sub}}(X)\big( \Phi_{\mathrm{sub}}(X)^{\top}\Phi_{\mathrm{sub}}(X)\big)^{-1} \Phi_{\mathrm{sub}}(X)^{\top}, \quad P_{Z,\mathrm{sub}} = \Psi_{\mathrm{sub}}(Z)\big( \Psi_{\mathrm{sub}}(Z)^{\top}\Psi_{\mathrm{sub}}(Z)\big)^{-1} \Psi_{\mathrm{sub}}(Z)^{\top}, \] \(\Phi_{\mathrm{sub}}(X) = \big[ \boldsymbol{\phi}(x_1), \ldots, \boldsymbol{\phi}(x_{N_0}) \big]^\top\), and \(\Psi_{\mathrm{sub}}(Z) = \big[ \boldsymbol{\psi}(z_1), \ldots, \boldsymbol{\psi}(z_{T_0}) \big]^\top\). In addition, let \(\delta_N = N_0/N\) and \(\delta_T = T_0/T\). The following theorem provides the convergence rate of our estimator.
Theorem (ref) highlights the advantage of our estimator. A discussion similar to that in Section (ref) applies. For simplicity, consider a typical case in which \(\kappa\) is bounded and \(N_0 \geq cN\) and \(T_0 \geq cT\) for some \(c>0\). First, consider the most favorable case, \(M=M_1\). In this case, the convergence rate of our estimator is \(O_p(J)\) provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth. By contrast, if we estimate the submatrices \(M_{\mathrm{tall}}\) and \(M_{\mathrm{wide}}\) using an MNAR low-rank method (e.g., bai2021matrix), the convergence rate would be \(O_p(\sqrt{K_1(N+T)})\), which is much larger than ours.
Next, consider the least favorable case, \(M=M_4\). In this case, the convergence rate of our estimator is \(O_p(\sqrt{K_4(N+T)})\), provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth and \(J \ll \sqrt{K_4(N+T)}\). This rate coincides with that obtained by applying an MNAR low-rank estimator to the submatrices \(M_{\mathrm{tall}}\) and \(M_{\mathrm{wide}}\).
In addition, if \(M=M_2\), the convergence rate of our estimator is \(O_p(\sqrt{K_2 T})\) provided that \(g_{\iota,k}(\cdot)\) and \(q_{\iota,k}(\cdot)\) are sufficiently smooth and \(J \ll \sqrt{K_2 T}\), whereas a low-rank approach yields the rate \(O_p(\sqrt{K_2(N+T)})\). Similarly, if \(M=M_3\), the convergence rate of our estimator is \(O_p(\sqrt{K_3 N})\), whereas a low-rank approach yields \(O_p(\sqrt{K_3(N+T)})\). Therefore, when \(M\) is at least partially explained by observable characteristics \(X\) and \(Z\), incorporating this information can substantially improve estimation accuracy.
To demonstrate the practical merits and finite-sample performance of our methodology, we conducted several sets of simulation experiments.
First, to study how the relative advantage of our estimator over existing methods varies with the contribution of each component, we change the component weights and compare the estimation performance across estimators. Specifically, we consider the model \[ M = \alpha_1 M_1 + \alpha_2 M_2 + \alpha_3 M_3 + \alpha_4 M_4, \] where \(\sum_{r=1}^4 \alpha_r = 1\) and \(\|M_r\|_F = 2\sqrt{NT}\) for all \(1 \le r \le 4\). We vary the values of \(\alpha_r\) and evaluate the mean squared error of the estimators.
\paragraph{Data generating process.}
We consider eight characteristics, with \(x_i = [x_{1,i}, x_{2,i}, x_{3,i}, x_{4,i}]^\top\) and \(z_t = [z_{1,t}, z_{2,t}, z_{3,t}, z_{4,t}]^\top\). We draw \(x_{1,i} \overset{\mathrm{i.i.d.}}{\sim} \mathrm{Unif}[-1,1]\), \(x_{2,i} \overset{\mathrm{i.i.d.}}{\sim} \mathrm{Unif}[-0.5,0.5]\), \(x_{3,i} \overset{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0,0.2^2)\), and \(x_{4,i} \overset{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0,0.3^2)\). We generate \(z_{1,t}, z_{2,t}, z_{3,t}, z_{4,t}\) in the same way.
For the matrix \(M_1\), we set
and draw the coefficients \(b^{(1,k)}_{0}\) and \(b^{(1,k)}_{d,j}\) from the standard normal distribution. Similarly, we set
and draw the coefficients \(a^{(1,k)}_{0}\) and \(a^{(1,k)}_{d,j}\) from the standard normal distribution.
For the matrix \(M_2\), we set \(K_2=3\) and generate \(G_2(X)\) using the same specification as above. In addition, we generate \(v_{1,t} \in \mathbb{R}^{3}\) i.i.d.\ from \(\mathcal{N}\!\big(0,\textrm{diag}(0.5,1,1.5)\big)\), and stack them into \(V_1 = [v_{1,1}, \ldots, v_{1,T}]^\top\).
For the matrix \(M_3\), we set \(K_3=3\) and generate \(Q_2(Z)\) using the same specification as above. We generate \(w_{1,i} \in \mathbb{R}^{3}\) i.i.d.\ from \(\mathcal{N}\!\big(0,\textrm{diag}(0.5,1,1.5)\big)\), and stack them into \(W_1 = [w_{1,1}, \ldots, w_{1,N}]^\top\).
For the matrix \(M_4\), we draw \(w_{2,i} \in \mathbb{R}^{3}\) i.i.d.\ from \(\mathcal{N}\!\big(0,\textrm{diag}(0.5,1,1.5)\big)\) and \(v_{2,t} \in \mathbb{R}^{3}\) i.i.d.\ from \(\mathcal{N}\!\big(0,1.5^2 I_3\big)\). Lastly, we normalize all four matrices so that \(\|M_r\|_F = 2\sqrt{NT}\) for \(r=1,2,3,4\). We generate the noise entries i.i.d.\ from \(\mathcal{N}(0,0.5^2)\).
For the fully observed case, we set \(N=T=200\). For the MAR (missing at random) case, we set \(N=T=400\) and the observation probability \(p=0.6\). For the MNAR (missing not at random) case, we set \(N=T=400\) and \(N_0=T_0=200\).
\paragraph{Results.}
Here, we use a polynomial sieve with \(J=5\), and we set the number of iterations to \(100\). We vary \(\alpha_r\) under the restrictions \(\sum_{r=1}^4 \alpha_r = 1\) and \(\alpha_2 = \alpha_3\). In addition, to keep the rank of \(M\) constant, we restrict attention to cases in which \(\alpha_r \geq 0.01\) for all \(r=1,2,3,4\).
We first study the fully observed case. Figure (ref) reports the AMSE (average mean squared error) of the estimators. Here, “Oracle” denotes the spectral estimator with known \(K\); “Nuclear norm” denotes the plain nuclear-norm-penalized estimator; “Double projection” denotes \(P_X Y P_Z\); and “Spectral” denotes the spectral estimator with an estimated \(K\). For rank estimation, we use the eigenvalue-ratio method of ahn2013eigenvalue.
From the left panel, we see that, in general, the double projection estimator and the spectral estimator with an estimated rank perform poorly relative to the other estimators. From the right panel, we find that our method performs better than the spectral estimator with known \(K\). The AMSE of the oracle estimator is quite stable and is not sensitive to changes in \(\alpha_r\). In contrast, the AMSEs of our method and the nuclear-norm-penalized estimator are strongly affected by \(\alpha_r\). Overall, our method outperforms nuclear-norm penalization except when \(\alpha_1\) is very small and \(\alpha_4\) is large. When \(\alpha_1\) is large and \(\alpha_4\) is small, our estimator performs particularly well.
To assess the relative advantage of our estimator over others, Figure (ref) plots \((AMSE_{\mathrm{other}} - AMSE_{\mathrm{our}})/AMSE_{\mathrm{our}}\). Relative to nuclear-norm penalization, the advantage of our estimator increases as \(\alpha_1\) increases. Roughly speaking, the advantage also becomes larger as \(\alpha_4\) decreases. In particular, when \(\alpha_1\) is close to zero (e.g., \(\alpha_1 = 0.01\)), the relative performance improves as \(\alpha_2=\alpha_3\) increases and \(\alpha_4\) decreases. In the right panel, which compares our estimator with the oracle estimator, the dependence on \(\alpha_4\) is less clear; nevertheless, we still observe that the relative advantage increases with \(\alpha_1\).
Next, we study the MAR (missing at random) case. Figure (ref) reports the AMSE (average mean squared error) of the estimators as well as their relative performance. Here, we include the double projection method \(\big(p^{-1} P_X(\Omega \circ Y) P_Z\big)\) and the nuclear-norm-penalized estimator, which is a standard approach in the MAR setting. We find that the double projection method performs very poorly except when \(\alpha_1 = 1\). Relative to nuclear-norm penalization, the advantage of our method increases as \(\alpha_1\) increases. However, the pattern with respect to \(\alpha_4\) is less clear than in the fully observed case. This may be because, in the MAR setting, we cannot separately estimate \(M_2\), \(M_3\), and \(M_4\). When \(\alpha_1\) is very small and \(\alpha_4\) is relatively large, nuclear-norm penalization performs better than our method.
Lastly, we study the MNAR (missing not at random) case. Figure (ref) reports the AMSE of the estimators and the relative advantage of our estimator. Here, we assume the rank is known and compare our method with the standard spectral-based estimator for the MNAR setting in bai2021matrix,yan2024entrywise. We find that when \(\alpha_1\) is large and/or \(\alpha_4\) is small, our method generally outperforms the spectral-based estimator. Conversely, when \(\alpha_1\) is small and/or \(\alpha_4\) is large, the spectral-based estimator typically performs better. However, the relative advantage when \(\alpha_1\) is large (or \(\alpha_4\) is small) is substantially greater than the relative disadvantage when \(\alpha_1\) is small (or \(\alpha_4\) is large). Moreover, even when \(\alpha_1\) is close to zero, our method can still perform better when \(\alpha_2\) and \(\alpha_3\) are large.
In summary, across all settings, our method performs substantially better than the competing estimators when \(\alpha_1\) is large and/or \(\alpha_4\) is small. When \(\alpha_1\) is small and/or \(\alpha_4\) is large, the disadvantage of our method is relatively mild compared to the gains achieved when \(\alpha_1\) is large and/or \(\alpha_4\) is small.
Additionally, to examine how our estimator’s relative advantage over existing methods varies with the rank of each component, we vary the ranks and compare estimation performance. Overall, our method performs markedly better than competing estimators when $K_1$ is large. In contrast, when $K_1$ is small and $K_4$ is large, the advantage is more modest and performance is comparable to that of other estimators. For details, please refer to Section (ref) in the Appendix.
In this section, we conduct a real-data experiment using the tobacco sales data in abadie2010synthetic, which is widely used in the literature. In 1988, California introduced the first major anti-tobacco legislation in the United States (Proposition 99). To study the effect of this legislation on tobacco sales, abadie2010synthetic used per-capita cigarette sales data collected from 1970 to 2000 across 38 U.S.\ states with no anti-tobacco legislation prior to 2000 (\(N=38, T=31\)). We encode these data into a \(38 \times 31\) matrix \(Y\), where the entry \(y_{it}\) represents the “potential” outcome of per-capita cigarette sales (in packs) for state \(i\) in year \(t\) under “control,” i.e., in the absence of any intervention. To generate missing entries, we artificially introduce interventions (i.e., missingness) for a subset of states: in each iteration, we randomly select 8 states to adopt an intervention (e.g., a tobacco control program) starting from period \(T_0+1\). After rearranging the matrix, this yields the block-missing pattern shown in Figure (ref), with an \(8 \times (T-T_0)\) missing submatrix.
For state-level characteristics, we use the time-averaged retail price of cigarettes, log per-capita state personal income, the percentage of the population aged 18--24, the percentage of adults completing four years of college or more, and per-capita beer consumption. Most of these variables are averaged over the 1970--2000 period. In addition, as a proxy for a state's general preference for tobacco, we use per-capita cigarette sales in 2001. For year-level characteristics, we use log per-capita real GDP and the state-average retail price of cigarettes in each year. We also use the average per-capita cigarette sales of Florida and Michigan as a proxy for general tobacco preference in each year. Although Florida and Michigan are not included among the above 38 states because they had interventions before 2000, the effects of those interventions were relatively mild compared to other treated states. Appendix Section (ref) provides additional details on the construction of these characteristics.
We compare the performance of our estimator with that of the spectral-based estimator in bai2021matrix,yan2024entrywise, which is a standard method for block-missing patterns. This approach estimates the “tall” and “wide” submatrices using a spectral estimator. For rank estimation, we use the eigenvalue-ratio method of ahn2013eigenvalue. For the projection step in our method, we use a second-order polynomial sieve (\(J=2\)). We set the number of iterations to 100.
We first compare the AMSE (average mean squared error) over the missing entries. Specifically, in each iteration, we sum the squared estimation errors over all \(8(T-T_0)\) missing entries and divide by \(8(T-T_0)\), and then average this quantity across iterations. We compute the estimation error as the difference between the estimated value and the observed per-capita cigarette sales for each missing entry. The first two rows of Table (ref) report the results for different adoption times \(T_0\). We find that our method outperforms the spectral-based estimator in all cases. In particular, when the number of observed periods is relatively small (i.e., \(T_0\) is small), the performance gap is larger.
Moreover, we consider the AMSE of (i) the average of the missing entries in each year and (ii) the average of all missing entries. For the AMSE of the year-by-year averages, in each iteration we compute the average of the missing entries in each post-intervention year, compute the squared estimation error for each such average, sum these squared errors, and divide by \(T-T_0\). We then average this quantity across iterations. For the AMSE of the overall average, in each iteration we compute the average of all missing entries, compute its squared estimation error, and then average it across iterations. These average-type targets have the advantage that the noise in outcomes is averaged out; therefore, averages of \(y_{it}\) are close to averages of \(m_{it}\). The last four rows of Table (ref) report the results for different adoption times \(T_0\). We find that our method outperforms the spectral-based estimator, and the performance gap increases when the number of observed periods is small (i.e., when \(T_0\) is small).
As an alternative, we also consider different proxies for tobacco preference. For the proxy of each state's general preference for tobacco, we use the time average of per-capita cigarette sales over the full sample period when estimating the “wide” submatrix, and we use the average over the pre-intervention period \(1,\ldots,T_0\) when estimating the “tall” submatrix. Similarly, for the proxy of each year's general preference for tobacco, we use the average per-capita cigarette sales across the 30 control states in each year when estimating the “wide” submatrix, and we use the average across all 38 states in each year when estimating the “tall” submatrix.
Table (ref) reports the results when we use these alternative proxies as characteristics. We find that the performance of our method improves in most cases, and its relative advantage becomes larger. This type of proxy is not fully consistent with the theory because it may violate the exogeneity condition; however, because the proxy averages out outcome noise, the resulting endogeneity may be negligible in practice. In our experiment, the results using these proxies are indeed favorable.
In summary, the empirical results suggest that incorporating side information can improve estimation of the control potential outcomes in causal panel settings, relative to typical low-rank methods.
This paper proposes a flexible framework for high-dimensional matrix estimation that systematically incorporates rich side information on both rows and columns. By decomposing the signal into components explained jointly by $(X, Z)$, by $X$ alone, by $Z$ alone, and by neither, and by estimating these components using sieve projection combined with nuclear-norm penalization, our approach accommodates nonlinear covariate effects, avoids explicit rank selection for each component, and automatically thresholds weak or negligible signals. We establish convergence rates that demonstrate robustness across diverse model configurations, matching specialized procedures in favorable settings while remaining competitive when side information is uninformative. We further extend the method to partially observed matrices under both MAR and MNAR mechanisms, including block-missing patterns motivated by causal panel data, and show through simulations and a tobacco-sales application that leveraging side information can substantially improve imputation accuracy and enhance treatment-effect estimation.