Matrix Completion under Low-Rank Missing Mechanism

TL;DR

This paper introduces a novel approach for matrix completion that accounts for a low-rank missing data mechanism, estimating observation probabilities via high-dimensional low-rank matrix estimation and inverse probability weighting, with proven optimal convergence rates.

Contribution

It develops a new framework for matrix completion under low-rank missing mechanisms, extending beyond uniform missing assumptions with theoretical guarantees.

Findings

Derived optimal convergence rates for probability and matrix estimators

Proposed a high-dimensional low-rank estimation procedure for missing data

Analyzed the effects of extreme missing probabilities on estimation accuracy

Abstract

Matrix completion is a modern missing data problem where both the missing structure and the underlying parameter are high dimensional. Although missing structure is a key component to any missing data problems, existing matrix completion methods often assume a simple uniform missing mechanism. In this work, we study matrix completion from corrupted data under a novel low-rank missing mechanism. The probability matrix of observation is estimated via a high dimensional low-rank matrix estimation procedure, and further used to complete the target matrix via inverse probabilities weighting. Due to both high dimensional and extreme (i.e., very small) nature of the true probability matrix, the effect of inverse probability weighting requires careful study. We derive optimal asymptotic convergence rates of the proposed estimators for both the observation probabilities and the target matrix.

Tables

Table 1. Table 1: Root mean squared errors RMSE ​ ( 𝑴 ^ , 𝑴 ⋆ ) RMSE ^ 𝑴 subscript 𝑴 ⋆ \text{RMSE}(\widehat{\bm{M}},\bm{M}_{\star}) , Hellinger distance d H 2 ​ ( 𝚯 ^ , 𝚯 ⋆ ) superscript subscript 𝑑 𝐻 2 ^ 𝚯 subscript 𝚯 ⋆ d_{H}^{2}(\widehat{\bm{\Theta}},\bm{\Theta}_{\star}) , rank of linear predictor 𝑴 ^ ^ 𝑴 \widehat{\bm{M}} and estimated 𝚯 ^ ^ 𝚯 \widehat{\bm{\Theta}} and their standard errors (in parentheses) under the low rank missing observation mechanism, with ( n 1 , n 2 subscript 𝑛 1 subscript 𝑛 2 n_{1},n_{2} ) = (600, 600), (800, 800), (1000, 1000), (1200, 1200) and r 𝑴 ⋆ = 11 subscript 𝑟 subscript 𝑴 ⋆ 11 r_{\bm{M}_{\star}}=11 , for the proposed estimators 𝚯 ^ α subscript ^ 𝚯 𝛼 \widehat{\bm{\Theta}}_{\alpha} , 𝚯 ^ 1-bit , α subscript ^ 𝚯 1-bit 𝛼 \widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha} and the two existing estimators ( 𝚯 ^ NW subscript ^ 𝚯 NW \widehat{\bm{\Theta}}_{\textsf{NW}} and 𝚯 ^ UNI subscript ^ 𝚯 UNI \widehat{\bm{\Theta}}_{\textsf{UNI}} ).

600	${\widehat{𝚯}}_{\alpha }$	${\widehat{𝚯}}_{\text{1-bit},\alpha }$	${\widehat{𝚯}}_{\text{NW}}$	${\widehat{𝚯}}_{\text{UNI}}$
RMSE($\widehat{𝑴},{𝑴}_{\star }$)	2.6923 (0.0342)	2.9155 (0.0295)	-	-
$d_H^2(\widehat{𝚯},{𝚯}_{\star })$	0.0369 (0.0015)	0.0450 (0.0016)	0.1233 (1e-04)	0.1729 (1e-04)
$r_{\widehat{𝑴}}$	12.45 (0.50)	12.69 (0.46)	-	-
$r_{\widehat{𝚯}}$	600.00 (0.00)	600.00 (0.00)	-	-
800	${\widehat{𝚯}}_{\alpha }$	${\widehat{𝚯}}_{\text{1-bit},\alpha }$	${\widehat{𝚯}}_{\text{NW}}$	${\widehat{𝚯}}_{\text{UNI}}$
RMSE($\widehat{𝑴},{𝑴}_{\star }$)	2.5739 (0.0116)	2.7796 (0.0033)	-	-
$d_H^2(\widehat{𝚯},{𝚯}_{\star })$	0.0317 (5e-04)	0.0379 (1e-04)	0.1219 (1e-04)	0.1767 (1e-04)
$r_{\widehat{𝑴}}$	12.04 (0.20)	12.03 (0.17)	-	-
$r_{\widehat{𝚯}}$	800.00 (0.00)	800.00 (0.00)	-	-
1000	${\widehat{𝚯}}_{\alpha }$	${\widehat{𝚯}}_{\text{1-bit},\alpha }$	${\widehat{𝚯}}_{\text{NW}}$	${\widehat{𝚯}}_{\text{UNI}}$
RMSE($\widehat{𝑴},{𝑴}_{\star }$)	2.4870 (0.0212)	2.7731 (0.0015)	-	-
$d_H^2(\widehat{𝚯},{𝚯}_{\star })$	0.0266 (8e-04)	0.0351 (1e-04)	0.1246 (1e-04)	0.1767 (1e-04)
$r_{\widehat{𝑴}}$	12.68 (0.53)	12.00 (0.00)	-	-
$r_{\widehat{𝚯}}$	1000.00 (0.00)	1000.00 (0.00)	-	-
1200	${\widehat{𝚯}}_{\alpha }$	${\widehat{𝚯}}_{\text{1-bit},\alpha }$	${\widehat{𝚯}}_{\text{NW}}$	${\widehat{𝚯}}_{\text{UNI}}$
RMSE($\widehat{𝑴},{𝑴}_{\star }$)	2.3809 (0.0018)	2.6470 (0.0012)	-	-
$d_H^2(\widehat{𝚯},{𝚯}_{\star })$	0.0242 (1e-04)	0.0314 (1e-04)	0.1211 (1e-04)	0.1761 (1e-04)
$r_{\widehat{𝑴}}$	12.00 (0.00)	12.00 (0.00)	-	-
$r_{\widehat{𝚯}}$	1200.00 (0.00)	1200.00 (0.00)	-	-

Table 2. Table 2: Root mean squared errors, test errors, estimated ranks r 𝑨 ^ β subscript 𝑟 subscript ^ 𝑨 𝛽 r_{\widehat{\bm{A}}_{\beta}} and their standard deviations (in parentheses) under the low rank missing observation mechanism, for three existing methods and ten versions of the proposed methods where Proposed indicates the estimators are obtained by solving problem ( 3 ), while 𝚯 ^ β subscript ^ 𝚯 𝛽 \widehat{\bm{\Theta}}_{\beta} , 𝚯 ^ Win , β subscript ^ 𝚯 Win 𝛽 \widehat{\bm{\Theta}}_{\textsf{Win},\beta} , 𝚯 ^ α subscript ^ 𝚯 𝛼 \widehat{\bm{\Theta}}_{\alpha} , 𝚯 ^ 1-bit , β subscript ^ 𝚯 1-bit 𝛽 \widehat{\bm{\Theta}}_{\textsf{1-bit},\beta} , 𝚯 ^ 1-bit , Win , β subscript ^ 𝚯 1-bit Win 𝛽 \widehat{\bm{\Theta}}_{\textsf{1-bit},\textsf{Win},\beta} and 𝚯 ^ 1-bit , α subscript ^ 𝚯 1-bit 𝛼 \widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha} represents the probability estimators used in ( 3 ), as described in Section 6.1 , and t = 0.05 𝑡 0.05 t=0.05 or 0.1 0.1 0.1 denote the winsorized proportion for which β 𝛽 \beta is chosen.

($n_1,n_2$) = (600, 600)	$\text{RMSE}({\widehat{𝑨}}_{\beta },{𝑨}_{\star })$	Test Error	$r_{{\widehat{𝑨}}_{\beta }}$
Proposed_${\widehat{𝚯}}_{\text{Win},\beta }\_0.05$	1.5615 (0.0147)	0.3005 (0.0062)	65.28 (5.72)
Proposed_${\widehat{𝚯}}_{\beta }\_0.05$	1.5548 (0.0085)	0.2996 (0.0034)	54.98 (3.01)
Proposed_${\widehat{𝚯}}_{\text{Win},\beta }\_0.1$	1.5621 (0.0111)	0.3013 (0.0046)	63.68 (5.36)
Proposed_${\widehat{𝚯}}_{\beta }\_0.1$	1.5509 (0.0085)	0.2983 (0.0034)	53.13 (2.72)
Proposed_${\widehat{𝚯}}_{\alpha }$	1.5637 (0.0147)	0.3010 (0.0061)	65.63 (5.89)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\text{Win},\beta }\_0.05$	1.5664 (0.0093)	0.3028 (0.0037)	62.76 (5.96)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\beta }\_0.05$	1.5573 (0.0089)	0.2996 (0.0036)	61.80 (5.34)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\text{Win},\beta }\_0.1$	1.5669 (0.0092)	0.3032 (0.0037)	62.78 (2.68)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\beta }\_0.1$	1.5540 (0.0089)	0.2987 (0.0036)	60.79 (3.01)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\alpha }$	1.5612 (0.0097)	0.3005 (0.0040)	62.12 (4.76)
NW	1.9896 (0.2814)	0.4676 (0.1341)	167.67 (54.78)
KLT	2.2867 (0.0073)	0.5951 (0.0026)	1.00 (0.00)
MHT	1.6543 (0.0097)	0.3432 (0.0041)	51.20 (2.61)
($n_1,n_2$) = (800, 800)	$\text{RMSE}({\widehat{𝑨}}_{\beta },{𝑨}_{\star })$	Test Error	$r_{{\widehat{𝑨}}_{\beta }}$
Proposed_${\widehat{𝚯}}_{\text{Win},\beta }\_0.05$	1.4754 (0.0107)	0.2669 (0.0041)	88.58 (10.81)
Proposed_${\widehat{𝚯}}_{\beta }\_0.05$	1.4797 (0.0080)	0.2714 (0.0030)	71.79 (4.12)
Proposed_${\widehat{𝚯}}_{\text{Win},\beta }\_0.1$	1.4724 (0.0108)	0.2664 (0.0042)	86.25 (10.34)
Proposed_${\widehat{𝚯}}_{\beta }\_0.1$	1.4763 (0.0082)	0.2704 (0.0031)	67.08 (4.22)
Proposed_${\widehat{𝚯}}_{\alpha }$	1.4783 (0.0115)	0.2676 (0.0041)	88.92 (11.70)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\text{Win},\beta }\_0.05$	1.4917 (0.0078)	0.2743 (0.0030)	83.51 (1.45)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\beta }\_0.05$	1.4804 (0.0080)	0.2705 (0.0031)	82.60 (3.47)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\text{Win},\beta }\_0.1$	1.4972 (0.0080)	0.2765 (0.0031)	81.64 (7.23)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\beta }\_0.1$	1.4800 (0.0078)	0.2708 (0.0030)	74.89 (3.54)
Proposed_${\widehat{𝚯}}_{\text{1-bit},\alpha }$	1.4790 (0.0099)	0.2685 (0.0039)	88.57 (9.56)
NW	1.9515 (0.3625)	0.4585 (0.1593)	215.61 (82.24)
KLT	2.3447 (0.0064)	0.6081 (0.0020)	1.00 (0.00)
MHT	1.6067 (0.0086)	0.3245 (0.0036)	63.68 (3.02)

Table 3. Table 3: Root mean squared prediction errors based on Yahoo! Webscope dataset for the ten versions of the proposed method and the three existing methods proposed respectively in Negahban and Wainwright ( 2012 ) ( NW ), Koltchinskii et al. ( 2011 ) ( KLT ) and Mazumder et al. ( 2010 ) ( MHT ).

	Proposed_${\widehat{𝚯}}_{\text{Win},\beta }$_0.05	Proposed_${\widehat{𝚯}}_{\beta }$_0.05	Proposed_${\widehat{𝚯}}_{\text{Win},\beta }$_0.1
RMSPE	1.0396	1.0381	1.0476
	Proposed_${\widehat{𝚯}}_{\beta }$_0.1	Proposed_${\widehat{𝚯}}_{\alpha }$	Proposed_${\widehat{𝚯}}_{\text{1-bit},\text{Win},\beta }$_0.05
RMSPE	1.0490	1.0383	1.0831
	Proposed_${\widehat{𝚯}}_{\text{1-bit},\beta }$_0.05	Proposed_${\widehat{𝚯}}_{\text{1-bit},\text{Win},\beta }$_0.1	Proposed_${\widehat{𝚯}}_{\text{1-bit},\beta }$_0.1
RMSPE	1.1091	1.0760	1.0523
	Proposed_${\widehat{𝚯}}_{\text{1-bit},\alpha }$	NW	KLT
RMSPE	1.1065	1.7068	3.6334
	MHT
RMSPE	1.3821

Equations

$$y_{ij}=a_{\star ,ij}+{ϵ}_{ij},\text{for i=1,…,n1;j=1,…,n2,}$$

$$\bm{M}_{\star}=\mu_{\star}\bm{J}+\bm{Z}_{\star}\quad\hbox{ where}\quad\mbox{$\bm{1}_{n_{1}}^{{\mathrm{\scriptscriptstyle T}}}\bm{Z}_{\star}\bm{1}_{n_{2}}=0$}$$

$$R\left(\mathbf{A}\right)=\frac{1}{n_1{n}_2}\text{E}\left({\Vert \mathbf{A}-\mathbf{Y}\Vert }_F^2\right)=\frac{1}{n_1{n}_2}\text{E}\left({\Vert \mathbf{W}\circ {\mathbf{\Theta }}_{\star }^{‡}\circ \left(\mathbf{A}-\mathbf{Y}\right)\Vert }_F^2\right).$$

$$\hat{R}\left(\mathbf{A}\right)=\frac{1}{n_1{n}_2}{\Vert \mathbf{W}\circ {\mathbf{\Theta }}_{\star }^{‡}\circ \left(\mathbf{A}-\mathbf{Y}\right)\Vert }_F^2,$$

$${\ell }_{\mathbf{W}}\left(\mathbf{M}\right)=\sum _{i,j}\left[\text{mathds}{I}_{\left[w_{ij}=1\right]}\log \left\{f\left(m_{ij}\right)\right\}+\text{mathds}{I}_{\left[w_{ij}=0\right]}\log \left\{1-f\left(m_{ij}\right)\right\}\right],$$

$$\left(\hat{\mu },\hat{\mathbf{Z}}\right)=$$

$${\mathcal{C}}_{n_1,n_2}\left({\alpha }_1,{\alpha }_2\right)=\{\left(\mu ,\mathbf{Z}\right)\in \mathbb{R}\times {\mathbb{R}}^{n_1\times n_2}:\left|\mu \right|\le {\alpha }_1,{\Vert \mathbf{Z}\Vert }_{\infty }\le {\alpha }_2,1_{n_1}^{\mathrm{T}}\mathbf{Z}{1}_{n_2}=0\}.$$

$$P_L\left\{\left(\mu ,\mathbf{Z}\right),\left({\mu }_{\text{old}},{\mathbf{Z}}_{\text{old}}\right)\right\}=$$

$$+\left(\mu -{\mu }_{\text{old}}\right)1_{n_1}^{\mathrm{T}}\left\{-{\nabla }_{\mu }{\ell }_{\mathbf{W}}\left({\mu }_{\text{old}}\mathbf{J}+{\mathbf{Z}}_{\text{old}}\right)\right\}{1}_{n_2}+\frac{L{n}_1{n}_2}{2}{\left(\mu -{\mu }_{\text{old}}\right)}^2$$

$$+⟨\mathbf{Z}-{\mathbf{Z}}_{\text{old}},-{\nabla }_{\mathbf{Z}}{\ell }_{\mathbf{W}}\left({\mu }_{\text{old}}\mathbf{J}+{\mathbf{Z}}_{\text{old}}\right)⟩+\frac{L}{2}{\Vert \mathbf{Z}-{\mathbf{Z}}_{\text{old}}\Vert }_F^2+\lambda {\Vert \mathbf{Z}\Vert }_{\ast },$$

$$\underset{\left(\mu ,\mathbf{Z}\right)\in {\mathcal{C}}_{n_1,n_2}\left({\alpha }_1,{\alpha }_2\right)}{\mathrm{argmin}}{P}_L\left\{\left(\mu ,\mathbf{Z}\right),\left({\mu }_{\text{old}},{\mathbf{Z}}_{\text{old}}\right)\right\},$$

$$\min \left[{\alpha }_1,\max \left[-{\alpha }_1,{\mu }_{\text{old}}+{\left(L{n}_1{n}_2\right)}^{-1}{1}_{n_1}^{\mathrm{T}}\left\{-{\nabla }_{\mu }{\ell }_{\mathbf{W}}\left({\mu }_{\text{old}}\mathbf{J}+{\mathbf{Z}}_{\text{old}}\right)\right\}{1}_{n_2}\right]\right].$$

$$\underset{{\Vert \mathbf{Z}\Vert }_{\infty }\le {\alpha }_2,1_{n_1}^{\mathrm{T}}\mathbf{Z}{1}_{n_2}=0}{\mathrm{argmin}}⟨\mathbf{Z}-{\mathbf{Z}}_{\text{old}},-{\nabla }_{\mathbf{Z}}{\ell }_{\mathbf{W}}\left({\mu }_{\text{old}}\mathbf{J}+{\mathbf{Z}}_{\text{old}}\right)⟩+\frac{L}{2}{\Vert \mathbf{Z}-{\mathbf{Z}}_{\text{old}}\Vert }_F^2+\lambda {\Vert \mathbf{Z}\Vert }_{\ast },$$

$$\underset{{\Vert \mathbf{Z}\Vert }_{\infty }\le {\alpha }_2,1_{n_1}^{\mathrm{T}}\mathbf{Z}{1}_{n_2}=0}{\mathrm{argmin}}\frac{1}{2}{\Vert \mathbf{Z}-{\mathbf{Z}}_{\text{old}}-\frac{1}{L}{\nabla }_{\mathbf{Z}}{\ell }_{\mathbf{W}}\left({\mu }_{\text{old}}\mathbf{J}+{\mathbf{Z}}_{\text{old}}\right)\Vert }_F^2+\frac{\lambda }{L}{\Vert \mathbf{Z}\Vert }_{\ast }.$$

$$\underset{\mathbf{Z}={\mathbf{G}}_1={\mathbf{G}}_2,1_{n_1}^{\mathrm{T}}{\mathbf{G}}_1{1}_{n_2}=0,{\Vert {\mathbf{G}}_2\Vert }_{\infty }\le {\alpha }_2}{\mathrm{argmin}}\frac{\lambda }{L}{\Vert \mathbf{Z}\Vert }_{\ast }+\frac{1}{2}{\Vert {\mathbf{G}}_2-{\mathbf{Z}}_{\text{old}}-\frac{1}{L}{\nabla }_{\mathbf{Z}}{\ell }_{\mathbf{W}}\left({\mu }_{\text{old}}\mathbf{J}+{\mathbf{Z}}_{\text{old}}\right)\Vert }_F^2.$$

$${\mathcal{L}}_u\left(\mathbf{Z},{\mathbf{G}}_1,{\mathbf{G}}_2;\mathbf{H}\right)=$$

$$-⟨{\mathbf{H}}_1,\mathbf{Z}-{\mathbf{G}}_1⟩-⟨{\mathbf{H}}_2,\mathbf{Z}-{\mathbf{G}}_2⟩+\frac{u}{2}{\Vert \mathbf{Z}-{\mathbf{G}}_1\Vert }_F^2+\frac{u}{2}{\Vert \mathbf{Z}-{\mathbf{G}}_2\Vert }_F^2$$

$$+\text{mathds}{I}_{\left[1_{n_1}^{\mathrm{T}}{\mathbf{G}}_1{1}_{n_2}=0\right]}+\text{mathds}{I}_{\left[{\Vert {\mathbf{G}}_2\Vert }_{\infty }\le {\alpha }_2\right]},$$

$${\hat{\mathbf{Z}}}_{\beta }=\underset{\mathbf{Z}\in {\mathbb{R}}^{n_1\times n_2}}{\mathrm{argmax}}{\ell }_{\mathbf{W}}\left(\hat{\mu }\mathbf{J}+\mathbf{Z}\right)-{\lambda }^{\prime }{\Vert \mathbf{Z}\Vert }_{\ast }\mathrm{subject}\mathrm{to}{\Vert \mathbf{Z}\Vert }_{\infty }\le \beta .$$

$$\tilde{R}\left(\mathbf{A}\right)=\frac{1}{n_1{n}_2}{\Vert \mathbf{W}\circ {\hat{\mathbf{\Theta }}}_{\beta }^{‡}\circ \left(\mathbf{A}-\mathbf{Y}\right)\Vert }_F^2,$$

$$\tilde{R}\left(\mathbf{A}\right)+\tau {\Vert \mathbf{A}\Vert }_{\ast },$$

$${\hat{\mathbf{A}}}_{\beta }=\underset{{\Vert \mathbf{A}\Vert }_{\infty }\le a}{\mathrm{argmin}}\left\{\frac{1}{n_1{n}_2}{\Vert \mathbf{W}\circ {\hat{\mathbf{\Theta }}}_{\beta }^{‡}\circ \left(\mathbf{A}-\mathbf{Y}\right)\Vert }_F^2+\tau {\Vert \mathbf{A}\Vert }_{\ast }\right\},$$

$$\left(\hat{\mu },\hat{\mathbf{Z}}\right)=$$

$${\tilde{\mathcal{C}}}_{n_1,n_2}\left({\alpha }_1,{\alpha }_2\right)=$$

$${\Vert \mathbf{Z}\Vert }_{\ast }\le {\alpha }_2\sqrt{r_{{\mathbf{Z}}_{\star }}{n}_1{n}_2},1_{n_1}^{\mathrm{T}}\mathbf{Z}{1}_{n_2}=0\}.$$

$${\left({\mu }_{\star }-\hat{\mu }\right)}^2\le C_2\left({\alpha }_1^2\wedge {\Gamma }_n\right),\frac{1}{n^2}{\Vert \hat{\mathbf{Z}}-{\mathbf{Z}}_{\star }\Vert }_F^2\le C_2\left({\alpha }_2^2\wedge {\Gamma }_n\right)$$

$$\text{mbox}and\frac{1}{n^2}{\Vert \hat{\mathbf{M}}-{\mathbf{M}}_{\star }\Vert }_F^2\le C_2\left({\alpha }_0^2\wedge {\Gamma }_n\right).$$

$${\hat{\mathbf{Z}}}_{\beta }=\underset{\mathbf{Z}\in {\mathbb{R}}^{n_1\times n_2}}{\mathrm{argmax}}{\ell }_{\mathbf{W}}\left(\hat{\mu }\mathbf{J}+\mathbf{Z}\right)$$

$$\mathrm{subject}\mathrm{to}{\Vert \mathbf{Z}\Vert }_{\infty }\le \beta ,{\Vert \mathbf{Z}\Vert }_{\ast }\le \beta \sqrt{r_{{\mathcal{T}}_{\beta }({\mathbf{Z}}_{\star })}{n}_1{n}_2},.$$

$$T_{\beta }(b_{ij})=b_{ij}\text{mathds}{I}_{[-\beta \le b_{ij}\le \beta ]}+\beta \text{mathds}{I}_{[b_{ij}>\beta ]}-\beta \text{mathds}{I}_{[b_{ij}<-\beta ]}\text{for any }\beta \ge 0.$$

$$d^2\left\{{\hat{\mathbf{Z}}}_{\beta },{\mathcal{T}}_{\beta }\left({\mathbf{Z}}_{\star }\right)\right\}\le C_3{\Lambda }_n,d^2\left({\hat{\mathbf{M}}}_{\beta },{\mathbf{M}}_{\star }\right)\le C_2\left({\alpha }_1^2\wedge {\Gamma }_n\right)+C_3{\Lambda }_n+\frac{2({\alpha }_2-\beta {)}_{+}^2{N}_{\beta }}{n^2}$$

$$\text{mbox}and{d}^2\left({\hat{\mathbf{\Theta }}}_{\beta }^{†},{\mathbf{\Theta }}_{\star }^{†}\right)\le \frac{C_2}{h_{{\alpha }_1,\beta }^2}\left({\alpha }_1^2\wedge {\Gamma }_n\right)+\frac{C_3{\Lambda }_n}{h_{{\alpha }_1,\beta }^2}+\frac{8N_{\beta }}{n^2{\theta }_L^2}.$$

$$k_{{\alpha }_1,{\alpha }_2,\beta ,n}\asymp \min \left[{\beta }^2,h_{{\alpha }_1,\beta }^{-1}\beta \left\{8N_{\beta }+\left(n^2-N_{\beta }\right)k_{{\alpha }_1,{\alpha }_2,n}^{\prime 1/2}\right\}{n}^{-2}+h_{{\alpha }_1,\beta }^{-1}{n}^{-1/2}({\alpha }_1+\beta r_{{\mathcal{T}}_{\beta }({\mathbf{Z}}_{\star })}^{1/2})\right].$$

$$\Delta =\max \left\{\frac{\left(c_{\sigma }\vee a\right)e^{-{\mu }_{\star }/2+{\alpha }_2-\beta +|{\alpha }_2/2-\beta |}{\left(n\log n\right)}^{1/2}}{n^2},\frac{\eta e^{{\mu }_{\star }/2+{\alpha }_1+|{\alpha }_2/2-\beta |}{k}_{{\alpha }_1,{\alpha }_2,\beta ,n}^{1/2}{\log }^{3/2}n}{h_{{\alpha }_1,\beta }n}\right\}.$$

Full text

keywords: Low-rank; Missing; Nuclear-norm; Regularization.

1 Introduction

The problem of recovering a high-dimensional matrix $\bm{A}_{\star}\in\mathbb{R}^{n_{1}\times n_{2}}$ from very few (noisy) observations of its entries is commonly known as matrix completion, whose applications include, collaborative filtering, computer visions and positioning. From a statistical viewpoint, it is a high-dimensional missing data problem where a high percentage of matrix entries are missing. As in many missing data problems, the underlying missing mechanism plays an important role. Most existing work (e.g., Candès and Recht, 2009; Keshavan et al., 2009; Recht, 2011; Rohde and Tsybakov, 2011; Koltchinskii et al., 2011) adopt a uniform observation mechanism, where each entry has the same marginal probability of being observed. This leads to significant simplifications, and enables the domain to move forward rapidly with various theoretical breakthroughs in the last decade. However, the uniform mechanism is often unrealistic. Recent works (Foygel et al., 2011; Negahban and Wainwright, 2012; Klopp, 2014; Cai and Zhou, 2016; Cai et al., 2016; Bi et al., 2017; Mao et al., 2019) have been devoted to relaxing such an restrictive assumption by adopting other missing structures. The usage of these settings hinges on strong prior knowledge of the underlying problems. At a high level, many of them utilize some special forms of low-rank structures for missing mechanism. For instance, Foygel et al. (2011) and Negahban and Wainwright (2012) both adopt a rank-1 structure based on the estimated marginal probabilities. In this paper, we aim at recovering the target matrix $\bm{A}_{\star}$ under a flexible high-dimensional low-rank sampling structure. This is achieved by a weighted empirical risk minimization, with application of inverse probability weighting (e.g., Schnabel et al., 2016; Mao et al., 2019) to adjust for the effect of non-uniform missingness.

Data arising in many applications of matrix completion, such as recommender systems, usually possesses complex “sampling” structure which is largely unknown. For a movie recommender system, users tend to rate movies that they prefer or dislike most, while often remain “silent” to other movies. Another example of the complex sampling regime is in the online merchandising, where some users may purchase certain items regularly without often rating them, but often evaluate products that they rarely buy. Similar to the widely adopted model that ratings are generated from a small number of hidden factors, it is reasonable to believe that the missingness is also governed by a small and possibly different set of hidden factors, which leads to a low-rank model the missing structure.

Inspired by generalized linear models, we model the probabilities of observation $\bm{\Theta}_{\star}=(\theta_{\star,ij})^{n_{1},n_{2}}_{i,j=1}\in(0,1)^{n_{1}\times n_{2}}$ by a high-dimensional low-rank matrix $\bm{M}_{\star}=(m_{\star,ij})_{i,j=1}^{n_{1},n_{2}}\in\mathbb{R}^{n_{1}\times n_{2}}$ through a known function $f$. That means, on the entry-wise level, we have $\theta_{\star,ij}=f(m_{\star,ij})$. In generalized linear models, the linear predictor $m_{\star,ij}$ is further modeled as a linear function of observed covariates. However, to reflect difficulties to attain (appropriate and adequate) covariate information and the complexity in the modeling of $\bm{\Theta}_{\star}$ in some situations of the matrix completion, the predictor matrix $\bm{M}_{\star}$ is assumed completely hidden in this study. Despite $\bm{M}_{\star}$ being hidden, as demonstrated in this work, the low-rankness of $\bm{M}_{\star}$ together with the high dimensionality of $\bm{M}_{\star}$ allows both identification and consistent estimation of $\bm{\Theta}_{\star}$, which facilitates inverse probability weighting based matrix completion. Motivated by the nature of matrix completion, we propose a novel parametrization $\bm{M}_{\star}={\mu}_{\star}\bm{1}_{n_{1}}\bm{1}_{n_{2}}^{{\mathrm{\scriptscriptstyle T}}}+\bm{Z}_{\star}$ where $\bm{Z}_{\star}$ satisfies $\bm{1}_{n_{1}}^{{\mathrm{\scriptscriptstyle T}}}\bm{Z}_{\star}\bm{1}_{n_{2}}=0$. Our proposal extends the work of Davenport et al. (2014), which aims to solve a binary matrix completion problem and pursues a different goal. Compared with Davenport et al. (2014), the proposed method does not regularize the estimation of $\mu_{\star}$, but only regularize the nuclear norm of the estimation of $\bm{Z}_{\star}$, which require different algorithmic treatment to avoid bias caused by the nuclear-norm penalty.

There are three fundamental aspects that set our work aside from the existing works of matrix completion as we consider: (i) the high-dimensional probability matrix $\bm{\Theta}$, whose dimensions $n_{1},n_{2}$ go to infinity in our asymptotic setting; (ii) the diminishing lower bound of the observation probabilities (as $n_{1},n_{2}$ go to infinity), and added issue to the instability of inverse probability weighting; (iii) the effects of estimation error in inverse probability weighting to the matrix completion procedure. Aspects (i) and (ii) are unique to our problem, and not found in the literature of missing data. The work related to Aspect (iii) is sparse in the literature of matrix completion. It is noted that Mao et al. (2019) focused on a low-dimensional parametric modeling of inverse probability weighting with observable covariates.

We develop non-asymptotic upper bounds of the mean squared errors for the proposed estimators of the observation probabilities and the target matrix. To sustain the convergence rate of the target matrix under the high-dimensionality of $\bm{M}_{\star}$ and low levels of observation probabilities, we propose to re-estimate $\bm{Z}_{\star}$ by constraining the magnitude of its entries to a smaller threshold. Our analysis shows that the proposed constrained inverse probability weighting estimator achieves the optimal rate (up to a logarithmic factor in estimation of target matrix). We also compare the inverse probability weighting based completion based on the proposed constrained estimation, with the completion based on direct weight trimming (or winsorization), a known practice in the conventional missing value literature (e.g., Rubin, 2001; Kang and Schafer, 2007; Schafer and Kang, 2008) and show that the constrained estimation has both theoretical and empirical advantages.

2 Model and Method

2.1 General Setup

Let $\bm{A}_{\star}=(a_{\star,ij})_{i,j=1}^{n_{1},n_{2}}\in\mathbb{R}^{n_{1}\times n_{2}}$ be an unknown high-dimensional matrix of interest, and $\bm{Y}=(y_{ij})_{i,j=1}^{n_{1},n_{2}}$ be a contaminated version of $\bm{A}_{\star}$ according to the following additive noise model:

[TABLE]

where $\{\epsilon_{ij}\}$ are independently distributed random errors with zero mean and finite variance. In the setting of matrix completion, only a portion of $\{y_{ij}\}$ is observed. For the $(i,j)$th entry, define the sampling indicator $w_{ij}=1$ if $y_{ij}$ is observed, and [math] otherwise, and assume $\{\epsilon_{ij}\}$ are independent of $\{w_{ij}\}$.

As for the sampling mechanism, we adopt a Bernoulli model where $\{w_{ij}\}$ are independent Bernoulli random variables with observation probabilities $\{\theta_{\star,ij}\}$, collectively denoted by a matrix $\bm{\Theta}_{\star}=(\theta_{\star,ij})_{i,j=1}^{n_{1},n_{2}}\in(0,1)^{n_{1}\times n_{2}}$. Similar to generalized linear models, the observation probabilities can be expressed in terms of an unknown matrix $\bm{M}_{\star}=(m_{\star,ij})_{i,j=1}^{n_{1},n_{2}}\in\mathbb{R}^{n_{1}\times n_{2}}$ and a pre-specified monotone and differentiable function $f:\mathbb{R}\to[0,1]$, i.e., $\theta_{\star,ij}=f(m_{\star,ij})$ for all $i,j$. The matrix $\bm{M}_{\star}$ plays the same role as a linear predictor in the generalized linear model, while the function $f$ is an inverse link function. Two popular choices of $f$ are inverse logit function $g(m)={e^{m}}/{(1+e^{m})}$ (logistic model) and the standard normal cumulative distribution function (probit model).

2.2 Low-rank Modeling of $\bm{A}_{\star}$ and $\bm{M}_{\star}$

The above setup is general. Without additional assumption, it is virtually impossible to recover the hidden feature matrix $\bm{M}_{\star}$ and also the target matrix $\bm{A}_{\star}$. A common and powerful assumption is that $\bm{A}_{\star}$ is a low-rank matrix, i.e., $\mathrm{rank}(\bm{A}_{\star})\ll\min\{n_{1},n_{2}\}$. Take the Yahoo! Webscope data set (to be analyzed in Section 7) as an example. This data set contains a partially observed matrix of ratings from 15,400 users to 1000 songs, and the goal is to complete the rating matrix. The low-rank assumption reflects the belief that users’ ratings are generated by a small number of factors, representing several standard preference profiles for songs. This viewpoint has been proven useful in the modeling of recommender systems (e.g., Candès and Plan, 2010; Cai et al., 2010).

The same idea could be adapted to the missing pattern, despite that the factors that induce the missingness may be different from those that generate the ratings. To this end, we assume $\bm{M}_{\star}$ is also low-rank. Next, we consider decomposing $\bm{M}_{\star}$ as

[TABLE]

with $\bm{1}_{n}$ being an $n$-vector of ones, and $\bm{J}=\bm{1}_{n_{1}}\bm{1}_{n_{2}}^{{\mathrm{\scriptscriptstyle T}}}$. Here $\mu_{\star}$ is the mean of $\bm{M}_{\star}$, i.e., $\mu_{\star}=\bm{1}_{n_{1}}^{{\mathrm{\scriptscriptstyle T}}}\bm{M}_{\star}\bm{1}_{n_{2}}/(n_{1}n_{2})$. Note that this decomposition holds for any matrix $\bm{M}$ by setting $\mu=(n_{1}n_{2})^{-1}\bm{1}_{n_{1}}^{\rm T}\bm{M}\bm{1}_{n_{2}}$ and $\bm{Z}=\bm{M}-\mu\bm{J}$. Moreover, the decomposition is unique due to the constraint that $\bm{1}_{n_{1}}^{{\mathrm{\scriptscriptstyle T}}}\bm{Z}_{\star}\bm{1}_{n_{2}}=0$. The key here is to reparametrize $\bm{M}_{\star}$ in terms of $\mu_{\star}$ and $\bm{Z}_{\star}$, which require different treatments in their estimation. See Section 3 for details. Further, the low-rankness of $\bm{M}_{\star}$ can be translated to the low-rankness of $\bm{Z}_{\star}$.

We note that the rank of $\bm{M}_{\star}$ is not the same as that of $\bm{\Theta}_{\star}$ due to the nonlinear transformation $f$. In general, the low-rank structure of $\bm{M}_{\star}$ implies a specific low-dimensional nonlinear structure of $\bm{\Theta}_{\star}$. For a common high missingness scenario, most entries of $\bm{M}_{\star}$ are significantly negative, where many common choices of the inverse link function can be well-approximated by a linear function. So our modeling can be regarded as a low-rank modeling of $\bm{\Theta}_{\star}$.

There are some related but more specialized models. Srebro and Salakhutdinov (2010) and Negahban and Wainwright (2012) utilize an independent row and column sampling mechanism, leading to a rank-1 structure for $\bm{\Theta}_{\star}$. Cai et al. (2016) consider a matrix block structure for $\bm{\Theta}_{\star}$ and hence $\bm{M}_{\star}$, which can be regarded as a special case of the low-rank modeling. Mao et al. (2019) considered the case when the missingness is depends on observable covariates, and adopted a low-rank modeling with a known row space of $\bm{M}_{\star}$. The proposal in this paper is for the situation when the missingness is dependent on some hidden factors, which reflects situations when obvious covariates are unknown or not available.

2.3 Inverse Probability Weighting Based Matrix Completion: Motivations and Challenges

Write the Hadamard product as $\circ$ and the Frobenius norm as $\|\cdot\|_{F}$. To recover the target matrix $\bm{A}_{\star}$, many existing matrix completion techniques assume uniform missing structure and hence utilize an unweighted/uniform empirical risk function $\widehat{R}_{\text{UNI}}(\bm{A})=(n_{1}n_{2})^{-1}\|\bm{W}\circ(\bm{A}-\bm{Y})\|_{F}^{2}$ (e.g., Candès and Plan, 2010; Koltchinskii et al., 2011; Mazumder et al., 2010), which is an unbiased estimator of the risk $R(\bm{A})=\textsf{E}(\|\bm{A}-\bm{Y}\|^{2}_{F})/(n_{1}n_{2})$ (up to a multiplicative constant) under uniform missingness. The work of Klopp (2014) is a notable exception that considers the use of $\widehat{R}_{\text{UNI}}$ under non-uniform missingness.

For any matrix $\bm{B}=(b_{ij})_{i,j=1}^{n_{1},n_{2}}$, we denote $\bm{B}^{\dagger}=(b_{ij}^{-1})_{i,j=1}^{n_{1},n_{2}}$ and $\bm{B}^{\ddagger}=(b_{ij}^{-1/2})_{i,j=1}^{n_{1},n_{2}}$. Under general missingness (uniform or non-uniform), one can show that, for any $\bm{A}\in\mathbb{R}^{n_{1}\times n_{2}}$,

[TABLE]

Clearly, $\bm{A}_{\star}$ uniquely minimizes $R$. If $\bm{\Theta}_{\star}$ were known, an unbiased estimator of $R$ would be

[TABLE]

which motivates the use of inverse probability weighting in matrix completion as in Mao et al. (2019). In addition, our theoretical analysis shows that the nuclear-norm-regularized empirical risk estimator (to be defined in details later) based on $\widehat{R}$ (assuming the use of true observation probabilities) improves upon existing error upper bound of corresponding estimator based on $\widehat{R}_{\text{UNI}}$ achieved by Klopp (2014) as shown in Section 5.3. However, the inverse probability weights $\bm{\Theta}^{\ddagger}_{\star}$ are often unknown and have to be estimated, which has to be conducted carefully in the context of matrix completion.

Despite the popularity of inverse probability weighting in missing data literature, it is known to produce unstable estimation due to occurrences of small probabilities (e.g., Rubin, 2001; Kang and Schafer, 2007; Schafer and Kang, 2008). This problematic scenario is common for matrix completion problems for recovering a target matrix from very few observations. Theoretically, a reasonable setup should allow some $\theta_{\star,ij}$ to go to zero as $n_{1},n_{2}\rightarrow\infty$, leading to diverging weights and a non-standard setup of inverse probability weighting. Due to these observations, a careful construction of the estimation procedure is required.

For uniform sampling ($\theta_{\star,ij}\equiv\theta_{0}$ for some probability $\theta_{0}$), one only has to worry about a small common probability $\theta_{0}$ (or that $\theta_{0}$ diminishes in an asymptotic sense.) Although small $\theta_{0}$ increases the difficulty of estimation, $\widehat{R}(\bm{A})$ changes only up to a multiplicative constant. However, for non-uniform setting, it is not as straightforward due to the heterogeneity among $\{\theta_{\star,ij}\}$. To demonstrate the issue, we now briefly look at the Yahoo! Webscope dataset described in Section 7. A sign of the strong heterogeneity in $\{\theta_{\star,ij}\}$ is a large $\theta_{U}/\theta_{L}$, where $\theta_{L}=\min_{i,j}\theta_{\star,ij}$ and $\theta_{U}=\max_{i,j}\theta_{\star,ij}$. We found that the corresponding ratio of estimated probabilities $\widehat{\theta}_{U}/\widehat{\theta}_{L}$ based on the rank-1 structure of Negahban and Wainwright (2012) was $25656.2$, and that based on our proposed method (without re-estimation, to be described below) was $23988.0$. This strong heterogeneity can jeopardize the convergence rate of our estimator, which will be properly addressed in our framework.

In the following section, we propose an estimation approach for $\bm{\Theta}_{\star}$ in Section 3.1 and an appropriate modification in Section 3.3 which, when substituted into the empirical risk $\widehat{R}$, allows us to construct a stable estimator for $\bm{A}_{\star}$.

3 Estimation of $\bm{\Theta}_{\star}$

3.1 Regularized Maximum Likelihood Estimation

We develop the estimation of $\bm{\Theta}_{\star}$ based upon the framework of regularized maximum likelihood. Given the inverse of link function $f$, the log-likelihood function with respect to the indicator matrix $\bm{W}=(w_{ij})\in\mathbb{R}^{n_{1}\times n_{2}}$ is

[TABLE]

for any $\bm{M}=(m_{ij})_{i,j=1}^{n_{1},n_{2}}\in\mathbb{R}^{n_{1}\times n_{2}}$, where $\mathds{I}_{\mathcal{A}}$ is the indicator of an event $\mathcal{A}$. Due to the low-rank assumption of $\bm{M}_{\star}$, one natural candidate of estimators to $\bm{M}_{\star}$ is the maximizer of the regularized log-likelihood $\ell_{\bm{W}}(\bm{M})-\lambda\|\bm{M}\|_{*}$, where $\|\cdot\|_{*}$ represents the nuclear norm and $\lambda>0$ is a tuning parameter. It is also common to enforce an additional max-norm constraint $\|\bm{M}\|_{\infty}\leq\alpha$ for some $\alpha>0$ in the maximization (e.g., Davenport et al., 2014). Note that the nuclear norm penalty flavors $\bm{M}=0$, corresponding to that $\Pr(w_{ij}=1)=0.5$ for all $i,j$. Nevertheless, this would not align well with common settings of matrix completion under which the average probability of observations is quit small, and hence would result in a large bias. In view of this, we instead adopt a parametrization $\bm{M}_{\star}=\mu_{\star}\bm{J}+\bm{Z}_{\star}$ and consider the following estimator of $(\mu_{\star},\bm{Z}_{\star})$:

[TABLE]

[TABLE]

Note that the mean $\mu$ of the linear predictor $\mu\bm{J}+\bm{Z}$ is not penalized. The constraint $\bm{1}_{n_{1}}^{{\mathrm{\scriptscriptstyle T}}}\bm{Z}\bm{1}_{n_{2}}=0$ ensures the identifiability of $\mu$ and $\bm{Z}$. Apparently, the constraints in $\mathcal{C}_{n_{1},n_{2}}\left(\alpha_{1},\alpha_{2}\right)$ are analogous to $\|\bm{M}\|_{\infty}\leq\alpha_{0}$, where $\alpha_{0}=\alpha_{1}+\alpha_{2}$, but on the parameters $\mu$ and $\bm{Z}$ respectively. With $(\widehat{\mu},\widehat{\bm{Z}})$, the corresponding estimator of $\bm{M}_{\star}$ is $\widehat{\bm{M}}=\widehat{\mu}\bm{J}+\widehat{\bm{Z}}$.

Davenport et al. (2014) considered a regularized maximum likelihood approach for a binary matrix completion problem. Their goal was different, as they aimed at recovering a binary rating matrix in lieu of the missing structure, and considered a regularization on $\bm{M}$ (instead of $\bm{Z}$) via $\|\bm{M}\|_{\ast}\leq\alpha^{\prime}\{\mathrm{rank}(\bm{M}_{\star})n_{1}n_{2}\}^{1/2}$. As for the scaling parameter $\alpha^{\prime}$, Davenport et al. (2014) considered an $\alpha^{\prime}$ independent of the dimensions $n_{1}$ and $n_{2}$ to restrict the “spikiness” of $\bm{M}$. As explained earlier, in our framework, $\theta_{L}$ should be allowed to go to zero as $n_{1},n_{2}\rightarrow\infty$. To this end, we allow $\alpha_{1}$ and $\alpha_{2}$ to depend on the dimensions $n_{1}$ and $n_{2}$. See more details in Section 5.

3.2 Computational algorithm and tuning parameter selection

To solve the optimization (1), we begin with the observation that $\ell_{\bm{W}}$ is a smooth concave function, which allows the usage of an iterative algorithm called accelerated proximal gradient algorithm (Beck and Teboulle, 2009). Given a pair $(\mu_{\text{old}},\bm{Z}_{\text{old}})$ from a previous iteration, a quadratic approximation of the objective function $-\ell_{\bm{W}}(\mu\bm{J}+\bm{Z})+\lambda\|\bm{Z}\|_{\ast}$ is formed:

[TABLE]

where $L>0$ is an algorithmic parameter determining the step size of the proximal gradient algorithm, and is chosen by a backtracking method (Beck and Teboulle, 2009). Here $\langle\bm{B},\bm{C}\rangle=\sum_{i,j}b_{ij}c_{ij}$ for any matrices $\bm{B}=(b_{ij})$ and $\bm{C}=(c_{ij})$ of same dimensions.

In this iterative algorithm, a successive update of $(\mu,\bm{Z})$ can be obtained by

[TABLE]

where the optimization with respect to $\mu$ and $\bm{Z}$ can be performed separately. For $\mu$, one can derive a closed-form update

[TABLE]

As for $\bm{Z}$, we need to perform the minimization

[TABLE]

which is equivalent to

[TABLE]

We apply a three-block extension of the alternative direction method of multipliers (Chen et al., 2016) to an equivalent form of (2):

[TABLE]

Write $\bm{H}=(\bm{H}_{1},\bm{H}_{2})$. The augmented Lagrangian for (3) is

[TABLE]

where $u>0$ is an algorithmic parameter and, $\mathds{I}_{\mathcal{A}}=0$ if the constraint $\mathcal{A}$ holds and $\infty$ otherwise. The detailed algorithm to solve (3) is summarized in Algorithm LABEL:appalg:ADMM in the supplementary material. It is noted that, in general, the multi-block alternative direction method of multipliers may fail to converge for some $u>0$ (Chen et al., 2016). In those cases, an appropriate selection of $u$ is crucial. However, we are able to show that the form of our algorithm belongs to a special class (Chen et al., 2016) in which convergence is guaranteed for any $u>0$. Therefore, we simply set $u=1$. We summarize the corresponding convergence result in the following theorem whose proof is provided in the supplementary material.

Theorem 1**.**

The sequence $\{\bm{Z}^{(k)},\bm{G}_{1}^{(k)},\bm{G}_{2}^{(k)}\}$, generated by Algorithm LABEL:appalg:ADMM in the supplementary material, converges to the global optimum of (3).

Notice that the alternative direction method of multipliers algorithm is nested within the proximal gradient algorithm. But, from our practical experiences, both the number of inner iterations (alternative direction method of multipliers) and outer iterations (proximal gradient) are small, usually less than twenty in our numerical experiments. Similarly, we summarize the corresponding convergence result of the overall proximal gradient algorithm in the following theorem whose proof is provided in the supplementary material.

Theorem 2**.**

The estimator $(\widehat{\mu},\widehat{\bm{Z}})$ generated by the proximal gradient algorithm, converges to the global optimum of (1).

The tuning parameters $\alpha_{1}$ and $\alpha_{2}$ can be chosen according to prior knowledge of the problem setup, if available. When a prior knowledge is not available, one can choose large values for these parameters. Once these parameters are large enough, our method is not sensitive to their specific values. A more principled way to tune $\alpha_{1}$ and $\alpha_{2}$ is a challenging problem and beyond the scope of this work. As for $\lambda$, we adopt Akaike information criterion (aic) where the degree of freedom is approximated by $r_{\widehat{\bm{M}}}(n_{1}+n_{2}-r_{\widehat{\bm{M}}})$.

3.3 Constrained estimation

To use $\widehat{R}$ of (3), a naive idea is to obtain $\widehat{\bm{\Theta}}=(\hat{\theta}_{ij})_{i,j=1}^{n_{1},n_{2}}=\mathcal{F}({\widehat{\bm{M}}})$, where $\mathcal{F}(\bm{M})=(f(m_{ij}))_{i,j=1}^{n_{1},n_{2}}\in\mathbb{R}^{n_{1}\times n_{2}}$ for any $\bm{M}=(m_{ij})_{i,j=1}^{n_{1},n_{2}}\in\mathbb{R}^{n_{1}\times n_{2}}$, and then replace $\bm{\Theta}^{\ddagger}_{\star}$ by $\widehat{\bm{\Theta}}^{\ddagger}=(\hat{\theta}_{ij}^{-1/2})_{i,j=1}^{n_{1},n_{2}}$. However, this direct implementation is not robust to extremely small probabilities of observation, and our theoretical analysis shows that this could lead to a slower convergence rate in the estimation of $\bm{A}_{\star}$. In the literature of missing data, a simple solution to robustify is winsorizing the small probabilities (Potter, 1990; Scharfstein et al., 1999).

In the estimation of $\widehat{\bm{\Theta}}$ defined in (1) that assumes $\|\bm{Z}_{\star}\|_{\infty}\leq\alpha_{2}$, a large $\alpha_{2}$ has an adverse effect on the estimation. In the setting of diverging $\alpha_{2}$ (due to diminishing $\theta_{L}$), the convergence rate of $\widehat{\bm{Z}}$ becomes slower and the estimator obtained after direct winsorization will also be affected. That is, even though the extreme probabilities could be controlled by winsorizing, the unchanged entries of $\widehat{\bm{Z}}$ (in the procedure of winsorizing) may already suffer from a slower rate of convergence. This results in a larger estimation error under certain settings of missingness, as revealed in Section 5.

A seemingly better strategy is to impose a tighter constraint directly in the minimization problem (1). That is to adopt the constraint $\|\bm{Z}\|_{\infty}\leq\beta$ where $0\leq\beta\leq\alpha_{2}$. Theoretically, one can better control the errors on those entries of magnitude smaller than $\beta$. However, the mean-zero constraint of $\bm{Z}$ no longer makes sense as the constraint $\|\bm{Z}\|_{\infty}\leq\beta$ may have shifted the mean.

We propose a re-estimation of $\bm{Z}_{\star}$ with a different constraint level $\beta$:

[TABLE]

Note that we only re-compute $\bm{Z}$ but not $\mu$, which allows us to drop the mean-zero constraint. Thus, $\widehat{\bm{M}}_{\beta}=\widehat{\mu}\bm{J}+\widehat{\bm{Z}}_{\beta}$. The corresponding algorithm for optimization (4) can be derived similarly as in Davenport et al. (2014), and is provided in the supplementary material. In what follows, we write $\widehat{\bm{\Theta}}=\mathcal{F}(\widehat{\bm{M}})$ and $\widehat{\bm{\Theta}}_{\beta}=\mathcal{F}(\widehat{\bm{M}}_{\beta})$.

4 Estimation of $\bm{A}_{\star}$

Now, we come back to (3) and replace $\bm{\Theta}^{\ddagger}_{\star}$ by $\widehat{\bm{\Theta}}^{\ddagger}_{\beta}$ to obtain a modified empirical risk:

[TABLE]

where $\widehat{\bm{\Theta}}^{\ddagger}_{\beta}=(\widehat{\theta}_{ij,\beta}^{-1/2})\in\mathbb{R}^{n_{1}\times n_{2}}$. Since $\bm{A}$ is a high-dimensional parameter, a direct minimization of $\hat{R}^{\ast}$ often results in over-fitting. To circumvent it, we consider a regularized version:

[TABLE]

where $\tau>0$ is a regularization parameter. Again, the nuclear norm regularization encourages low-rank solution. Based on (2), our estimator of $\bm{A}_{\star}$ is

[TABLE]

where $a$ is an upper bound on $\|\bm{A}_{\star}\|_{\infty}$. The above $\widehat{\bm{A}}_{\beta}$ contains as special cases (i) the matrix completion $\widehat{\bm{A}}_{\alpha_{2}}$, with unconstrained probability estimator $\widehat{\bm{\Theta}}$, by setting $\beta=\alpha_{2}$ and (ii) the estimator $\widehat{\bm{A}}_{\beta}$, with constrained probability estimator $\widehat{\bm{\Theta}}_{\beta}$, when $\beta<\alpha_{2}$.

We use an accelerated proximal gradient algorithm (Beck and Teboulle, 2009) to solve (3). For the choice of tuning parameter $\tau$ in (3), we adopt a $5$-fold cross-validation to select the remaining tuning parameters. Due to the non-uniform missing mechanism, we use a weighted version of the validation errors. The specific details are given in Algorithm LABEL:appalg:A in the supplementary material.

5 Theoretical Properties

5.1 Probabilities of observation

Let $\|\bm{B}\|=\sigma_{\max}(\bm{B})$, $\|\bm{B}\|_{\infty}=\max_{i,j}|b_{ij}|$ and $\|\bm{B}\|_{\infty,2}=(\max_{i}\sum_{j}b_{ij}^{2})^{1/2}$ be the spectral norm, the maximum norm and $l_{\infty,2}$-norm of a matrix $\bm{B}$ respectively. We use the symbol $\asymp$ to represent the asymptotic equivalence in order, i.e., $a_{n}\asymp b_{n}$ if $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. The average squared distance between two matrices $\bm{B},\bm{C}\in\mathbb{R}^{n_{1}\times n_{2}}$ is $d^{2}(\bm{B},\bm{C})=\|\bm{B}-\bm{C}\|_{F}^{2}/(n_{1}n_{2})$. The average squared errors of $\widehat{\bm{M}}_{\beta}$ and $\widehat{\bm{\Theta}}^{\dagger}_{\beta}$ are then $d^{2}(\widehat{\bm{M}}_{\beta},\bm{M}_{\star})$ and $d^{2}(\widehat{\bm{\Theta}}^{\dagger}_{\beta},\bm{\Theta}^{\dagger}_{\star})$ respectively. We adopt the Hellinger distance for any two matrices $S,T\in[0,1]^{n_{1}\times n_{2}}$, $d_{H}^{2}(S,T)=(n_{1}n_{2})^{-1}\sum_{i,j=1}^{n_{1},n_{2}}d_{H}^{2}(s_{ij},t_{ij})$ where $d_{H}^{2}(s,t)=(s^{1/2}-t^{1/2})^{2}+\{(1-s)^{1/2}-(1-t)^{1/2}\}^{2}$ for $s,t\in[0,1]$. In the literature of matrix completion, most discussions related to optimal convergence rate are only up to certain polynomial orders of $\log n$. For convenience, we use $\mathrm{polylog}(n)$ for polynomials of $\log n$.

To investigate the asymptotic properties of $\widehat{\bm{M}}_{\beta}$ and $\widehat{\bm{\Theta}}^{\dagger}_{\beta}$ defined in Section 3, we introduce the following conditions on the missing structure.

C1. The indicators $\{w_{ij}\}_{i,j=1}^{n_{1},n_{2}}$ are mutually independent, and independent of $\{\epsilon_{ij}\}_{i,j=1}^{n_{1},n_{2}}$. For $i=1,\dots,n_{1}$ and $j=1,\dots,n_{2}$, $w_{ij}$ follows a Bernoulli distribution with probability of success $\theta_{\star,ij}=f(m_{\star,ij})\in(0,1)$. Furthermore, $f$ is monotonic increasing and differentiable.

C2. The hidden feature matrix $\bm{M}_{\star}=\mu_{\star}\bm{J}+\bm{Z}_{\star}$ where $\bm{1}_{n_{1}}^{T}\bm{Z}_{\star}\bm{1}_{n_{2}}=0$, $|\mu_{\star}|\leq\alpha_{1}<\infty$ and $\|\bm{Z}_{\star}\|_{\infty}\leq\alpha_{2}<\infty$. Here $\alpha_{1}$ and $\alpha_{2}$ are allowed to depend on the dimensions $n_{1}$ and $n_{2}$. This also implies that there exists a lower bound $\theta_{L}\in(0,1)$ (allowed to depend on $n_{1},n_{2}$) such that $\underset{i,j}{\min}\{\theta_{ij}\}\geq\theta_{L}\geq f(-\alpha_{1}-\alpha_{2})>0$.

For the convenience of the theoretical analysis, we consider an equivalent estimator of $(\mu_{\star},\bm{Z}_{\star})$ defined by the constrained maximization problem (1) instead of the Lagrangian form (1). For $r_{\bm{Z}_{\star}}\leq\min\{n_{1},n_{2}\}$ and $\alpha_{1},\alpha_{2}\geq 0$, let

[TABLE]

[TABLE]

It is easy to see that we have $(\mu_{\star},\bm{Z}_{\star})\in\widetilde{\mathcal{C}}_{n_{1},n_{2}}(\alpha_{1},\alpha_{2})$ once $(\mu_{\star},\bm{Z}_{\star})\in\mathcal{C}_{n_{1},n_{2}}(\alpha_{1},\alpha_{2})$ holds. For the ease of presentation, we assume $n_{1}=n_{2}=n$ and choose the logit function as the inverse link function $f$ in the rest of Section 5, while corresponding results under general settings of $n_{1}$, $n_{2}$ and $f$ are delegated to Section LABEL:appsec:gentheory in the supplementary material. We first establish the convergence results for $\widehat{\mu}$, $\widehat{\bm{Z}}$ and $\widehat{\bm{M}}$, respectively. To simplify notations, let $\alpha_{0}=\alpha_{1}+\alpha_{2}$, $h_{\alpha_{1},\beta}=(1+e^{\alpha_{1}+\beta})^{-1}$ and $\Gamma_{n}=e^{\alpha_{0}}(\alpha_{1}+\alpha_{2}r^{1/2}_{\bm{Z}_{\star}})n^{-1/2}$.

Lemma 1**.**

Suppose Conditions C1-C2 hold, and $(\mu_{\star},\bm{Z}_{\star})\in\mathcal{C}_{n_{1},n_{2}}(\alpha_{1},\alpha_{2})$. Consider $\widehat{\bm{M}}=\widehat{\mu}\bm{J}+\widehat{\bm{Z}}$ where $(\widehat{\mu},\widehat{\bm{Z}})$ is the solution to (1). There exist positive constants $C_{1},C_{2}$ such that we have with probability at least $1-C_{1}/n$,

[TABLE]

The upper bounds in (1) all consist of trivial bounds $\alpha_{j}^{2}$ and a more dedicated bound $\Gamma_{n}$. The trivial upper bounds $\alpha_{1}^{2}$, $\alpha_{2}^{2}$ and $\alpha_{0}^{2}$ can be easily derived from the constraint set $\mathcal{C}_{n_{1},n_{2}}(\alpha_{1},\alpha_{2})$. For extreme settings of increasing $\alpha_{0}$, the more dedicated bound $\Gamma_{n}$ is diverging and the trivial bounds may provide better control. The term $\Gamma_{n}$ can be controlled by the rank of $\bm{Z}_{\star}$. For a range of non-extreme scenarios, i.e., $\alpha_{0}\leq 1/2\log n$ or $\theta_{L}\geq n^{-1/2}$, the second term in $\Gamma_{n}$ attains the convergence order once $r_{\bm{Z}_{\star}}=O(1)$.

Similarly, we study the theoretical results of the re-estimation of $\bm{Z}_{\star}$ in terms of the constrained optimization:

[TABLE]

We now consider the constrained estimation for $\bm{Z}_{\star}$, $\bm{M}_{\star}$ and $\bm{\Theta}^{\dagger}_{\star}$. For any matrix $\bm{B}=(b_{ij})_{i,j=1}^{n_{1},n_{2}}$, define the winsorizing operator $\mathcal{T}_{\beta}$ by $\mathcal{T}_{\beta}(\bm{B})=(T_{\beta}(b_{ij}))$ where

[TABLE]

Write $\bm{M}_{\star,\beta}=\mu_{\star}\bm{J}+\mathcal{T}_{\beta}(\bm{Z}_{\star})$ and $\widehat{\bm{M}}_{\star,\beta}=\widehat{\mu}\bm{J}+\mathcal{T}_{\beta}(\bm{Z}_{\star})$, and $\bm{\Theta}_{\star,\beta}=\mathcal{F}(\bm{M}_{\star,\beta})$ and $\widehat{\bm{\Theta}}_{\star,\beta}=\mathcal{F}(\widehat{\bm{M}}_{\star,\beta})$ respectively. It is noted that $\widehat{\bm{M}}_{\star,\beta}$ serves as a “bridge” between the underlying $\bm{M}_{\star,\beta}$ and the empirical $\widehat{\bm{M}}_{\beta}$. Write $N_{\beta}=\sum_{i,j}(\mathds{I}_{[z_{\star,ij}>\beta]}+\mathds{I}_{[z_{\star,ij}<-\beta]})$ as the number of extreme values in $\bm{Z}_{\star}$ at level $\beta$. The convergence rates of $d^{2}(\widehat{\bm{M}}_{\beta},\bm{M}_{\star})$ and $d^{2}(\widehat{\bm{\Theta}}^{\dagger}_{\beta},\bm{\Theta}^{\dagger}_{\star})$ are investigated in the next theorem. Define $\Lambda_{n}=\min[\beta^{2},\widetilde{\Gamma}_{n}+h_{\alpha_{1},\beta}^{-1}n^{-2}\beta\{8N_{\beta}+(n^{2}-N_{\beta})|\mu_{\star}-\widehat{\mu}|\}]$ where $\widetilde{\Gamma}_{n}=h_{\alpha_{1},\beta}^{-1}(\alpha_{1}+\beta r^{1/2}_{\mathcal{T}_{\beta}(\bm{Z}_{\star})})n^{-1/2}$.

Theorem 3**.**

Assume that Conditions C1-C2 hold, and $(\mu_{\star},\bm{Z}_{\star})\in\widetilde{\mathcal{C}}_{n_{1},n_{2}}(\alpha_{1},\alpha_{2})$. Consider $\widehat{\bm{M}}_{\beta}=\widehat{\mu}\bm{J}+\widehat{\bm{Z}}_{\beta}$ where $\widehat{\bm{Z}}_{\beta}$ is the solution to (5.1) and $\beta\geq 0$, there exist some positive constants $C_{1}$, $C_{2}$ and $C_{3}$ such that we have with probability at least $1-2C_{1}/n$,

[TABLE]

[TABLE]

We can derive an upper bound $4\beta^{2}$ for $d^{2}(\widehat{\bm{Z}}_{\textsf{Win},\beta},\mathcal{T}_{\beta}(\bm{Z}_{\star}))$ from the second term in Theorem 1 where $\widehat{\bm{Z}}_{\textsf{Win},\beta}=\mathcal{T}_{\beta}(\widehat{\bm{Z}})$ is directly winsorized from $\widehat{\bm{Z}}$. Obviously, the order of this upper bound is larger than or equal to $\Lambda_{n}$. Moreover, there are scenarios where $\Lambda_{n}$ is a smaller order of $\beta^{2}$. To illustrate, assume that both $\alpha_{1}\asymp 1$ and $\beta\asymp 1$, we have $h_{\alpha_{1},\beta}\asymp 1$. Once we have $N_{\beta}=o(n)$, $r_{\mathcal{T}_{\beta}(\bm{Z}_{\star})}=o(n)$ and $|\widehat{\mu}-\mu_{\star}|=o(1)$, then $\Lambda_{n}=o(\beta^{2})$.

With a more dedicated investigation of (5), one can derive an upper bound for $d^{2}(\widehat{\bm{\Theta}}^{\dagger}_{\beta},\widehat{\bm{\Theta}}^{\dagger}_{\star,\beta})$, which will be used in Section 5.2. Denote $k_{\alpha_{1},\alpha_{2},n}^{\prime}=\min\{\alpha_{1}^{2},e^{\alpha_{0}}(\alpha_{1}+\alpha_{2}r^{1/2}_{\bm{Z}_{\star}})n^{-1/2}\}$, such an upper bound is of order $k_{\alpha_{1},\alpha_{2},\beta,n}h_{\alpha_{1},\beta}^{-2}$ where

[TABLE]

5.2 Target matrix

To study the convergence of $d^{2}(\widehat{\bm{A}}_{\beta},\bm{A}_{\star})$, we require the following conditions regarding the random errors $\bm{\epsilon}$ and the target matrix $\bm{A}_{\star}$. Recall that $\widehat{\bm{A}}_{\beta}$ includes both the estimations obtained with the unconstrained estimator $\widehat{\bm{\Theta}}$ and the constrained estimator $\widehat{\bm{\Theta}}_{\beta}$ as $\widehat{\bm{A}}(\widehat{\bm{\Theta}})=\widehat{\bm{A}}_{\alpha_{2}}$ with $\beta=\alpha_{2}$.

C3. (a) The random errors $\{\epsilon_{ij}\}$ in Model (1) are independently distributed random variables such that $\textsf{E}(\epsilon_{ij})=0$ and $\textsf{E}(\epsilon^{2}_{ij})=\sigma_{ij}^{2}<\infty$ for all $i,j$. (b) For some finite positive constants $c_{\sigma}$ and $\eta$, $\underset{i,j}{\max}\,\textsf{E}|\epsilon_{ij}|^{l}\leq\frac{1}{2}l!c_{\sigma}^{2}\eta^{l-2}$ for any positive integer $l\geq 2$.

C4. There exists a positive constant $a_{0}$ such that $\|\bm{A}_{\star}\|_{\infty}\leq a_{0}$.

Denote $h_{(1),\beta}=\underset{i,j}{\max}(\theta_{\star,ij}^{-1}\theta_{\star,ij,\beta})$ and

[TABLE]

The following theorem established a general upper bound for $d^{2}(\widehat{\bm{A}}_{\beta},\bm{A}_{\star})$.

Theorem 4**.**

Assume Conditions C1-C4 hold. For $\beta\geq 0$, there exist some positive constants $C_{4}$, $C_{5}$, $C_{6}$ and $C_{7}$, both independent of $\beta$, such that for $h_{(1),\beta}\tau\geq C_{4}\Delta$, we have with probability at least $1-3/(2n)$,

[TABLE]

As for the estimator of the target matrix based on direct winsorization $\widehat{\bm{\Theta}}_{\textsf{Win},\beta}=\mathcal{F}(\widehat{\mu}\bm{J}+\widehat{\bm{Z}}_{\textsf{Win},\beta})$ where $\widehat{\bm{Z}}_{\textsf{Win},\beta}=\mathcal{T}_{\beta}(\widehat{\bm{Z}})$, an upper bound can be derived using Theorem 3. As noted in a remark after Theorem 3, $d^{2}(\widehat{\bm{Z}}_{\textsf{Win},\beta},\mathcal{T}_{\beta}(\bm{Z}_{\star}))$ converges at a slower rate $\beta^{2}$ which will cause a larger error bound for the target matrix.

Now, we discuss the rates of $d^{2}(\widehat{\bm{A}}_{\beta},\bm{A}_{\star})$ under various missing structures. For simplicity, the following discussion focuses on the low-rank linear predictor ($\bm{M}_{\star}$) setting such that $r_{\bm{M}_{\star}}\asymp 1$.

Uniform missingness. Under the uniform missingness, i.e., $\theta_{ij}\equiv\theta_{0}$, it has been shown in Koltchinskii et al. (2011) that $\theta_{0}^{-1}n^{-1}\mathrm{polylog}(n)$ is the optimal rate for $d^{2}(\widehat{\bm{A}}_{\beta},\bm{A}_{\star})$. Therefore it is reasonable to require $\alpha_{1}+\alpha_{2}=\alpha_{0}=O(\mathrm{polylog}(n))$ for the convergence of $d^{2}(\widehat{\bm{A}}_{\beta},\bm{A}_{\star})$. Under the uniform missingness, we have $\alpha_{2}=0$, $\alpha_{0}=\alpha_{1}$ and $e^{\mu_{\star}}\asymp\theta_{0}$. For $\beta=0$, our estimator $\widehat{\bm{A}}_{\beta}$ degenerates to the estimator based on the unweighted empirical risk function. Theorem 4 shows that $\widehat{\bm{A}}_{\beta}$ achieves the optimal rate $\theta_{0}^{-1}n^{-1}\mathrm{polylog}(n)$. As for $\beta>0$, by taking $\beta\to 0$ such that $k_{\alpha_{1},\alpha_{2},\beta,n}=O(e^{\mu_{\star}-2\alpha_{1}-2\beta}n^{-1}\log^{-2}n)$, the estimator can also reach the optimal rate. Of interest here is that $\beta$ is allowed to be strictly positive to achieve the same rate.

Non-uniform missingness. Under the non-uniform missingness, suppose the lower and upper bounds of observation probability satisfy $\theta_{L}\asymp e^{\mu_{\star}-\alpha_{2}}$ and $\theta_{U}\asymp e^{\mu_{\star}+\alpha_{2}}$. For the non-constrained case of $\beta=\alpha_{2}$ and $h_{\alpha_{1},\beta}\asymp e^{-\alpha_{1}-\alpha_{2}}$, the second term of $\Delta$ in (6) dominates due to the fact that

[TABLE]

Thus, the convergence rate of $d^{2}(\widehat{\bm{A}}_{\beta},\bm{A}_{\star})$ is $e^{\mu_{\star}+5\alpha_{1}+3\alpha_{2}}n^{-1/2}\log^{3}n$. To guarantee convergence, as $e^{\mu_{\star}/2+5\alpha_{1}/2+3\alpha_{2}/2}\leq e^{3\alpha_{1}+3\alpha_{2}/2}$, it requires that $\alpha_{1}+\alpha_{2}/2<(1/12)\log n$ which implies that $\theta_{L}^{-1}=O(n^{1/6})$.

However, the above range of $\theta_{L}=O(n^{-1/6})$ excludes $\theta_{L}\equiv(n^{-1}\mathrm{polylog}(n))$, the case that results in the number of the observed matrix entries at the order of $n\,\mathrm{polylog}(n)$ which represents the most sparse case of observation where the matrix can still be recovered (Candès and Recht, 2009; Candès and Plan, 2010; Koltchinskii et al., 2011; Negahban and Wainwright, 2012). We will show in the following that with an appropriately chosen $\beta$, the constrained estimator $\widehat{\bm{\Theta}}_{\beta}$ can accommodate the case of $\theta_{L}^{-1}=O(n\log^{-1}n)$.

Case (I): $\beta=0$. To demonstrate this, we start with the absolute constrained case, i.e., $\beta=0$, which forces the estimated probabilities to be uniform and implies $e^{-\mu_{\star}/2+\alpha_{2}-\beta+|\alpha_{2}/2-\beta|}=e^{-\mu_{\star}/2+3\alpha_{2}/2}\asymp\theta_{U}^{1/2}\theta_{L}^{-1}$. Then, according to Theorem 4, $d^{2}(\widehat{\bm{A}}_{\beta},\bm{A}_{\star})$ attains the convergence rate $\theta_{U}\theta_{L}^{-2}n^{-1}\log(n)$, which converges to 0 provided $\theta_{U}\theta_{L}^{-2}=o(n\log^{-1}n)$. Obviously, the condition $\theta_{U}\theta_{L}^{-2}=o(n\log^{-1}n)$ includes the extreme case of $\theta_{L}^{-1}=O(n\log^{-1}n)$ and $n\,\mathrm{polylog}(n)$ observations.

Case (II): $\beta>0$. For the more interesting setting $\beta>0$, to simplify the discussion, we concentrate on the case when the first term in $k_{\alpha_{1},\alpha_{2},\beta,n}$ is of a smaller order, which can be achieved by choosing $\beta=O(e^{-\mu_{\star}-2\alpha_{1}+\alpha_{2}}n^{-1/2}\log^{-1}n)$. Then, according to Theorem 4,

[TABLE]

since $e^{-\mu_{\star}/2+\alpha_{2}-\beta+|\alpha_{2}/2-\beta|}\leq e^{\alpha_{1}/2+3\alpha_{2}/2}$. In the following we consider two further cases: (i) $\alpha_{2}=O((\log\log n)^{-1}\alpha_{1})$ and (ii) $\alpha_{1}=o(\alpha_{2}\log\log n)$. Note that for either cases, $e^{-\mu_{\star}+2\alpha_{2}-2\beta+2|\alpha_{2}/2-\beta|}\asymp\theta_{U}\theta_{L}^{-2}$ which leads to

[TABLE]

If $\alpha_{2}=O\{(\log\log n)^{-1}\alpha_{1}\}$, we require $\alpha_{1}<(1+3\log\log n)^{-1}(\log n-\log\log n)$ to guarantee convergence, which implies that $\theta_{L}=O(n^{-1})$. Thus, we only lose a $\mathrm{polylog}(n)$ factor when compared with the most extreme but feasible setting of $\theta_{L}^{-1}=O[n\{\mathrm{polylog}(n)\}^{-1}]$. Also $\beta=O(e^{-\mu_{\star}-2\alpha_{1}+\alpha_{2}}n^{-1/2}\log^{-1}n)$ implies that $\beta=O(n^{-1/2}\log^{-1}n)$.

If $\alpha_{1}=o\{(\log\log n)\alpha_{2}\}$, we require that $\alpha_{2}<\{3+(\log\log n)^{-1}\}^{-1}(\log n-\log\log n)$ which leads to $\theta_{L}^{-1}=O(n^{1/3})$. Also $\beta=O(e^{-\mu_{\star}-2\alpha_{1}+\alpha_{2}}n^{-1/2}\log^{-1}n)$ implies that $\beta=O(n^{-1/6}\log^{-1}n)$. However, to make $d^{2}(\widehat{\bm{A}}_{\beta},\bm{A}_{\star})$ convergent, the attained rate for $\theta_{L}^{-1}$ has to be $O(n^{1/3})$, which excludes the most extreme heterogeneity case of $\theta_{L}^{-1}=O\{n(\mathrm{polylog}(n))^{-1}\}$. The reason for not being able to cover the most extreme case of $\theta_{L}^{-1}=O\{n(\mathrm{polylog}(n))^{-1}\}$ is that the current Case (ii) allows more heterogeneity in $\bm{Z}_{\star}$ as reflected by having a larger $\alpha_{2}$ than that prescribed under Case (i). As $\mu_{\star}$ is jointly estimated with $\bm{Z}_{\star}$ in the unconstrained estimation (Section 3.1), stronger heterogeneity slows down the convergence rate in the estimation of $\mu_{\star}$, which becomes a bottleneck for further improvement. If $\mu_{\star}$ was observable, the problem would not be as serious despite the adverse effect of stronger heterogeneity on the estimation of $\bm{Z}_{\star}$.

To summarize, under the uniform missing and Case (I), (II)(i) in the non-uniform missing, we can achieve the optimal rate up to a $\mathrm{polylog}(n)$ order. For Case (II)(ii), when the missingness is not extreme, with an appropriately chosen $\beta>0$, the proposed estimator can also attain the optimal rate up to the $\mathrm{polylog}(n)$ order.

5.3 Comparison with Uniform Objective Function

Recall that the unweighted empirical risk function $\widehat{R}_{\text{UNI}}(\bm{A})=n^{-2}\|\bm{W}\circ(\bm{A}-\bm{Y})\|_{F}^{2}$ is adopted by many existing matrix completion techniques (Klopp, 2014). An interesting question is whether there is any benefit in adopting the proposed weighted empirical risk function for matrix completion. In this subsection, we aim to shed some light on this aspect by comparing the non-asymptotic error bounds of the corresponding estimators. Due to the additional complication from the estimation error of the observation probability matrix, we only focus on the weighted empirical risk function with true inverse probability weighting in this section. We will demonstrate empirically in Sections 6 and 7 the benefits of the weighted objective function with estimated weights.

Most existing work with unweighted empirical risk function assume the true missingness is uniform (Candès and Plan, 2010; Koltchinskii et al., 2011). One notable exception is Klopp (2014), where unweighted empirical risk function is studied under possibly non-uniform missing structure. The estimator of Klopp (2014) is equivalent to our estimator when $\beta=0$, which is denoted by $\widehat{\bm{A}}^{\text{UNI}}$. Thus, according to Theorem 4, we have with probability at least $1-3/(2n)$,

[TABLE]

which is the same upper bound obtained in Klopp (2014). Define $\widehat{\bm{A}}^{\text{KNOWN}}$ as the estimator which minimizes the known weighted empirical risk function (3). Then,

[TABLE]

The improvement in the upper bounds of the weighted objective function $\widehat{R}$ lies in that, under non-uniform missingness, $\theta_{U}\theta_{L}^{-1}>1$ which implies that $U^{\text{KNOWN}}<U^{\text{UNI}}$ as summarized below.

Theorem 5**.**

Assume Conditions C1-C4 holds, and take $\tau_{\text{KNOWN}}=C_{3}\theta_{L}^{-1/2}n^{-3/2}\log^{1/2}n$ and $\tau_{\text{UNI}}=C_{3}\theta_{U}^{1/2}f^{-1}(\mu_{\star})n^{-3/2}\log^{1/2}n$. The upper bound of $d^{2}(\widehat{\bm{A}}^{\text{UNI}},\bm{A}_{\star})$ is the same as $U^{\text{UNI}}$ and the upper bound of $d^{2}(\widehat{\bm{A}}^{\text{KNOWN}},\bm{A}_{\star})$ is the same as $U^{\text{KNOWN}}$. In addition, $U^{\text{KNOWN}}\leq U^{\text{UNI}}$, and $U^{\text{KNOWN}}<U^{\text{UNI}}$ if $\theta_{U}>\theta_{L}$, i.e., the true missing mechanism is non-uniform.

Our approach draws inspiration from the missing value literature, for instance in Chen et al. (2008), which showed that using the estimated parameters in the inverse probability weighting can actually reduce the variance of the parameter of interest; see Theorem 1 of the paper. Given the results of Chen et al. (2008), we would expect using the estimated parameters $\widehat{\bm{\Theta}}_{\beta}$ in the weighting probability would not be inferior to the version with the true parameter $\widehat{\bm{\Theta}}_{\star}$.

6 Simulation Study

6.1 Missingness

This section reports results from simulation experiments which were designed to evaluate the numerical performance of the proposed methodologies. We first evaluate the estimation performances of the observation probabilities in Section 6.1 and then those of the target matrix in Section 6.2. In the simulation, the true observation probabilities $\bm{\Theta}_{\star}$ and the target matrix $\bm{A}_{\star}$ were randomly generated once and kept fixed for each simulation setting to be described below. To generate $\bm{\Theta}_{\star}$, we first generated $\bm{U}_{\bm{M}_{\star}}\in\mathbb{R}^{n_{1}\times(r_{\bm{M}_{\star}}-1)}$ and $\bm{V}_{\bm{M}_{\star}}\in\mathbb{R}^{(r_{\bm{M}_{\star}}-1)\times n_{2}}$ as random Gaussian matrices with independent entries each following $\mathcal{N}(-0.4,1)$. We then obtained $\bm{M}_{\star}=\bm{U}_{\bm{M}_{\star}}\bm{V}_{\bm{M}_{\star}}^{{\mathrm{\scriptscriptstyle T}}}-\bar{m}_{n_{1},n_{2},r_{\bm{M}_{\star}}}\bm{J}$ where $\bar{m}_{n_{1},n_{2},r_{\bm{M}_{\star}}}$ is a scalar chosen to ensure the average observation rate is $0.2$ in each simulation setting. We finally set $\bm{\Theta}_{\star}=\mathcal{F}(\bm{M}_{\star})$ where the inverse link function $f$ is a logistic function.

In our study, we set $r_{\bm{M}_{\star}}=11$, (or $r_{\bm{Z}_{\star}}=10$) and chose $n_{1}=n_{2}$ with four sizes: $600$, $800$, $1000$ and $1200$, and the number of simulation runs for each settings was $500$.

For the purpose of benchmarking, we compared various estimators of the missingness:

the non-constrained estimator $\widehat{\bm{\Theta}}_{\alpha}$ defined in (1); 2. 2.

the constrained estimator $\widehat{\bm{\Theta}}_{\beta}$ defined in (4); 3. 3.

the directly winsorized estimator $\widehat{\bm{\Theta}}_{\textsf{Win},\beta}=\mathcal{F}\{\widehat{\mu}\bm{J}+\mathcal{T}_{\beta}(\widehat{\bm{Z}})\}$; 4. 4.

the 1-bit estimator $\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$ proposed in Davenport et al. (2014) and its corresponding constrained and winsorized versions $\widehat{\bm{\Theta}}_{\textsf{1-bit},\beta}$ and $\widehat{\bm{\Theta}}_{\textsf{1-bit},\textsf{Win},\beta}$; (note that the 1-bit estimator $\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$ imposes the nuclear-norm regularization on the whole $\bm{M}$ instead of $\bm{Z}$, when compared to $\widehat{\bm{\Theta}}_{\alpha}$) 5. 5.

the rank-1 probability estimator $\widehat{\bm{\Theta}}_{\textsf{NW}}$ used in Negahban and Wainwright (2012) where $g_{i.}=n_{2}^{-1}\sum_{j=1}^{n_{2}}w_{ij}$, $g_{.j}=n_{1}^{-1}\sum_{i=1}^{n_{1}}w_{ij}$ and $\theta_{ij,\textsf{NW}}=g_{i.}g_{.j}$; 6. 6.

the uniform estimator $\widehat{\bm{\Theta}}_{\textsf{UNI}}=N/(n_{1}n_{2})\bm{J}$.

For the non-constrained estimator $\widehat{\bm{\Theta}}_{\alpha}$ and the 1-bit estimator $\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$, the parameter $\alpha$ is set according to the knowledge of the true $\bm{M}_{\star}$. For the constrained estimators $\widehat{\bm{\Theta}}_{\beta}$ and $\widehat{\bm{\Theta}}_{\textsf{Win},\beta}$, the constraint level $\beta$ was chosen so that either $5\%$ or $10\%$ of the elements in $\widehat{\bm{Z}}_{\alpha}$ were winsorized. Similarly for $\widehat{\bm{\Theta}}_{\textsf{1-bit},\beta}$ and $\widehat{\bm{\Theta}}_{\textsf{1-bit},\textsf{Win},\beta}$.

To quantify the estimation performance of linear predictor $\bm{M}_{\star}$ and observation probabilities $\bm{\Theta}_{\star}$, we considered the empirical root mean squared errors $\text{RMSE}(\bm{B},\bm{C})$ with respect to any two matrices $\bm{B}$ and $\bm{C}$ of dimension $n_{1}\times n_{2}$, and the Hellinger distance $d_{H}^{2}(\widehat{\bm{\Theta}},\bm{\Theta}_{\star})$ between $\widehat{\bm{\Theta}}$ and $\bm{\Theta}_{\star}$ defined as follows:

[TABLE]

As the estimators $\mathcal{F}^{-1}(\widehat{\bm{\Theta}}_{\alpha})$ and $\mathcal{F}^{-1}(\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha})$ are both low-rank, we also report their corresponding ranks.

Table 1 summarizes the simulation results for the missingness. The most visible aspect of the results is that the proposed estimators $\widehat{\bm{\Theta}}_{\alpha}$ and $\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$ both have superior performance than the two existing estimators $\widehat{\bm{\Theta}}_{\textsf{NW}}$ and $\widehat{\bm{\Theta}}_{\textsf{UNI}}$ by having smaller root mean square errors with respect to $\widehat{\bm{M}}$, Hellinger distances $d_{H}^{2}(\widehat{\bm{\Theta}},\bm{\Theta}_{\star})$ and more accuracy estimated rank of $\bm{M}_{\star}$. Without the separation of $\mu_{\star}$ from $\bm{M}_{\star}$, $\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$ has larger error and Hellinger distance than the proposed estimators. The performance of $\widehat{\bm{\Theta}}_{\textsf{NW}}$ is roughly between the proposed estimators and the uniform estimator $\widehat{\bm{\Theta}}_{\textsf{UNI}}$. Estimator $\widehat{\bm{\Theta}}_{\textsf{UNI}}$ is a benchmark which captures no variation of the observation probabilities.

6.2 Target matrix

To generate a target matrix $\bm{A}_{\star}$, we first generated $\bm{U}_{\bm{A}_{\star}}\in\mathbb{R}^{n_{1}\times(r_{\bm{A}_{\star}}-1)}$ and $\bm{V}_{\bm{A}_{\star}}\in\mathbb{R}^{(r_{\bm{A}_{\star}}-1)\times n_{2}}$ as random matrices with independent Gaussian entries distributed as $\mathcal{N}(0,\sigma_{\bm{A}_{\star}}^{2})$ and obtained $\bm{A}_{\star}=2.5\bm{J}+\bm{U}_{\bm{A}_{\star}}\bm{V}_{\bm{A}_{\star}}^{{\mathrm{\scriptscriptstyle T}}}$. Here we set the standard deviation of the entries in the matrix product $\bm{U}_{\bm{A}_{\star}}\bm{V}_{\bm{A}_{\star}}^{{\mathrm{\scriptscriptstyle T}}}$ to be $2.5$ to mimic the Yahoo! Webscope data set described in Section 7. To achieve this, $\sigma_{\bm{A}_{\star}}=(2.5^{2}/(r_{\bm{A}_{\star}}-1))^{1/4}$. The contaminated version of $\bm{A}_{\star}$ was then generated as $\bm{Y}=\bm{A}_{\star}+\bm{\epsilon}$, where $\bm{\epsilon}\in\mathbb{R}^{n_{1}\times n_{2}}$ has i.i.d. mean zero Gaussian entries $\epsilon_{ij}\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$. The $\sigma_{\epsilon}^{2}$ is chosen such that $\text{SNR}=(\textsf{E}\|\bm{A}_{\star}\|_{F}^{2}/\textsf{E}\|\bm{\epsilon}\|_{F}^{2})^{1/2}=1$, where $\textsf{E}\|\bm{A}_{\star}\|_{F}^{2}=n_{1}n_{2}(r_{\bm{A}_{\star}}-1+2.5^{2})$ implies $\sigma_{\epsilon}=0.5(r_{\bm{A}_{\star}}-1+2.5^{2})^{1/2}$.

For the estimation of the target matrix, we evaluated ten versions of the proposed estimators Proposed_$\widehat{\bm{\Theta}}_{\beta}$$t$, Proposed$\widehat{\bm{\Theta}}_{\textsf{Win},\beta}$$t$, Proposed$\widehat{\bm{\Theta}}_{\alpha}$, Proposed_$\widehat{\bm{\Theta}}_{\textsf{1-bit},\beta}$$t$, Proposed$\widehat{\bm{\Theta}}_{\textsf{1-bit},\textsf{Win},\beta}$$t$ and Proposed$\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$. Here Proposed indicates the estimators are obtained by solving problem (3), while $\widehat{\bm{\Theta}}_{\beta}$, $\widehat{\bm{\Theta}}_{\textsf{Win},\beta}$, $\widehat{\bm{\Theta}}_{\alpha}$, $\widehat{\bm{\Theta}}_{\textsf{1-bit},\beta}$, $\widehat{\bm{\Theta}}_{\textsf{1-bit},\textsf{Win},\beta}$ and $\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$ represents the probability estimators used in (3), as described in Section 6.1, and $t=0.05$ or $0.1$ denote the winsorized proportion for which $\beta$ is chosen. In addition, same as Mao et al. (2019), we also compared them with three existing matrix completion techniques: the methods proposed in Negahban and Wainwright (2012) (NW), Koltchinskii et al. (2011) (KLT) and Mazumder et al. (2010) (MHT). Among these three methods, NW is the only one that adjusts for non-uniform missingness. All three methods require tuning parameter selection, for which cross-validation is adopted. See Mao et al. (2019) for more details.

To quantify the performance of the matrix completion, in addition to the empirical root mean squared errors with respect to $\widehat{\bm{A}}_{\beta}$ and $\bm{A}_{\star}$, we used one more measure:

[TABLE]

where $\bm{W}^{\star}$ is the matrix of missing indicator with the $(i,j)$th entry being $(1-w_{ij})$. The test error measures the relative estimation error of the unobserved entries to their signal strength. The estimated ranks of $\widehat{\bm{A}}_{\beta}$ are also reported.

Tables 2 summarize the simulation results for different dimensions $n_{1}$=$n_{2}$ ranges from $600$ to $800$ and two different settings of $r_{\bm{A}_{\star}}=11$. The results of $r_{\bm{A}_{\star}}=11$ for different dimensions $n_{1}$=$n_{2}$ ranges from $1000$ to $1200$ are delegated to Table LABEL:apptab:Asimresr11n2 and the results of $r_{\bm{A}_{\star}}=31$ are delegated to Tables LABEL:apptab:Asimresr31n1-LABEL:apptab:Asimresr31n2 of Section LABEL:appsec:sim in the supplementary material. From the tables, we notice that the ten versions of the proposed methods possess superior performance than the three existing methods by having smaller root mean squared errors and Test Errors. Among the first five proposed methods in the tables, Proposed_$\widehat{\bm{\Theta}}_{\beta}$ is better than Proposed_$\widehat{\bm{\Theta}}_{\alpha}$ for most of the time. It is because that the constrained estimator $\widehat{\bm{\Theta}}_{\beta}$ has much smaller ratio $\widehat{\theta}_{U}/\widehat{\theta}_{L}$ than $\widehat{\bm{\Theta}}_{\alpha}$ which improve the stability of prediction and the accuracy. Another observation is that Proposed_$\widehat{\bm{\Theta}}_{\beta}$0.1 performs better than Proposed$\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$ at most times.

7 Real data application

In this section we demonstrate the proposed methodology by analyzing the Yahoo! Webscope dataset (ydata-ymusic-user-artist-ratings-v1_0) available at

http://research.yahoo.com/Academic_Relations. It contains (incomplete) ratings from 15,400 users on 1000 songs. The dataset consists of two subsets, a training set and a test set. The training set records approximately 300,000 ratings given by the aforementioned 15,400 users. Each song has at least 10 ratings. The test set was constructed by surveying 5,400 out of these 15,400 users, each rates exactly 10 songs that are not rated in the training set. The missing rates are $0.9763$ overall, $0.3520$ to $0.9900$ across users, and $0.6372$ to $0.9957$ across songs. The non-uniformity of the missingness is shown in Figure LABEL:appfig:webuser of Section LABEL:appsec:real in the supplementary material. In this experiment, we applied those methods as described in Section 6 to the training set and evaluated the test errors based on the corresponding test set. As there is no prior knowledge about true parameters $\alpha_{1}$ and $\alpha_{2}$, we suggest to choose $\alpha_{1}$ and $\alpha_{2}$ large enough, say $\alpha_{1}=100$ and $\alpha_{2}=100$, to ensure that the range covers all the missing probabilities. It was noted that $\widehat{\bm{\Theta}}_{\alpha}$ is not sensitive to larger $\alpha$.

Table 3 reports the root mean squared prediction errors, where $\text{RMSPE}=\|\bm{W}^{test}\circ(\widehat{\bm{A}}_{\beta}-\bm{Y})\|_{F}/(\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}w_{ij}^{test})^{1/2}$ and $\bm{W}^{test}$ is the indicator matrix of test set with the $(i,j)$th entry being $w_{ij}^{test}$. Note that Proposed_$\widehat{\bm{\Theta}}_{\beta}$0.05 performs the best among all ten versions of proposed methods. Besides, Proposed$\widehat{\bm{\Theta}}_{\alpha}$ also has much smaller root mean squared prediction error than the other eight versions of proposed methods. This may indicate that only slight constraint is required for the probabilities estimator for this dataset. Note that we cannot guarantee the optimal convergence rate or even asymptotic convergence in certain setting of missingness for Proposed_$\widehat{\bm{\Theta}}_{\alpha}$, see Section 5.2 for details.

With the separation of $\mu$, Proposed_$\widehat{\bm{\Theta}}_{\alpha}$ is better than Proposed_$\widehat{\bm{\Theta}}_{\textsf{1-bit},\alpha}$; analogously, Proposed_$\widehat{\bm{\Theta}}_{\beta}$$t$ is better than Proposed$\widehat{\bm{\Theta}}_{\textsf{1-bit},\beta}$$t$ with different constraint level $t$, same to Proposed$\widehat{\bm{\Theta}}_{\textsf{Win},\beta}$$s$ and Proposed$\widehat{\bm{\Theta}}_{\textsf{1-bit},\textsf{Win},\beta}$_$s$ with different winsorization level $s$.

As compared with the existing methods NW, KLT and MHT, our proposed methods perform significantly better in terms of root mean squared prediction errors, and achieve as much as 25% improvement when compared with Mazumder, Hastie and Tibshirani’s method (the best among the three existing methods). This suggests that a more flexible modeling of missing structure improves the prediction power.

8 Concluding Remarks

When the matrix entries are heterogeneously observed due to selection bias, this heterogeneity should be taken into account. This paper focuses on the problem of matrix completion under low-rank missing structure. In the recovery of probabilities of observation, we adopt a generalized linear model with a low-rank linear predictor matrix. To avoid unnecessary bias, we introduce a separation of the mean effect $\mu$. As the extreme values of probabilities may lead to unstable estimation of target matrix, we propose an inverse probability weighting based method with constrained probability estimates and demonstrate the improvements in empirical perspectives. Our theoretical result shows that the estimator of the high dimensional probability matrix can be embedded into the inverse probability weighting framework without compromising the rate of convergence of the target matrix (for an appropriately tuned $\beta>0$), and reveals a possible regime change in the tuning of the constraint parameter ($\beta>0$ vs. $\beta=0$). In addition, corresponding computational algorithms are developed, and a related algorithmic convergence result is established. Empirical studies show the attractive performance of the proposed methods as compared with existing matrix completion methods.