One-Sided Matrix Completion from Ultra-Sparse Samples

TL;DR

This work tackles one-sided matrix completion in ultra-sparse regimes where only $C n$ entries are observed from an $n\times d$ matrix with $n\ge d$, making full recovery of $M$ impossible but enabling recovery of the row-space via the second-moment $T = M^{\top}M / n$. The authors introduce the Hájek estimator for $T$, show it is exactly unbiased on the observed subset and offers variance reduction over Horvitz-Thompson, and couple it with gradient-descent imputation in a rank-$r$ incoherent model to recover $T$ with provable guarantees: if $n \ge c d r^5 \kappa^6 \mu^2 \log d /(C^2 \varepsilon^2)$, any local minimum satisfies $\|XX^{\top}-T\|_F^2 \le \varepsilon^2$. They further provide extensions to full matrix imputation and differential privacy, and validate the approach on synthetic and real-world datasets (MovieLens, Amazon Reviews, Genomes), showing substantial bias reductions and improved recovery with near-linear runtime scaling. The work offers a practical, scalable path to extracting row-space information from ultra-sparse panel data, with implications for privacy-preserving summaries and downstream imputation tasks. Overall, the Hájek-GD method advances one-sided matrix completion by delivering finite-sample unbiasedness, variance reduction, and tractable sample complexity in non-random missingness settings.

Abstract

Matrix completion is a classical problem that has received recurring interest across a wide range of fields. In this paper, we revisit this problem in an ultra-sparse sampling regime, where each entry of an unknown, $n\times d$ matrix $M$ (with $n \ge d$) is observed independently with probability $p = C / d$, for a fixed integer $C \ge 2$. This setting is motivated by applications involving large, sparse panel datasets, where the number of rows far exceeds the number of columns. When each row contains only $C$ entries -- fewer than the rank of $M$ -- accurate imputation of $M$ is impossible. Instead, we estimate the row span of $M$ or the averaged second-moment matrix $T = M^{\top} M / n$. The empirical second-moment matrix computed from observed entries exhibits non-random and sparse missingness. We propose an unbiased estimator that normalizes each nonzero entry of the second moment by its observed frequency, followed by gradient descent to impute the missing entries of $T$. The normalization divides a weighted sum of $n$ binomial random variables by the total number of ones. We show that the estimator is unbiased for any $p$ and enjoys low variance. When the row vectors of $M$ are drawn uniformly from a rank-$r$ factor model satisfying an incoherence condition, we prove that if $n \ge O({d r^5 ε^{-2} C^{-2} \log d})$, any local minimum of the gradient-descent objective is approximately global and recovers $T$ with error at most $ε^2$. Experiments on both synthetic and real-world data validate our approach. On three MovieLens datasets, our algorithm reduces bias by $88\%$ relative to baseline estimators. We also empirically validate the linear sampling complexity of $n$ relative to $d$ on synthetic data. On an Amazon reviews dataset with sparsity $10^{-7}$, our method reduces the recovery error of $T$ by $59\%$ and $M$ by $38\%$ compared to baseline methods.

Figures (7)

• Figure 1: An illustration of the variance reduction using the Hájek estimator vs. the Horvitz-Thompson (HT) estimator, measured on synthetic data with $n=10^4$ and $d=10^3$, by repeating the data sampling procedure $100$ times and calculating the variance across either uniform sampling with probability $p$, or sampling $C$ entries per row. In particular, $Var(\widehat{T})$ is generally $10^{-3}$ lower than $Var(\overline{T})$. For details regarding the simulation setup, see Section \ref{['sec_experiments']}.
• Figure 2: We illustrate the results from applying a first-order approximation to the variance of the diagonal entries of Hájek-GD, run on synthetic data with $n =10^4$ and $d = 10^3$. Figures \ref{['fig_taylor_iid_mean']} and \ref{['fig_taylor_iid_var']}: Sample each entry with probability $p$. Figures \ref{['fig_taylor_fix_mean']} and \ref{['fig_taylor_fix_var']}: Sample $C$ entries per row without repetition. In particular, the approximation errors incurred by both first-order Taylor's expansion and the variance approximation are generally less than $10^{-6}$.
• Figure 3: Illustration of the proposed one-sided matrix completion pipeline.Input: A partially observed matrix $\hat{M}\in\mathbb{R}^{n\times d}$. Each row represents a user's data (e.g., movie ratings), and each column represents an item. Thus, $M$ is a tall-and-skinny matrix, with $n$ being larger than $d$. Step 1: Apply the Hájek estimator to correct for non-random missingness in the empirical second-moment matrix $\hat{M}^{\top} \hat{M}$ over the observed index set $\Omega$. Step 2: Impute the remaining unobserved entries of $T$ by running gradient descent on a low-rank reconstruction loss between the estimated and true second-moment matrices. The final $T$ combines the observed entries from $\widehat{T}$ with the missing entries from $X_t X_t^{\top}$, where $X_t$ denotes the $t$-th gradient descent iterate.Extension: Perform full matrix imputation by projecting onto the recovered low-rank subspace $U_r$ (obtained via rank-$r$ SVD of the estimated $T$) and solving a least-squares regression problem to estimate the missing entries of $M$.
• Figure 4: Comparing the bias of the Hájek estimator vs. the Horvitz-Thompson estimator, measured as the mean squared error between $\widehat{T}$ (or $\overline T$) and $T$ on the index set $\Omega$. Figure \ref{['fig_T_hat_bias']}: Reporting the bias on synthetic matrix data with $n = 10^4$ and $d = 10^3$. Figure \ref{['fig_T_hat_ml']}: Reporting the bias on MovieLens datasets. We consider both uniform sampling with probability $p$ as well as sampling $C$ entries from each row. In Figure \ref{['fig_T_hat_ml']}, we further consider a non-uniform sampling setting, where the sampling probability at each row is proportional to the number of nonzero entries in that row.
• Figure 5: Illustration of the recovery error of Algorithm \ref{['alg_ipw']}. Figure \ref{['fig_sample_complexity_fix_rank']}: We vary the sampling probability and find that the error roughly stays constant after we also fix the number of samples. Figure \ref{['fig_sample_complexity_fix_p']}: We see similar results as we fix the number of samples but vary the rank of the underlying matrix. Figure \ref{['fig_sample_complexity_sigma']}: We gradually increase the noise level and find that the estimation errors also increase. The reported results are aggregated based on five independent runs.

Theorems & Definitions (35)

• Claim 2.1
• Lemma 3.1
• Theorem 3.2
• Theorem 3.4
• Remark 3.5
• Lemma 4.1
• Lemma 4.2
• Remark 4.3
• proof : Proof of Lemma \ref{['prop_consistency']}
• proof : Proof of Theorem \ref{['prop_var_est']}