A Pairwise Pseudo-likelihood Approach for Matrix Completion with Informative Missingness

Authors: jiangyuan  li
jiayi  wang
raymond k w wong
kwun  chuen
gary  chan
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 08:46:00

Abstract: While several recent matrix completion methods are developed to deal with non-uniform observation probabilities across matrix entries, very few allow the missingness to depend on the mostly unobserved matrix measurements, which is generally ill-posed. We aim to tackle a subclass of these ill-posed settings, characterized by a flexible separable observation probability assumption that can depend on the matrix measurements. We propose a regularized pairwise pseudo-likelihood approach for matrix completion and prove that the proposed estimator can asymptotically recover the low-rank parameter matrix up to an identifiable equivalence class of a constant shift and scaling, at a near-optimal asymptotic convergence rate of the standard well-posed (non-informative missing) setting, while effectively mitigating the impact of informative missingness. The efficacy of our method is validated via numerical experiments, positioning it as a robust tool for matrix completion to mitigate data bias.

Main paper: 
 introduction: 
 The goal of matrix completion is to recover a target matrix from its noisy and incomplete measurements. It is a modern high-dimensional missing data problem. Despite various significant advances made in the last two decades [e.g., 7, 17, 20, 18] , many works on matrix completion still focus on the missing-at-random mechanism. Although such an assumption is doubtful in many real-life applications, there are very few options available for more general missing data mechanism, especially those with theoretical guarantees. This work aims to provide a principled and theoretically well-supported alternative method for missing-not-at-random settings, where missingness could depend on the measurements that are mostly unobserved.A usual assumption to allow for succeeding matrix completion is that the unknown matrix is lowrank or approximately low-rank. The noiseless setting has been studied in [7] using nuclear norm minimization. The vast majority of existing theories on matrix completion assume that entries are revealed with a constant probability with respect to both entry location and measurement value.Recent approaches to handling entries revealed with nonuniform probabilities, depending on the entry location, have shown the strength to improve matrix completion with solid theoretical guarantees [e.g., 31, 11, 28, 25, 23, 35] . These works aim to mitigate the effects due to the non-uniformity in observation probabilities [30] . When additional row/column attributes are available, it is also possible to use this additional information for handling the non-uniform missing [e.g. , 24] . Although the non-uniform missing mechanism is quite flexible, it is fundamentally different from the missing-notat-random mechanism. The key difference is whether the missing probability depends on the possibly unobserved measurement, which we will highlight below. In a missing-not-at-random setting, the methods developed for the non-uniform missing mechanism could still be biased in matrix recovery. See Section 6 for a numerical example. The method we propose in this work not only deals with (a flexible subclass of) missing-not-at-random settings, but it is also applicable in the non-uniform missing settings mentioned above.In the missing data literature, likelihood-based methods for missing-not-at-random settings commonly involve specifying a parametric distribution of the missing data mechanism. However, this assumption should be used with caution, as it is highly sensitive and may easily induce a misspecified model, resulting in biased estimation and inaccurate results. To circumvent such issues, it is preferable to adopt a missingness assumption as flexible and generally applicable as possible. This type of assumption, often referred to as an unspecified missing data mechanism [37] , avoids explicitly specifying a parametric model. Instead of using the full likelihood for estimation, certain unspecified missing-not-at-random assumption allows for the derivation of a non-standard likelihood [22] , which serves as the foundation for subsequent estimation. Such non-standard likelihood approaches have been used in regression analysis [32] and variable selection when confronted with informative missing [37] . One disadvantage of this approach is that not all the unknown parameters are estimable due to certain non-identification issues [16, 22] .In this work, we extend the pairwise pseudo-likelihood approach [22] to matrix completion with a mild separable informative missingness assumption (see Assumption 2.1 in Section 2), which is very flexible and generally applicable. While not all the parameters are estimable, we can identify the dispersion-scaled matrix up to a constant shift without suffering from bias due to informative missingness. This shows great promise to be applied in practice, for example, in recommendation systems where the rankings of entries are of interest.Apart from the informative missing mechanism, our matrix completion method is based on the exponential family model, which has received extensive attention within the matrix completion literature for its efficacy in modeling non-Gaussian data, particularly discrete data. Notably, researchers have investigated its application in specific scenarios such as one-bit matrix completion [10] and multinomial matrix completion [4, 19] . The application of the exponential family model also extends to accommodating unbounded non-Gaussian observations, including Poisson matrix completion [8] and exponential family matrix completion [13, 21] .Overall, the combination of the separable missing-not-at-random mechanism and the exponential family model allows the proposed method to be applicable in a wide range of settings. We summarize the major contributions of this work as follows.1. We formulate the pairwise pseudo-likelihood approach for matrix completion under informative missingness and exponential family model. To the best of our knowledge, the pairwise pseudo-likelihood approach has never been adopted in the matrix completion setup before. As opposed to the classical applications of pairwise pseudo-likelihood that assume i.i.d. sampling, matrix completion problems exhibit a non-identical and high-dimensional sampling structure. 2. We investigate the identifiability issues of the crucial separable missingness structures (Assumption 2.1) which lies at the core of the pairwise pseudo-likelihood approach. 3. We provide a non-trivial convergence analysis of the proposed estimator up to an identifiable equivalence class. Such analysis involves a novel concentration analysis of matrix-valued U -statistics where existing works on this type of concentration is sparse.To the best of our knowledge, we are one of the first works that consider the missing not at random (MNAR) setting in matrix completion and provide solid theoretical guarantees. [2] claims that they can deal with the MNAR setting. However, they assume selections and noise are independent conditioned on latent factors, as shown in their Assumption 2. On the contrary, our setting allows missingness to depend on noise. [15] also addresses informative missingness in matrix completion. However, they require additional covariate information to complete the matrix. Compared to the above two works, our setting is more general as we do not require independence between selections and noise given the true matrix, and we do not need additional covariate information. However, we do require that ϵ i is independent of ϵ j for i ̸ = j or |i -j| > 1, while [2] allows selection to be dependent among different entries. Regarding to the theoretical bounds, [2] requires additional technical conditions to develop finite sample error bounds, and their bound is point-wise, i.e., the bound is for a given location i. [15] also requires additional conditions on the likelihood and restricted eigenvalues to obtain convergence. Our error bound is developed under relatively weak conditions and achieves the minimax convergence rate.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: The main contributions are adequately described in the abstract and introduction.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
models: 
 Let A ⋆ = (A ⋆,ij ) m1,m2 i,j=1 ∈ R m1×m2 be the matrix of interest, which is related to the observation through a generalized linear model. More specifically, we posit that the measurements Y ij of the (i, j)-th entry possesses a probability density/mass function of the exponential family form:EQUATIONwhere h : R → R + and G : R → R are the base measure function and log-partition function associated with the canonical representation, and ϕ ⋆ > 0 is the dispersion parameter. Note that this family of distributions covers a wide variety of popular distributions including Bernoulli, Gaussian, and Poisson distributions. For matrix completion problems, we do not have measurements from every single entry. Let T ij be the observation indicator variable of the (i, j)-th entry, with value 1 if Y ij is observed and 0 otherwise. We assume that {((Y ij , T ij ) : i = 1, . . . , m 1 ; j = 1, . . . , m 2 } are independent.Uniform-sampling-at-random (USR) mechanism is regarded as one of the simplest missing structures for matrix completion. Under USR, Pr(T ij = 1|Y ij ) is a constant across all i, j, which implies that the observation indicator T ij is independent of the measurement Y ij . While this has been widely used to simplify theoretical analyses in many prior ground works [e.g., 7, 6, 17] , USR is a strong assumption that can be violated in many applications. To address this issue, a few analyses and methods [e.g., 31, 11, 28, 18, 4, 25, 23, 35] have been developed based on the non-uniform missing structures, where Pr(T ij = 1|Y ij ) = t ij for 0 < t ij ≤ 1.Here the observation probabilities are allowed to differ across i, j, but the missingness remains independent of the measurement Y ij . In this paper, we relax this restriction and allow whether an entry is observed or not to depend on the corresponding possibly unobserved measurement, leading to a challenging missing-not-at-random (MNAR) setup.Matrix completion under general MNAR is ill-posed, leading to non-identifiability of A ⋆ (even under standard low-rank assumption). Indeed, general MNAR is ill-posed [33] not only in matrix completion, but also in regression [28, 31] and statistical inference [26] in general. However, some additional structure imposed within the MNAR setting can ensure identifiability. To proceed, we make the following assumption, which corresponds to a flexible subclass of MNAR settings. This assumption makes it possible to identify A ⋆ up to some equivalence relations (see Section 3). Assumption 2.1. The observation probability is separable in the following sense:Pr(T ij = 1|Y ij ) = t ij s(Y ij ), for some t ij ∈ (0, 1] and some non-negative function s(•) : R → R + .As will be made clear later, the proposed technique does not require the knowledge of s(•) and {t ij }.A similar condition has been widely used in various regression problems [e.g., 22, 29, 37] . These works posit an i.i.d. setup with additional covariates, while our setup does not imply an identically distributed assumption across locations and has no covariates. Moreover, in our setup, it is not possible to observe replicates in the same location, while the i.i.d. setup generally allows replicates. Assumption 2.1 is flexible and widely applicable. Not only does it accomodate USR, it also includes non-uniform missing mechanism as a special case, where we can set s(•) ≡ 1 and leave {t ij } variable to account for the non-uniform missing. Obviously, as the observation probability is allowed to depend on possibly unobserved Y ij , it also includes many MNAR settings.Clearly, we only have access to the observed data, i.e., Y ij conditional on T ij = 1. To estimate A ⋆ , we first look at the observed data likelihood of the (i, j)-th entry:Pr(Y ij |T ij = 1; A, ϕ) forA ∈ R m1×m2 and ϕ > 0. By the Bayes' Theorem and Assumption 2.1,Pr(Y ij |T ij = 1; A, ϕ) = Pr(T ij = 1|Y ij )f ij (Y ij ; A, ϕ) Pr(T ij = 1|Y ij )f ij (Y ij ; A, ϕ)dY ij = s(Y ij ) 1 s(y)f ij (y; A, ϕ)dy f ij (Y ij ; A, ϕ) = s(Y ij )b ij (A, ϕ)f ij (Y ij |X; A, ϕ),where b ij (A, ϕ) = 1/ s(y)f ij (y|A, ϕ)dy. We see that the conditional likelihood involves unknown functions s(•) and b ij (•), which makes the estimation of A ⋆ difficult. To address this issue, we adopt a pseudo-likelihood approach [22] based on local ranks.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: The limitation of this work such as identifiability issue is discussed in Section 4.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
pseudo-likelihood approach: 
 Let E = (e 1 , . . . , e n ) ⊆ {1, . . . , m 1 } × {1, . . . , m 2 } be a lexicographically ordered set of n unique locations {(i, j) : T ij = 1}. (Indeed, the specific choice of ordering does not matter.) Let the corresponding measurements be Y = ( Ỹ1 , . . . , Ỹn ) := (Y e1 , . . . , Y en ) and the observation indicator be T = ( T1 , . . . , Tn ) := (T e1 , . . . , T en ). We also write Ãk = Ãe k for k = 1, . . . , n. We decompose the vector Y into two vectors: the order statistics Y (•) = ( Ỹ(1) , . . . , Ỹ(n) ) and the rank statistics R = (R 1 , . . . , R n ). Precisely, Ỹ(j) is the j-th smallest entry in Y and R k is the rank of the k-th entry in Y. To motivate the proposed pseudolikelihood in (3), we first consider the conditional likelihood based on the full rank statistics given the observed data:EQUATIONwhere Ξ is the set of all one-to-one maps from {1, . . . , n} to {1, . . . , n}, i.e., permutations. We notice that (2) does not involve unknown components s(•) and t ij due to the separable assumption (Assumption 2.1), and does not depend on the base measure h(•) and the log-partition function G(•). However, (2) is computationally infeasible due to the summation over all permutations. The proposed pairwise pseudo-likelihood consider local ranks for pairs of observations. For any k and k ′ , let R L kk ′ denote the local rank statistic of Ỹk and Ỹk ′ among the pair ( Ỹk , Ỹk ′ ). We denote Y L (k,k ′ ) as the local order statistics (min{ Ỹk , Ỹk ′ }, max{ Ỹk , Ỹk ′ }). Instead of the full conditional probability (2), we consider the product of all possible combinations of the local rank conditional probability on observations:EQUATIONSimilar to (2) , this pairwise pseudo-likelihood (3) (of A and ϕ) does not contain unknown functions and quantities. However, unlike (2), it does not involve all permutations and is therefore significantly easier to compute.The negative logarithm of the pairwise pseudo-likelihood readsEQUATIONwhereR kk ′ (ϕ -1 A) = exp{-( Ỹk -Ỹk ′ )(ϕ -1 Ãk -ϕ -1 Ãk ′ )} = exp{-( Ỹk -Ỹk ′ )(ϕ -1 Ãk - ϕ -1Ãk ′ )}. We notice two immediate issues with estimating A ⋆ (and ϕ ⋆ ) via minimizing (4).The values of (4) evaluated at any two pairs (A 1 , ϕ 1 ) and (A 2 , ϕ 2 ) are the same when ϕ -11 A 1 = ϕ -1 2 A 2 .Therefore, (4) does not have the ability to distinguish between these pairs.In other words, if ϕ ⋆ > 0 is unknown, (4) would not be able to identify elements in the equivalence class of A ⋆ under equivalence relation: A ∼ c 1 A for any c 1 > 0. Instead, we try to estimate the dispersion-scaled matrix ϕ -1 ⋆ A ⋆ . Therefore, we considerℓ(A) = 1≤k<k ′ ≤n log(1 + R kk ′ (A)).However, this does not solve all the identifiability issues, and, indeed, ℓ cannot identify a shiftequivalence class described below.Shift Equivalence: Let J be a matrix with all entries being one. Consider A + c 2 J for any c 2 ∈ R.Then ℓ(A + c 2 J) = R kk ′ ((A + c 2 J)) = exp{-( Ỹk -Ỹk ′ )( Ãk + c 2 -Ãk ′ -c 2 )} = exp{-( Ỹk -Ỹk ′ )( Ãk -Ãk ′ )} = ℓ(A).Combining the scale and shift equivalence, we can only estimate A ⋆ up to an equivalence relation A ∼ c 1 A + c 2 J for any c 1 > 0 and c 2 ∈ R, which we will refer to as scale-shift equivalence. We remark that the scale-shift equivalence still allows the identification of much useful information from A ⋆ , such as ranking an arbitrary set of entries of A ⋆ . For example, in recommender system applications, one is mostly interested in the ranking within each row/column. Among the elements in the scale-shift equivalence class, we choose to estimate the following representerEQUATIONby imposing the constraint ⟨J, A⟩ = 0 in the optimization. Here ⟨A,B⟩ = i,j A ij B ij for any matrices A, B ∈ R m1×m2Overall, we propose the following penalized pairwise pseudo-likelihood estimatorEQUATIONwhere ∥A∥ ⋆ and ∥A∥ ∞ (= max i,j A ij ) represent the nuclear norm and the entrywise max norm of a matrix A respectively, and a, λ ≥ 0 are tuning parameters. We also use ∥A∥ F to denote the Frobenius norm of a matrix A. Nuclear norm regularization has been commonly used to promote low-rankness in the estimation [25, 24, 14] . Since ℓ is convex, this optimization is convex. The discussion of the optimization algorithm is given in Appendix C. One natural question is whether there would be further hidden identifiability issues beyond scale-shift equivalence. In Section 5, we will provide a finite-sample error bound of the proposed estimator (6) based on the pairwise pseudo-likelihood, which indicates convergence to Ā⋆ , eliminating the possibility of additional identifiability issues.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: The data used in the paper has been cited and linked properly (Section 7). The implementation detail has been described in detail in Appendix C and the code will be made available on GitHub after acceptance.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
identifiability based on separable assumption: 
 One of the major difficulties associated with informative missing is non-identifiability. We first emphasize the non-identifiability for constant shift is not an artifact of the pseudo-likelihood approach. The root cause is the informative missingness (Assumption 2.1). Here is a simple univariate example inspired from [27] to illustrate this point. Suppose we observe from two data-generating models, whose observations are identical in distributions.Model I: Y 1 ∼ N (-1, 1) with observation probability Pr(T 1 = 1|Y 1 = y) = exp(y) 1+exp(y) , thenPr(T 1 = 1, Y 1 = y) = p N (y + 1) exp(y) 1 + exp(y),where p N (•) is the p.d.f. of standard normal distribution.Model II: Y 2 ∼ N (0, 1) with observation probability Pr(T 2 = 1|Y 2 = y) = exp(-1/2) exp(-y) 1+exp(-y) , thenPr(T 2 = 1, Y 2 = y) = p N (y + 1) exp(-1) exp(y + 1) exp(-y) 1 + exp(-y) =p N (y + 1) exp(y) 1 + exp(y) = Pr(T 1 = 1, Y 1 = y).Extending to the matrix form, the observation probabilities of the following two models, where Y 1 , Y 2 ∈ R m1×m2 , are exactly the same. Model I: vec(Y 1 ) ∼ N (-1, I), t 1,ij = 1 for any (i, j) and s 1 (y) = exp(y) 1+exp(y) . Model II: vec(Y 2 ) ∼ N (0, I), t 2,ij = 1 for any (i, j) and s 2 (y) = exp(-1/2) exp (-y) 1+exp(-y) . As such, under Assumption 2.1, we cannot identify the constant shift. We note that, low-rank assumption generally would not provide enough additional information to eliminate this identifiability issue, as constant shift corresponds to at most a rank-1 perturbation.The identification of the dispersion parameter is a difficult task because of the fact that at most one observation is available for each entry. Interestingly, as we have shown in the Appendix (Theorem B.2), under Assumption 2.1, the identification of the dispersion parameter is actually feasible in Gaussian distributions with replicates. However, it is unclear whether the dispersion parameter can be identified in a typical matrix completion setup, which often does not allow replicates. That said, previous works on exponential family matrix completion [21, 13] assume the dispersion parameter is known, under which there would not be a related identifiability issue.Lemma A.5. For any collection of individual index pairs {(j, j ′ ) : 1 ≤ j < j ′ ≤ n}, (a) (From Lemma S.4. in [36] ) When n is even, we can decompose it into (n -1) groups such that within each group, there are n/2 pairs and no repeated individuals.(b) When n is odd, we can decompose it into n groups such that within each group, there are (n -1)/2 pairs and no repeated individuals.Proof. The proof for part (a) is done in [36] .For part (b), when n is odd, we consider an extra index n + 1 and add all the pairs {(j, n + 1) : 1 ≤ j ≤ n} to the original collection. For the new collection of individual index pairs {(j, j ′ ) : 1 ≤ j < j ′ ≤ n + 1}, since n + 1 is even, we can apply part (a) and get n groups such that within each group, there are (n + 1)/2 pairs and no repeated individuals. Therefore, every index appears in each group exactly once. In each group, we remove the pair with n + 1 in it. We now obtain n groups such that within each group, there are (n -1)/2 pairs and no repeated individuals for the collection of individual index pairs {(j, j ′ ) : 1 ≤ j < j ′ ≤ n}. 
theoretical guarantee: 
 Recall that A ⋆ ∈ R m1×m2 . We denote some convenient notation for dimensions, i.e., m = min{m 1 , m 2 }, M = max{m 1 , m 2 }, d = m 1 + m 2 . We use the notation ≲ (≳) to denote less (greater) than up to an absolute multiplicative constant. We write a ≍ b if a ≲ b and b ≳ a. Furthermore, define π L = min i∈[m1],j∈[m2] Pr(T ij = 1) and π U = max i∈[m1],j∈[m2] Pr(T ij = 1). We use [n] to represent {1, . . . , n} for integer n. In this section, we derive the convergence of ∥ A -Ā⋆ ∥ F . Recall that Ā⋆ , defined in (5) , is the representer in the equivalence class of A ⋆ . Assumption 5.1. The following conditions hold.(C1) There exists an absolute constant ρ > 0 such that π U /π L ≤ ρ.(C2) There exists a constant B such that ∥Y∥ ∞ ≤ B almost surely.(C3) There exists some constant κ > 0 (where κ can depend on∥ Ā⋆ ∥ ∞ ) such that E(Z 2 ij,i ′ j ′ ) ≥ κ for any i, i ′ ∈ [m 1 ], j, j ′ ∈ [m 2 ],whereZ ij,i ′ j ′ = (Y ij -Y i ′ j ′ ) × exp((Y ij -Y i ′ j ′ )(A ij -A i ′ j ′ )/2) 1 + exp((Y ij -Y i ′ j ′ )(A ij -A i ′ j ′ )) .Condition (C1) is posited to avoid some specific entries being sampled with very low probability in a relative sense, where the trace-norm penalization fails to work [11, 31] . Note that both π U and π L are allowed to diminish to zero as m 1 , m 2 → ∞, but Condition (C1) implies that their diminishing orders are the same. Condition (C2) is a technical assumption for analyzing the concentration inequalities of the involved U-statistics in pairwise pseudolikelihood. Note that this does not violate the parametric assumption on the distribution of Y . For example, truncated normal distribution satisfies both. We leave the extension to a light-tail type of assumption for future work. Condition (C3) is a technical condition, and is used in deriving the (expected) Hessian of the loss function with respect to A (see (8) in Appendix A). Note that (expected) Hessian of the loss is often important for deriving the convergence rate and so it is reasonable that a related term shows up in our condition. Indeed, this assumption Here, we provide further discussion to show that it is indeed a mild condition.Intuitively, it posits a positive lower bound for an expectation of a squared random variable. This expectation is always non-negative and is zero only when Z ij,i ′ j ′ is exactly zero almost everywhere.For noisy matrix completion settings, this assumption is very mild because, when there are noises, this variable is not exactly zero almost surely. Next, we show that with the exponential family model, we can explicitly characterize κ. First note that when Ā⋆ ∞ ≤ a as assumed in Theorem 5.3 and∥Y∥ ∞ ≤ B as in Condition (C2), we have EZ 2 ij,i ′ j ′ ≥ exp(4aB) [1+exp(4aB)] 2 E Y 2 ij + Y 2 i ′ j ′ -2Y ij Y i ′ j ′ ≥ exp(4aB) [1+exp(4aB)] 2 [Var(Y ij ) + Var(Y i ′ j ′ )]. Recall the density of Y ij (1), one can derive Var(Y ij ) = G ′′ (A ⋆,ij )ϕ, where G ′′ (•) is nonnegative, from the well-known variance formula for exponential family. Therefore one can take κ = min i,jexp(4aB) [1+exp(4aB)] 2 2ϕ[G ′′ (A ⋆,ij )].Lemma 5.2. We have E{∇ℓ(A ⋆ )} = 0, where E(•) is the expectation under the true parameter A ⋆ .Theorem 5.3. Assume rank( Ā⋆ ) ≤ r and Ā⋆ ∞ ≤ a for some positive constants r, a > 0, under Assumptions 2.1 and 5.1, if we further assume mπ U ≳ log(d 2 ) and λ ≍ B 2 log(d) m 1 m 2 π U √ M π U , then with probability at least 1 -6/d, the following holds:EQUATIONOur result implies that the penalized pairwise pseudolikelihood approach can consistently estimate Ā⋆ . Note that the difference between rank( Ā⋆ ) and rank(A ⋆ ) is at most 1. So a low-rank assumption on A ⋆ automatically translates to a low-rank assumption on Ā⋆ . Most existing work present the upper bound concerning the number of observed entries n and treat the matrix completion as a trace regression problem [e.g. 28, 18, 5, 25] . One can take n as m 1 m 2 π L in their bound to compare their results with ours. Similar to the bound established in [18] , our bound has two components and matches with the rates in their upper bound (up to some constants and logarithmic factor). [17] and [28] show a bound that has the same order as the first term (up to some constants and logarithmic factor) with some additional assumptions. [17] adopts the uniform sampling and boundedness of the condition number for ∥ Ā⋆ ∥. [28] assumes that the sampling distribution follows a product distribution and the "spikiness ratio" (see α sp in [28] ) is bounded. Besides the above matrix completion methods that use the nuclear norm regularization, the estimators utilizing the max-norm regularization [e.g. 5, 35] establish the same bound as the second term(up to some constants and logarithmic factor) when they assume the max-norm of Ā⋆ is bounded. While the aforementioned methods address various missing mechanisms, it is important to emphasize that none of them can handle MNAR setting, where the missingness may depend on the observations. However, our method can tackle such informative missingness. It is interesting to see that our error bound resembles the same convergence rate as [18] (minimax optimal rate) up to a logarithmic order, despite that our setup allows MNAR mechanism.In terms of theoretical analysis, the most notable distinction between our estimator with other existing ones lies in the objective function. The pairwise pseudo-likelihood we employ imposes unique theoretical challenges. Firstly, the gradient and Hessian are no longer as straightforward as those in the the commonly used squared loss or negative log-likelihood loss (for exponential family). We carefully derive these two terms, expressing them as pairwise summations (see exact forms in Eq. ( 8) and Eq. ( 9)). Secondly, the elements in these pairwise summations are not mutually independent, posing difficulties in establishing the concentration inequality to bound them. Indeed, we need to develop corresponding theoretical tools for tackling the corresponding matrix concentration of a matrix-valued U-statistics. To address this challenge, we leverage the grouping lemma (Lemma A.5) to decouple these summations into different groups where mutual independence holds within each group. To obtain the efficient grouping, the decoupling is applied to those observed entries. Additionally, while the trace regression model provides a convenient tool for analyzing the sampling distribution, it implicitly assumes "sampling with replacement", i.e., every entry can be observed repeatedly. We adopt the framework of the Bernoulli model for the observation indicator to avoid the issue. However, theoretical analysis become more challenging. A conditional argument (see the conditional event E in ( 10)) is developed to address the discrepancy between these two frameworks. In addition, Lemma A.6 is established to marginalize the conditional event.Finally, we remark that, while pseudo-likelihood approaches have been applied in regression analysis [32] and variable selection [37] to deal with informative missingness, such analyses mainly focus on i.i.d. design and usually make direct restricted eigenvalue condition of the (high-dimensional) Hessian matrix. In our problem, the eigenvalue condition is related to the observation probabilities.As in typical analysis of matrix completion, one is interested in the dependence on these probabilities, as they are allowed to diminish as m 1 , m 2 → ∞. As such, we also analyze the corresponding restricted eigenvalue bound, under the complicated grouping nature and identifiability issue. By adapting the techniques aforementioned, we provide a rigorous convergence result in non-i.i.d. design, which involves analyzing the concentration of a matrix-valued U-statistics (i.e., the Hessian matrix). This analysis distinguishes our work from a mere application of standard pseudo-likelihood theory, and the techniques used in the proof contribute to the field on their own merit. 
numerical experiments: 
 We conduct the following simulation study to demonstrate the efficacy of the proposed method. We generate a 50 × 50 matrix A ⋆ with rank r = 5. The observations Y ij are generated from a Gaussian distribution with mean A ⋆,ij and variance σ 2 independently. In our study, we have settings with different variances σ 2 . The probability of each entry being observed is related to the value of the entry itself:P(T ij = 1|Y ij ) = 1/[1 + exp(3Y ij )].Since the observation probability is smaller for larger Y ij , there exists a distinctive distributional shift between the observed and unobserved entries, as shown in Figure We use the observed entries as training data and equally split the unobserved data as validation and test data. We compare our method with SoftImpute [26] , CZ [5] , MFW [35] and SNN [2] . The validation data is used for hyper-parameter tuning in each method. Since the objective function is convex in the proposed method, we only tune the regularization parameter λ, and fix the number of iterations as T = 100 and step size η = 1.0 in Algorithm 1. We use the output of SoftImpute [26] with the same regularization parameter λ as a warm-up initialization to shorten the training time. For SoftImpute [26] , CZ [5] and MFW [35] , we tune the hyper-parameters involved in the optimization and regularization as suggested. As for SNN [2] , we choose uniform weights and spectral threshold suggested in [12] , and choose the number of neighbors between 1 and 2. Due to the identifiability issue, the validation data is also used to learn a shift and scale parameter (via a simple linear regression) for the proposed method, which is then used in reporting error metrics on test data. Before getting into the error metric, a simple check on the bias of each method is via histograms. We pick one run with σ 2 = 1 and plot the distribution of recovered entries without any transformation for each method, as shown in Figure 3 . It shows that only the proposed method is able to mitigate the observational bias due to the underlying informative missing structure, and exhibits a symmetric distribution, while the distributions generated by other methods are left-skewed, due to the leftskewness of distribution of the observed entries.1We added comparisons of the computational time regarding the setting in Figure 2 (σ 2 = 1). The computational times are listed in Table 1 . While incorporating more complex missing mechanisms, our method and SNN also take the most time. One practical way to speed up the computation of our method is to use a stochastic version of Algorithm 1 (i.e., training in batches). The focus of this paper is more on the robust recovery when encountering informative missing, and less on the computational efficiency with the knowledge that it could be theoretically slower than other SVD-based methods. However, our method is still faster than SNN, where both methods consider more complex missing mechanisms. Given the promising statistical properties of the proposed method, a future direction is to develop scalable algorithms for the proposed estimator or its variants.To further validate the effectiveness of the proposed method, we vary the variances σ 2 in the simulation. This setting is designed to differentiate non-uniform missingness and informative missingness. When the variance is small, the informative missingness is less severe, and non-uniform missingness might be used to approximately describe the missing mechanism. When the variance is large, the observational probability is more affected by the outcome as in a typical informative missingness setting. We choose the variances σ 2 = 0.0, 0.2, 0.4, 0.6, 0.8, 1.0. For each setting, we repeat the simulation 9 times and report the average test root mean squared errors (TRMSE) with standard errors, shown in Figure 2 . We see that as the variance gets larger, there is a larger improvement in the proposed method with the design to account for informative missing over other methods. SNN [2] performs the worst when the variance is large, as it mainly borrows information on observed entries which introduces a substantial bias. It demonstrates the robustness of the proposed method in difficult settings where the missing structure is informative.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: All the training and test details can be found in Section 6&7 and Appendix C.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
real data application: 
 In this section, we use three data examples to illustrate the robust performance of the proposed method. These are the Tobacco Dataset [9] , Coat Shopping Dataset [30] and Yahoo! Webscope Dataset1 . These datasets have been used in prior works for the demonstration of matrix completion methods [e.g., 2, 35] . Due to space limitation, we refer the readers to Appendix C.2 for more detailed discussions of the datasets and our analyses. Following the details of the implementation in Section 6, we report the results in Table 2 . For the Coat Shopping Dataset and Yahoo! Webscope Dataset, the evaluations are based on associated test sets from the original data sources. As for the Tobacco Dataset, following [2] , the missing data are randomly generated 100 times according to cigarette sales. Here is a summary of the results.Tabacco Dataset. As we can see from Table 2 , our method only performs worse than SNN for this MNAR dataset, with significantly smaller TRMSE than the other three methods. Note that in this synthetic missing data, the way to generate missingness is adapted from the SNN paper. When one entry is missed in Tobacco dataset, the entries in the following period are also missed. This does not satisfy the assumption of our work. So it is not surprising to see our method perform sub-optimality. However, the performance of our method still remains strong. Overall, our method performs robustly well across all these three datasets. In our comparison, a few alternatives can perform very well in one example, but badly in another. For example, SNN has an excellent performance in the Tobacco Dataset while performing very poorly in the Coat Shopping Dataset. The robust performance of our method is appealing in practice, as the missing mechanism is often unknown.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: The error bars and statistical significance has been shown and discussed in Section 6 and Appendix C.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
conclusion: 
 In this paper, we tackle the matrix completion problem where missingness could depend on the possibly unobserved measurements, constituting a challenging missing-not-at-random setting. The proposed method is developed under a flexible separable missingness assumption, which allows us to develop a pairwise pseudo-likelihood approach. Corresponding identification is investigated. We also provide a non-trivial convergence analysis, as well as some numerical experiments to illustrate the efficacy of the proposed estimation. Due to the flexibility in both the missing structure (separable missingness) and measurement model (exponential family model), the proposed technique would be useful in a wide range of applications.The grouping nature of the proposed method poses an additional burden in computation, particularly when dealing with a large number of observed entries. For future works, we consider adapting the stochastic grouping idea to reduce the computational cost and exploring its application in large-scale recommender systems.A Proof of Theorem 5.3We provide the proof of our theoretical results and discussion about identifiability and pseudolikelihood approach below. The proof follows the roadmap from [18] to construct the convergence, whereas the details differ because we deal with a matrix-valued U-statistic type estimator. We mainly rely on Lemma S.4. from [36] to decouple the dependence in U-statistic structure. See Lemma A.5 in Section A.4.We start by rewriting the pairwise pseudo-likelihood asℓ(A) = 1≤k<k ′ ≤n {ψ( Ỹk\k ′ Ãk\k ′ ) -Ỹk\k ′ Ãk\k ′ } = 1≤i,i ′ ≤m1 1≤j,j ′ ≤m2 T ij T i ′ j ′ {ψ(Y ij\i ′ j ′ A ij\i ′ j ′ ) -Y ij\i ′ j ′ A ij\i ′ j ′ },where Therefore, the first and second order derivatives of ℓ(A) areỸk\k ′ := Ỹk -Ỹk ′ = Y ij -Y i ′ j ′ =: Y ij\i ′ j ′ ,EQUATIONandEQUATIONwhereE ij\i ′ j ′ = E ij -E i ′ j ′ , E ij ∈ R m1×m2is the canonical basis with value 1 at the (i, j)-th entry and 0 elsewhere, vec(X) is the standard vectorization of matrix X and v ⊗2 = vv ⊤ for v ∈ R m1m2 . 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA]Justification: The paper does not release new assets.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Note that∇ℓ( Ā⋆ ) = - 1≤k<k ′ ≤n R kk ′ ( Ā⋆ ) 1 + R kk ′ ( Ā⋆ ) ( Ỹk -Ỹk ′ )(E e k -E e k ′ ) = - (i,j),(i ′ ,j ′ )∈[m1]×[m2] (i,j)≺(i ′ ,j ′ ) R ij,i ′ j ′ ( Ā⋆ ) 1 + R ij,i ′ j ′ ( Ā⋆ ) ( Ỹij -Ỹi ′ j ′ )(E ij -E i ′ j ′ )T ij T i ′ j ′ = - 1 2 (i,j),(i ′ ,j ′ )∈[m1]×[m2] R ij,i ′ j ′ ( Ā⋆ ) 1 + R ij,i ′ j ′ ( Ā⋆ ) ( Ỹij -Ỹi ′ j ′ )(E ij -E i ′ j ′ )T ij T i ′ j ′ ,where (i, j) ≺ (i ′ , j ′ ) means (i, j) appears before (i ′ , j ′ ) with the same ordering rule within E, e.g., dictionary order.Proof of Lemma 5.2. Note thatE{∇ℓ( Ā⋆ )} = E          - (i,j),(i ′ ,j ′ )∈[m1]×[m2] (i,j)≺(i ′ ,j ′ ) R ij,i ′ j ′ ( Ā⋆ ) 1 + R ij,i ′ j ′ ( Ā⋆ ) ( Ỹij -Ỹi ′ j ′ )(E ij -E i ′ j ′ )T ij T i ′ j ′          = -E          (i,j),(i ′ ,j ′ )∈[m1]×[m2] (i,j)≺(i ′ ,j ′ ) E R ij,i ′ j ′ ( Ā⋆ ) 1 + R ij,i ′ j ′ ( Ā⋆ ) ( Ỹij -Ỹi ′ j ′ )(E ij -E i ′ j ′ ) T ij = T i ′ j ′ = 1          . The proof of E R ij,i ′ j ′ ( Ā⋆) 1+R ij,i ′ j ′ ( Ā⋆) ( Ỹij -Ỹi ′ j ′ )(E ij -E i ′ j ′ ) T ij = T i ′ j ′ = 1 = 0 directly follows Theorem 4.1 from [29].To simplify the notation, we denoteS ij,i ′ j ′ = - R ij,i ′ j ′ ( Ā⋆ ) 1 + R ij,i ′ j ′ ( Ā⋆ ) ( Ỹij -Ỹi ′ j ′ )(E ij -E i ′ j ′ )T ij T i ′ j ′ . Recall that m = min{m 1 , m 2 }, d = m 1 + m 2 , M = max{m 1 , m 2 }, |Y k | ≤ B almost surely. Lemma A.1.With the condition that mπ U ≳ log(d 2 ). We havePr ∇ℓ( Ā⋆ ) ≳ B 2 log(d) m 1 m 2 π U M π U ≤ 3 d . Denote C ij,i ′ j ′ = - R ij,i ′ j ′ ( Ā⋆ ) 1 + R ij,i ′ j ′ ( Ā⋆ ) ( Ỹij -Ỹi ′ j ′ ).First of all, we can verify thatE(C ij,i ′ j ′ T ij T i ′ j ′ ) = E {E[C ij,i ′ j ′ T ij T i ′ j ′ | T ij = 1, T i ′ j ′ = 1]} = 0.Define the eventEQUATION, where X k = E i k ,j k for some index (i k , j k ), k = 1, . . . , n. Without loss of generality, we consider the case when n is even. Therefore, by Lemma A.5, we have the following decomposition.∇ℓ(A) | E = n-1 g=1 (k,k ′ )∈Gg C k,k ′ [X k -X k ′ ],whereC k,k ′ = C i k ,j k ,i k ′ ,j k ′ .Within every group G g , there is no repeated index. Therefore, every element in the group G g is independent of other elements in G g conditioned on E.E(C k,k ′ | E) = 0, ∀k, k ′ = 1, . . . , n.Next, we use Matrix Bernstein's inequality to bound∥ (k,k ′ )∈Gg C k,k ′ [X k -X k ′ ]∥ conditioned on the event E. Take S k,k ′ = C k,k ′ [X k -X k ′ ] and we have E (k,k ′ )∈Gg S k,k ′ S ⊺ k,k ′ | E = E[ (k,k ′ )∈Gg C 2 k,k ′ [X k -X k ′ ][X k -X k ′ ] ⊺ | E] ≤ B 2 (k,k ′ )∈Gg [X k -X k ′ ][X k -X k ′ ] ⊺ ≤ B 2   (k,k ′ )∈Gg X k X ⊺ k + X k ′ X ⊺ k ′ + (k,k ′ )∈Gg X k X ⊺ k ′ + X k ′ X ⊺ k   .Take i k and j k as the corresponding row index and column index for X k such thatX k = E i k ,j k . (k,k ′ )∈Gg X k X ⊺ k + X k ′ X ⊺ k ′ = n k=1 E i k ,i k = 1.The last inequality is due to the fact that the diagonal matrices have the operator norm 1.EQUATIONThen by Theorem 6.1 from [34] , we havePr    (k,k ′ )∈g S k,k ′ ≥ c B 2 (1 + ξ n,g )[2 log(d) + 2 log(m 1 m 2 )] | E    ≤d exp{-[2 log(m 1 m 2 ) + 2 log(d)]} = d 1 (m 1 m 2 ) 2 d 2 = 1 (m 1 m 2 ) 2 d ,for some constant c > 0.(k,k ′ )∈g S k,k ′ . By applying the union bound over all the groups, we obtainPr ∇ℓ( Ā⋆ ) ≥ c n-1 g=1 B 2 (1 + ξ n,g )[2 log(d) + 2 log(m 1 m 2 )] | E ≤ n 1 (m 1 m 2 ) 2 d ≤ 1 (m 1 m 2 )d .Marginalize the event E and we havePr ∇ℓ( Ā⋆ ) ≥ c n-1 g=1 B 2 (1 + ξ n,g )[2 log(d) + 2 log(m 1 m 2 )] ≤ 1 (m 1 m 2 )d . Note that n-1 g=1 B 2 (1 + ξ n,g )[2 log(d) + 2 log(m 1 m 2 )] ≤ n-1 g=1 B 2 [2 log(d) + 2 log(m 1 m 2 )] + n-1 g=1 B 2 ξ n,g [2 log(d) + 2 log(m 1 m 2 )] ≤ (n -1) B 2 [2 log(d) + 2 log(m 1 m 2 )] + √ n -1 B 2 [2 log(d) + 2 log(m 1 m 2 )]( G g=1 ξ n,g ) ≤ n B 2 [2 log(d) + 2 log(m 1 m 2 )] + √ n B 2 [2 log(d) + 2 log(m 1 m 2 )]( G g=1ξ n,g ).Pr    ∇ℓ( Ā⋆ ) ≥ c   n B 2 [2 log(d) + 2 log(m 1 m 2 )] + √ n B 2 [2 log(d) + 2 log(m 1 m 2 )]( G g=1 ξ n,g )      ≤ 1 (m 1 m 2 )d.By Lemma A.6, we adopt the bound for n and n-1 g=1 ξ n,g and further obtain thatPr ∇ℓ( Ā⋆ ) ≥ c B 2 [2 log(d) + 2 log(m 1 m 2 )] 3m 1 m 2 π U + √ 3m 1 m 2 π U 2m 1 m 2 M π 2 U ≤ 1 (m 1 m 2 )d + exp(-m 1 m 2 π L ) + 2M exp(-m 1 m 2 π 2 U /2). Pr ∇ℓ( Ā⋆ ) ≳ B 2 [log(d) + log(m 1 m 2 )] m 1 m 2 π U (1 + M π U ) ≤ 1 (m 1 m 2 )d + exp(-m 1 m 2 π L ) + 2M exp(-m 1 m 2 π 2 U /2).With the condition that mπ U ≳ log(d 2 ). We havePr ∇ℓ( Ā⋆ ) ≳ B 2 [log(d) + log(m 1 m 2 )] m 1 m 2 π U M π U ≤ 3 d .The conclusion follows by assimilating universal constants in ≳.Recall that ψ(t) = log ( The Hessian matrix reads∇ 2 ℓ(A) = 1≤k<k ′ ≤n {ψ ′′ (Y ij\i ′ j ′ A ij\i ′ j ′ )Y 2 ij\i ′ j ′ vec(E ij\i ′ j ′ ) ⊗2 } = 1≤k<k ′ ≤n (Y k -Y k ′ ) 2 exp(( Ỹk -Ỹk ′ )( Ãk -Ãk ′ )) (1 + exp(( Ỹk -Ỹk ′ )( Ãk -Ãk ′ ))) 2 vec(E k -E k ′ ) ⊗2 = 1 2 (i,j),(i ′ ,j ′ )∈[m1]×[m2] T ij T i ′ j ′ (Y ij -Y i ′ j ′ ) 2 exp((Y ij -Y i ′ j ′ )(A ij -A i ′ j ′ )) (1 + exp((Y ij -Y i ′ j ′ )(A ij -A i ′ j ′ ))) 2 vec(E ij -E i ′ j ′ ) ⊗2 = 1 2 (i,j),(i ′ ,j ′ )∈[m1]×[m2] T ij T i ′ j ′ Z 2 ij,i ′ j ′ vec(E ij -E i ′ j ′ ) ⊗2 ,whereZ ij,i ′ j ′ = (Y ij -Y i ′ j ′ ) exp((Y ij -Y i ′ j ′ )(A ij -A i ′ j ′ )/2) 1 + exp((Y ij -Y i ′ j ′ )(A ij -A i ′ j ′ )) .Lemma A.2. For any U ∈ R m1×m2 such that ⟨J, U⟩ = 0, and assume thatE(Z 2 ij,i ′ j ′ ) ≥ κ, ∀(i, j), (i ′ , j ′ ) ∈ [m 1 ] × [m 2 ]when ∥A∥ ∞ ≤ a for some constant a > 0, we have thatEvec(U) ⊤ ∇ 2 ℓ(A)vec(U) ≥ κm 1 m 2 π 2 L ∥U∥ 2 F .Proof. We derive the quadratic form first and take the expectation to obtain thatEu ⊤ ∇ 2 ℓ(A)u = 1 2 (i,j),(i ′ ,j ′ )∈[m1]×[m2] E{T ij T i ′ j ′ Z 2 ij,i ′ j ′ (U ij -U i ′ j ′ ) 2 } ≥ 1 2 κ (i,j),(i ′ ,j ′ )∈[m1]×[m2] E{T ij T i ′ j ′ (U ij -U i ′ j ′ ) 2 } ≥ 1 2 κπ 2 L (i,j),(i ′ ,j ′ )∈[m1]×[m2] (U ij -U i ′ j ′ ) 2 = 1 2 κπ 2 L (i,j),(i ′ ,j ′ )∈[m1]×[m2] U 2 ij + U 2 i ′ j ′ -2U ij U i ′ j ′ = κm 1 m 2 π 2 L ∥U∥ 2 F .We now start to lower bound u ⊤ ∇ 2 ℓ(A)u. We consider the following constraint setC = U ∈ R m1×m2 ⟨J, U⟩ = 0 ,Still, we derive the argument by conditioning on the event E defined in (10) .vec(U) ⊤ ∇ 2 ℓ(A)vec(U) | E = 1 2 k,k ′ =1,...,n Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 .where (i k , j k ) is the entry where E i k ,j k = 1.By Lemma A.5, WLOG, we assume n is even. We decompose the summation intovec(U) ⊤ ∇ 2 ℓ(A)vec(U) | E = n-1 g=1 (k,k ′ )∈Gg Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 ,where within every group G g , there is no repeated index. DenoteΣ g = (k,k ′ )∈Gg ε k,k ′ Z i k j k ,i k ′ j k ′ (X k -X k ′ ).where ε k,k ′ are independent Radamacher variables. For convenience, we denote u = vec(U) ∈ R m1m2 for any U ∈ R m1×m2 . Lemma A.3. For all U ∈ C, the following holdsu ⊤ ∇ 2 ℓ(A)u -Eu ⊤ ∇ 2 ℓ(A)u ≳ B 2 ∥U ∥ * log(d) m 1 m 2 π U M π U +B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2with probability at least 1 -3/d.Proof. Denote C T := {U ∈ C : ∥U ∥ * ≤ T }. W T = sup U∈C(T ) u ⊤ ∇ 2 ℓ(A)u -Eu ⊤ ∇ 2 ℓ(A)u , Let v = mπ U /π L . By Lemma A.4, we have with probability at least 1 -1/d, u ⊤ ∇ 2 ℓ(A)u -Eu ⊤ ∇ 2 ℓ(A)u ≳ B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 .Next, we focus on the case when ∥U ∥ * ≥ v. We will show that the probability of the following bad event is smallB = ∃U ∈ C such that u ⊤ ∇ 2 ℓ(A)u -Eu ⊤ ∇ 2 ℓ(A)u ≳ B 2 ∥U ∥ * log(d) m 1 m 2 π U M π U , ∥U ∥ * ≥ v .We use a standard peeling argument. For l ∈ N setS l = U ∈ C(r) : α l ν ≤ ∥U ∥ * ≤ α l+1 ν .For each T > ν, define the following eventB l = ∃U ∈ S l such that u ⊤ ∇ 2 ℓ(A)u -Eu ⊤ ∇ 2 ℓ(A)u ≳ B 2 ∥U ∥ * log(d) m 1 m 2 π U M π U .By lemma A.4, with probability smaller than (m 1 m 2 ) exp -T 2 π L mπ U log(d 3 ) , we have thatW T ≥ B 2 T log(d) m 1 m 2 π U M π U .Therefore, we obtain thatPr(B l ) ≤ (m 1 m 2 ) exp - (α 2l v 2 )π L mπ U log(d 3 ) .Using the union bound, we havePr(B) ≤ ∞ l=0 Pr(B l ) ≤ ∞ l=1 (m 1 m 2 ) exp -(α 2l v 2 )π L mπ U log(d 3 ) + 1/d ≤ ∞ l=1 (m 1 m 2 ) exp -(v 2 )π L mπ U log(d 3 )2l log(α) + 1/d ≤ ∞ l=1 (m 1 m 2 ) exp -(v 2 )π L mπ U log(d 3 )2l log(α) + 1/d ≤ (m 1 m 2 ) exp{-2 log(α) log(d 3 )} + 1/d ≤ 2/d.Combine two cases, and the conclusion follows.We now make up the concentration property on W T . Lemma A.4. Conditioned on event E, For any T > ν, we denoteC(T ) := {U ∈ C(r) : ∥U ∥ * ≤ T } and W T := sup U∈C(T ) u ⊤ ∇ 2 ℓ(A)u -Eu ⊤ ∇ 2 ℓ(A)u .When T ≤ mπ U /π L , we haveP W T ≳ B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 ≤ 1/d. When T ≥ mπ U /π L , By taking t = T √ 3 log d √ M π U m1m2π U , we have P W T ≳ B 2 T log(d) m 1 m 2 π U M π U ≤ (m 1 m 2 ) exp - -T 2 π L mπ U log(d 3 ) .Proof. Still, without generality, we focus on the case when n is even. Recall that conditioned on event E, we haveW T = sup U∈C(r,T ) u ⊤ ∇ 2 ℓ(A)u -Eu ⊤ ∇ 2 ℓ(A)u = sup U∈C(r,T ) g (k,k ′ )∈Gg Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 -E g (k,k ′ )∈Gg Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 ≤ g sup U∈C(r,T ) (k,k ′ )∈Gg Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 -E (k,k ′ )∈Gg Z 2 k,k ′ (U i k j k -U i k ′ j k ′ ) 2 ,andZ g,T = sup U∈C(r,T ) (k,k ′ )∈Gg Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 -E (k,k ′ )∈Gg Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 .Note that conditioned on the event E, within each grouping G g , every term 2 is independent of the others. The standard symmetrization trick still applies here,Z 2 i k j k ,i k ′ j k ′ (U i k j k - U i k ′ j k ′ )EZ g,T | E ≤ 2E sup U∈C(T ) (k,k ′ )∈Gg ε k,k ′ Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 . Since |Z i k i k ′ | ≤ 2B and ∥U∥ ∞ = 1, Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 ≤ 16B 2 for every (k, k ′ ). Therefore, ϕ(u) = u 2 , |ϕ(u) -ϕ(v)| ≤ |u + v||u -v| ≤ 8B|u -v|. The contraction inequality yields EZ g,T | E ≤ 2E sup U∈C(T ) (k,k ′ )∈Gg ε k,k ′ Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 ≤ 16BE sup U∈C(T ) (k,k ′ )∈Gg ε k,k ′ Z i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) ≤ 16BE   sup U∈C(T ) (k,k ′ )∈Gg ε k,k ′ Z i k j k ,i k ′ j k ′ ⟨X k -X k ′ , U⟩ | E   ≤ 16BE   sup U∈C(T ) (k,k ′ )∈Gg ε k,k ′ Z i k j k ,i k ′ j k ′ (X k -X k ′ ), U | E   ≤ 16BE   (k,k ′ )∈Gg ε k,k ′ Z i k j k ,i k ′ j k ′ (X k -X k ′ ) | E   sup U∈C(T ) ∥U∥ ⋆ ≤ 16BT E (∥Σ g ∥ | E) . Take γ k,k ′ (U ) = Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 -EZ 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 Note that |Z 2 i k j k ,i k ′ j k ′ (U i k j k -U i k ′ j k ′ ) 2 | ≤ 16B 2 . Then Eγ k,k (U ) = 0 sup (k,k ′ ) sup U ∈C(T ) γ k,k ′ (U )/(32B 2 ) ≤ 1By Massart's concentration inequality (e.g., Theorem 14.2 in [3] ). We have that conditioning on the event E, for any t > 0,Pr   sup U ∈C(T ) 1 n/2 (k,k ′ )∈Gg γ k,k ′ (U )/(32B 2 ) > E sup U ∈C(T ) 1 n/2 (k,k ′ )∈Gg γ k,k ′ (U )/(32B 2 ) + t   ≤ exp(-(n/2)t 2 /8) Pr   sup U ∈C(T ) 1 n/2 (k,k ′ )∈Gg γ k,k ′ (U ) > E sup U ∈C(T ) 1 n/2 (k,k ′ )∈Gg γ k,k ′ (U ) + 32B 2 t   ≤ exp(-(n/2)t 2 /8) Pr Z g,T > EZ g,T + 16nB 2 t ≤ exp(-nt 2 /16).Pr W T > n-1 g=1 EZ g,T + 16n(n -1)B 2 t | E ≤ (n -1) exp(-nt 2 /16).Next, we marginalize the event E.Pr W T > n-1 g=1 EZ g,T + 16n 2 B 2 t ≤ E[n exp(-nt 2 /16)] ≤ (m 1 m 2 )E[exp(-nt 2 /16)] ≤ (m 1 m 2 ) exp(-E(n)t 2 /16) ≤ (m 1 m 2 ) exp(-m 1 m 2 π L t 2 /16).n-1 g=1 EZ g,T . Note that n-1 g=1 EZ g,T ≤ 16BT E E n-1 g=1 ∥Σ g ∥ | E.By using a similar argument in Lemma A.1, we are able to show that for any x > 0,Pr n-1 g=1 ∥Σ g ∥ ≳ x B 2 log(d) m 1 m 2 π U M π U ≤ 3 d exp{-x log(d 2 )}. Therefore, n-1 g=1 EZ g,T ≲ 16B 2 T log(d) m 1 m 2 π U M π U .Combining the result from Lemma A.6, we havePr W T ≳ B 2 T log(d) m 1 m 2 π U M π U + (m 1 m 2 π U ) 2 B 2 t ≤ (m 1 m 2 ) exp(-m 1 m 2 π L t 2 /16).When T ≤ mπ U /π L , we havePr W T ≳ B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 ≤ 1/d. When T ≥ mπ U /π L , By taking t = T √ 3 log d √ M π U m1m2π U , we have Pr W T ≳ B 2 T log(d) m 1 m 2 π U M π U ≤ (m 1 m 2 ) exp -T 2 π L mπ U log(d 3 ) .A.3 Proof of Theorem 5.3Proof. It follows from the definition of the estimator A that ℓ( A) + λ A From Lemma A.3 and Lemma A.2, we haveu ⊤ ∇ 2 ℓ(A)u ≳ E u ⊤ ∇ 2 ℓ(A)u -B 2 ∥U ∥ * log(d) m 1 m 2 π U M π U -B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 ≳ κm 1 m 2 π 2 L ∥U ∥ 2 F -B 2 ∥U ∥ * log(d) m 1 m 2 π U M π U -B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 .with probability at most 1 -3/d.Therefore with probability at most 1 -3/d,κm 1 m 2 π 2 L ∥∆∥ 2 F ≲ B 2 ∥∆∥ * log(d) m 1 m 2 π U M π U + B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 +vec( A -Ā⋆ ) ⊤ ∇ 2 ℓ( A)vec( A -Ā⋆ ) ≲ B 2 ∥∆∥ * log(d) m 1 m 2 π U M π U + B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 +λ Ā⋆ ⋆ -A ⋆ + ∇ℓ( Ā⋆ ) ∥∆∥ ⋆ .Let {u k ∈ R m1 } and {v k ∈ R m2 } be the left and right singular vectors of Ā⋆ respectively. For any matrix A ∈ R m1×m2 , we let row(A) ⊆ R m2 and col(A) ⊆ R m1 denote its row space and column space respectively. Let the column and row span of Ā⋆ be U ⋆ = col( Ā⋆ ) = span{u k } and V ⋆ = row( Ā⋆ ) = span{v k } respectively. DefineM := {A : row(A) ⊆ V ⋆ , col(A) ⊆ U ⋆ }, M ⊥ := {A : row(A) ⊆ V ⊥ ⋆ , col(A) ⊆ U ⊥ ⋆ }. It is easy to see that M ⊆ M, but M ̸ = M. The subspace compatibility of M is upper bounded by √ 2r, i.e., sup A∈M\{0} ∥A∥ ⋆ ∥A∥ F ≤ √ 2r.We observe thatA ⋆ = Ā⋆ + ∆ M + ∆ M ⊥ ⋆ ≥ Ā⋆ + ∆ M ⊥ ⋆ -∥∆ M ∥ ⋆ = Ā⋆ ⋆ + ∆ M ⊥ ⋆ -∥∆ M ∥ ⋆ . By choosing λ ≳ 2( ∇ℓ( Ā⋆ ) + B 2 log(d) m 1 m 2 π U √ M π U ), we have κm 1 m 2 π 2 L ∥∆∥ 2 F ≲ ∇ℓ( Ā⋆ ) + B 2 log(d) m 1 m 2 π U M π U ( ∆ M ⊥ ⋆ + ∥∆ M ∥ ⋆ ) + λ(∥∆ M ∥ ⋆ -∆ M ⊥ ⋆ ) + B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 ≲ 3λ ∥∆ M ∥ ⋆ + B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 ≤ 3λ √ 2r∥∆∥ F + B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2Then, we could derive∥∆∥ 2 F ≲ 3λ √ 2r κm 1 m 2 π 2 L ∥∆∥ F + B 2 log(d)(m 1 m 2 π U ) 2 (m 1 m 2 π L ) -1/2 κm 1 m 2 π 2 L 1 m 1 m 2 ∥∆∥ F ≲ max λ 2 r m 1 m 2 (κm 1 m 2 π 2 L ) 2 , B 2 log(d) κ π U π L 2 1 m 1 m 2 π L Due to Lemma A.1, Pr ∥∇ℓ(A)∥ ≳ B 2 [log(d) + log(m 1 m 2 )] m 1 m 2 π U M π U ≤ 3 d . we can take λ ≍ (B 2 log(d) m 1 m 2 π U √ M π U ),and with probability at least 1 -6/d, we have1 m 1 m 2 ∥∆∥ 2 F ≲ max B 4 [log d] 2 κ 2 π U π L 3 M r m 1 m 2 π L , B 2 log(d) κ π U π L 2 1 m 1 m 2 π L A.Z ij,i ′ j ′ = (Y ij -Y i ′ j ′ ) 2 exp((Y ij -Y i ′ j ′ )(A ij -A i ′ j ′ )/2) 1 + exp((Y ij -Y i ′ j ′ )(A ij -A i ′ j ′ )) , Σ g = (ij,i ′ j ′ )∈g ε ij,i ′ j ′ T ij T i ′ j ′ Z ij,i ′ j ′ (E ij -E i ′ j ′ ),for any collection of non-overlap index pairs.We analyze the upper bound for E ∥Σ g ∥ to prove Corollary 5.3. Lemma A.6. Take n = i,j T i,j and ξ n,g defined in (11) . We havePr   n ≥ e( i,j π i,j )   ≤ exp(-m 1 m 2 π L ). Pr n-1 g=1 ξ n,g > 2m 1 m 2 M π 2 U ≤ 2M exp -m 1 m 2 π 2 U /2 .Proof. Note that T i,j are independent Bernoulli random variables and En = i,j π i,j . We apply Chernoff's inequality and obtainPr (n ≥ t) ≤ exp{- i,j π i,j } e( i,j π i,j ) t t ,for any t > En. Take t = e( i,j π i,j ) and we havePr   n ≥ e( i,j π i,j )   ≤ exp{- i,j π i,j } ≤ exp(-m 1 m 2 π L ). Note that n-1 g=1 ξ n,g = max    n k̸ =k ′ 1{i k = i k ′ }, n k̸ =k ′ 1{j k = j k ′ }    = max    m1 i=1 m2 j̸ =j ′ T ij T ij ′ , m2 j=1 m1 i̸ =i ′ T ij T ij ′    .We consider bounding n i=1 m2 j̸ =j ′ T i,j T i,j ′ . Similarly, without loss of generality, we assume j is even, and by Lemma A.5, we can decomposem1 i=1 m2 j̸ =j ′ T i,j T i,j ′ = m2-1 g=1 ( m1 i=1 (j,j ′ )∈G ′ g T ij T ij ′ ).Within every group G ′ g , every pair T i,j T i,j ′ is independent of others. Then we apply Bernstein inequality to boundEQUATIONTake a union bound over g = 1, . . . , m 2 -1, we havePr   m1 i=1 m2 j̸ =j ′ T ij T ij ′ > 2m 1 m 2 (m 2 -1)π 2 U   ≤ m 2 exp -m 1 m 2 π 2 U /3 . Pr   m1 i=1 m2 j̸ =j ′ T ij T ij ′ > 2m 1 m 2 M π 2 U   ≤ M exp -m 1 m 2 π 2 U /3 .With the same argument, we havePr   m2 j=1 m1 i̸ =i ′ T ij T ij ′ > 2m 1 m 2 M π 2 U   ≤ M exp -m 1 m 2 π 2 U /3 .Prn-1 g=1 ξ n,g > 2m 1 m 2 M π 2 U ≤ 2M exp -m 1 m 2 π 2 U /3 .Assume Y ij ∼ N (µ ij , σ 2 ) with missing mechanism:EQUATIONwherec ij , s(•) ∈ [0, 1]. Lemma B.1. Assume that lim y→∞ c ij s(y) c ′ ij s ′ (y) = b ij or lim y→-∞ c ij s(y) c ′ ij s ′ (y) = b ′ ij ,where b ij , b ′ ij can not be both 0 or ∞ simultaneously.Initialize: Initialize A(0) randomly, set learning rate η.for t = 0 to T do K(t) = A(t) -η∇ℓ(A(t)) Q(t) = S λ (K(t)) A(t + 1) = POCS(Q(t)) end for Algorithm 2 POCS Initialize: Input matrix Q ∈ R m1×m2 . t = 0. while Q ′ ̸ = Q do Q = Q ′ Q = Q - 1 m1m2 m1 i=1 m2 j=1 Q i,j J Q ′ i,j = 1 | Qi,j | ≤ α Qi,j + 1 | Qi,j | > α sign Qi,j α end whileDue to space constraints, we present the test mean absolute errors (TMAE) curve (Figure 4 ) for simulation setting in Section 6. Note that the trend is consistent on both TRMSE and TMAE with multiple runs.We conduct more simulation studies in higher dimensions with a different missing mechanism where the observation probability readsPr(T ij = 1|Y ij ) = 1 1 + exp(-3(Y ij -2)), which means larger potential observations imply higher observation probabilities, as shown in Figure 5 . We generate a 100 × 100 matrix A ⋆ with rank r = 5. The observations Y ij are generated from a Gaussian distribution with mean A ⋆,ij and variance 1 independently. Note that compared with experiments shown in Section 6, the dimension is larger and the observation probability is flipped.As shown in Figure 6 , other methods suffer from a severe observation bias, and the recovered entries are right-skewed, whereas the distribution of true entries is symmetric. Our method alleviates the observation bias and recovers the symmetric pattern of the distribution on recovered entries. The test root mean squared errors (TRMSE) and test mean absolute errors (TMAE) of recovered entries are reported in Table 3 . The experiments are repeated with 9 runs, and the standard deviation of both metrics is included. Our method shows a significant advantage when the observation bias persists. In this section, we use three data examples to illustrate the robust performance of the proposed method. Similar to the implementation procedures in Section 6, we equally separate half of the testing data points and use them as the validation dataset to tune the hyper parameters for all the methods mentioned in Section 6 and learn the shift and scale parameters for the proposed method. Then we used the remaining half of the test data points to evaluate their performance.This dataset is available in Table 11 in [9] and has been widely studied for synthetic control methods [e.g. 1]. From Table 11 in [9] , we can obtain the Tax-Paid Per Capita Sales in Number of Packs for 51 states across the US from year 1950 to year 2014. California implemented a large-scale tobacco control program in 1988 and we consider it as the "treated" state. We take the remaining 50 states as the control states. For the detailed background, we refer readers to [1] and [2] .We collect data from year 1970 to year 2000 and restrict our focus to the 50 control states, which results in a 50 by 31 matrix. Following the same experimental setup in [2] , we generate MNAR data. We introduce "interventions" to a subsets of states in 1989 based on their change in the mean of cigarette sales during 1989-2000 versus that during 1970-1988. See details of the intervention probabilities in Section 6.3 in [2] . As long as an intervention is adopted for state i, all sales under control after 1988 are unobserved, i.e. T i,j = 0 for j > 19. And we take Y i,j , j > 19 as the test data points.Table 2 provides the results for five methods under 100 randomizations on the intervention based on the intervention probability for every state. As we can see, our method only performs worse than SNN for this MNAR dataset, with significantly smaller TRMSE than other three methods. Note that in order to compare the errors, we need to perform a transformation on our estimated matrix. For this study, the untransformed data can also provide valuable information about the trend of sales change for every state during this 30 years. For example, as illustrated in Figure 7 , using untransformed estimated result from our method, we are able to capture the overall trend of the sales change for state KS across 30 years. Our method is able to capture the increasing trend of sales for state KS after year 1988. This dataset is available at https://www.cs.cornell.edu/~schnabts/mnar/. It contains ratings from 290 Turkers on an inventory of 300 items [30] . Training set contains 6960 self-selected ratings and the test set consists of 4640 entries. This dataset has been used as an illustration for the nonuniform missingness [e.g. 30, 35] .As Table 2 shows, SNN performs a lot worse than the remaining methods for this dataset. MFW has the smallest TRMSE. Our method has smaller errors than SoftImpute and has comparable performance to CZ.This dataset is available at https://webscope.sandbox.yahoo.com/catalog.php?datatype= r&did=3, which contains ratings from 15,400 users on 1000 songs. The IDs for users and songs are randomly assigned. The training set includes approximately 300,000 ratings from these 15,400 users. The ratings from the training set are collected during normal use of Yahoo! Music services for each user. The test set was constructed by surveying the first 5,400 users. And each surveyed-user provides ratings for exactly 10 additional songs.Due to its large size and to simplify the computation, we conducted a selection procedure to reduce the size of the matrix. First, we focus on the users with ID 1-50, 5401-5550 and songs with ID 1-250, which results in a matrix where we have 50 surveyed users and 150 unserveyed users. Next, we construct the training (test) matrix from the original training (test) dataset with selected user IDs and selected song IDs. Then, we remove those users and songs who have no single observation in the training matrix. In the end, we obtain a matrix with 199 users and 219 songs.From Table 2 , we can see that the two methods (SNN and our method) that are designed for MNAR have better performance than the remaining methods, and our method has the smallest TRMSE.Justification: There is no societal impact of the work performed.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).https://webscope.sandbox.yahoo.com/catalog.php?datatype=r&did=3⋆ ≤ ℓ( Ā⋆ ) + λ Ā⋆ ⋆ , equivalently ℓ( A)ℓ( Ā⋆ ) ≤ λ Ā⋆ ⋆ -A ⋆ which implies ⟨∇ℓ( Ā⋆ ), vec( A -Ā⋆ )⟩ + vec( A -Ā⋆ ) ⊤ ∇ 2 ℓ( A)vec( A -Ā⋆ ) ≤ λ Ā⋆ ⋆ -A ⋆ .whereA = t Ā⋆ + (1t) A for some t ∈ [0, 1].Let's denote ∆ = A -Ā⋆ . 


Decision: Good paper.

Reviews: review: paper_summary: paper_summary: Current matrix completion approaches mostly deal with non-uniform sampling scenarios. However, few  approaches can allow the missingness to depend on the mostly unobserved matrix measurements. In this paper, the authors propose a regularized pairwise pseudo-likelihood approach for matrix completion in this scenario. The proposed estimator can recover the low-rank parameter matrix asymptotically up to an identifiable equivalence class of a constant shift and scaling, at a near-optimal asymptotic convergence rate of the standard well-posed setting, while effectively mitigating the impact of informative missingness.
strength_weakness: strength_weakness: However, there is a lack of literature review related to the informative missing setting and a lack of comparison with related methods aimed at the informative missing setting, making it difficult to accurately evaluate the contribution. 

Also, I am unsure if these scenarios are widely applicable. In the experimental section, the proposed method does not achieve the best results on some datasets. For instance, on the Tabacco Dataset, the method does not outperform the SNN method. 

Lastly, in line 350, the author says the proposed method's grouping nature adds an extra computational burden, thus this approach seems not to be computationally efficient.
questions: questions: 1)	I still am unclear about the meaning of “informative missingness”, could the author give a more detailed definition and a concrete example? 

2)  Are the authors the first to consider this informative missingness setting in matrix completion? And are there any other methods to address informative missingness aimed at matrix completion? And How to evaluate the strengths and weaknesses of these methods ( including the method proposed in this paper).

3)	What is the difference between the pseudo-likelihood method and Bayes modeling?

4)	Please provide a comparison of the computational time between the proposed method and the previous algorithms for matrix completion tasks. This is crucial as people often deal with high-dimensional and massive datasets, where computational efficiency is a key factor. If possible, please also provide the errors (TRMSE) versus the iterations/time during the trajectories.


5)	It might be challenging to effectively identify whether a scenario qualifies as "informative missingness" and apply this approach. For instance, on the Tabacco Dataset, the method does not outperform the SNN method. Is this result arise from this phenomenon/dataset that doesn’t belong to the “informative missingness” scenario? How to identify a phenomenon/dataset belonging to the “informative missingness” scenario?
limitations: limitations: 1) I'm uncertain about the article's claim that "informative missingness" is widespread and significant in real-world situations. It might be challenging to effectively identify whether a scenario qualifies as "informative missingness" and apply this approach. For instance, on the Tabacco Dataset, the method does not outperform the SNN method. This dataset may not fall into the “informative missingness” scenario. How to identify whether a phenomenon is "informative missingness"?

2) There is not an extensive literature review on informative missingness in matrix completion or other related problems, please provide a detailed literature review about informative missingness and related approach comparisons. 

3)  In line 350, the author says the proposed method's grouping nature adds an extra computational burden and I am not sure whether this approach is computationally efficient in practice.
paper_summary: paper_summary: The paper proposes a novel matrix completion method where the entries are missing not at random. The author(s) tackle the problem using a penalized pairwise pseudo-likelihood approach. The key contributions of the paper are (a) a flexible separable observation probability assumption that depend on the measurements and (b) to provide a novel asymptotic result concerning the concentration bound on the proposed estimator. The success of the proposed method is validated via numerical experiments and on three real data analysis.
strength_weakness: strength_weakness: Following are some of the weaknesses of the paper

(a)  A bit more discussion on how the constants $\kappa$ and $\rho$ influence the right hand bound on equation (7) would be insightful.

(b) Some comparison of the computational times of the competing approaches vs the proposed approach would be helpful.

(c) In the simulations, apart from Gaussian distribution mechanism, it would have been interesting to see situations where observations have been generated from a skewed distribution.
questions: questions: I have the following queries for the authors-

(a) In assumption 2.1 is it required to have $s(Y_{ij})$ distinct at least for one pair $(i,j)$ and $(i^{'},j{'})$ ? Otherwise we may end up with a situation where the assumption 2.1 is another variant of non-uniform missing structure scenario.

(b) Why do we need the constraint $\ || A || \_{\infty} \leq a $ in the optimization described in (6) ?

(c) Computational cost vs performance gain should be examined in both real data and simulations to understand the efficacy of the proposed method.

(d) One thing that comes out in Table 1 tobacco dataset is that the standard error of pairwise pseudo-likelihood is the second highest among all the methods. It may be interesting to check through the simulations if the standard error of pairwise pseudo-likelihood  is in general remain high?
limitations: limitations: The authors have adequately addressed the limitations of the work.
paper_summary: paper_summary: This paper describes a method for estimating missing values of a given non-completed matrix in a situation called MNAR, where the probability distribution for the observation depends on both the matrix index and the matrix value. 

They propose a nonparametric method to prevent bias due to model selection. While missing value estimation in MNAR situations is generally an ill-posed problem, they removed the ill-posedness by imposing appropriate assumptions. This assumption decouples the probability of observation of the $(i,j)$ element of the matrix into an index depending part $t_{ij}$ and a non-negative function $s(Y_{ij})$ of the matrix value $Y_{ij}$, that is, $Pr(i,j)=t_{ij} s(Y_{ij})$, which makes $s$ and $t$ not appear explicitly in the conditional likelihood. 

The paper carefully discusses the difficulties of scale equivalence and shift equivalence that arise when formulating the optimization problem naively under the assumption. The authors also obtain an upper bound on the error of the missing estimation in the Frobenius norm under reasonable assumptions. According to the author's description in L19, methods for estimating missingness in MNAR with such theoretical guarantees are very limited. Finally, they experimentally verified the usefulness of the proposed method on real and synthetic data.
strength_weakness: strength_weakness: -	The biggest concern is that the source code is not available. For example, the paper does not explain how to implement POCS in Algorithm 1. In this situation, there is still a concern about reproducibility without the source code. According to the checklist, the computations can be performed on a personal computer, yet it is better to describe the experimental environment, preferably identifying the operating system, RAM, and programming language.
-	It is worth highlighting that their contribution is to remove ill-posedness through the appropriate assumption. However, it is necessary to define the term "ill-posed" properly and to prove why assumption 1 makes the missing value estimation a well-posed problem. Please let me know if I miss something at this point.
-	Although the experimental results with real data are not promising, this should not be a reason for rejection, given the theoretical guarantees and other strong contributions.

Citations:
- L166: please cite something to support your claim “Nuclear norm regularization has been commonly used to promote low-rankness in the estimation”.
- L100: please cite something to support your claim “Matrix completion under general MNAR is ill-posed”.

Notation:
- For probability distribution, the authors use multiple symbols such as $P$, $Pr$, $p$, and $\mathbb{P}$. I did not get how the authors make the difference among them. 
- Some symbols are not defined clearly. 
   - L161: I could not find the definition of the inner product between two matrices.
   - L177: What is the symbol “N”? Also, in L178, the authors said, “phi is the pdf,” but I think phi is the dispersion parameter. Please let me know if I misunderstood something.
   - L199: It is better to state clearly that $[n]$ is the set of natural numbers less than or equal to a natural number $n$.

Minor comments:
- What is the definition of the non-negative function in L106? $s:\mathbb{R}\rightarrow\mathbb{R}^+$? If so, please write it clearly. 
- The computational complexity of the proposed method should be described.
- L82: In my understanding, $Y_{ij}$ is assumed to be sampled from $f_{ij}(y \mid A, \phi)$. So, if they write $Y_{ij} \sim f_{ij}(y \mid A, \phi) = \dots $ in Equation (1), it is very kind to readers.
- L196: It would be kind to note again that $m_1$ and $m_2$ are the sizes of the matrix A.
- L656: $T_{i,j}$ → $T_{ij}$, $Y_{i,j}$ → $Y_{ij}$
questions: questions: - L88: I would like to confirm the meaning of “independent” here. You mean $Y_{ij}$ and $T_{ij}$ are not independent of each other, but the tuple $(Y_{ij}, T_{ij})$ is independent of $(Y_{kl}, T_{kl})$ for all $k$ and $l$. Is that right?

- In Appendix B, why do they use the symbol $a(y)$ instead of $s(y)$?

[Edit: after rebuttal]
As the authors addressed my comments appropriately, I raised my score from 6 to 7.
limitations: limitations: The authors adequately addressed the limitations related to shift equivalence and scale equivalence in Section 3.

Comments: title: Comments to the rebuttal
comment: Thank you for your feedback and support! We appreciate your acknowledgment of our efforts in rebuttal. Your constructive feedback is invaluable in refining our work.
I am happy with the authors' detailed responses to my queries. Most of the theoretical issues have been clarified and the authors provided a detailed numerical comparison with the existing methods in terms of computational times.

I think they have done a great job in answering in detail all the other queries raised by the reviewers. I am happy to keep the score at 7 and would support the "acceptance" of the paper.
We are glad we could provide the clarity you needed. Thank you again for your comments and feedback, which helped improve the quality of our work!
I am grateful for the meticulous reply from the author. Consequently, I have decided to upgrade my rating from 6 to 7
Dear Reviewer CJYv,

Thanks for your appreciation! We will be sure to clarify ill-posedness and clean up all the notations in the revision.

Best,

Authors of Submission 14957
Dear authors, 

I appreciate your informative responses. It makes sense that the problem is still ill-posed. As I understand it, a problem is ill-posed when either (i) there is no solution, (ii) the solution is not unique, or (iii) the solution is not continuous to the initial value. However, the usage of the term varies from community to community. Although ill-posedness seems to imply non-identifability in statistics, I would encourage the authors to write more clearly that ill-posed means non-identifability in the main text.

Now, the source code is submitted, and the provided insights for the experiments are reasonable. I am totally satisfied with your response. The paper is well-written and will be more readable after you will fix some notations. Thus, I raise my score from 6 to 7.

Thanks again for your kind responses and for considering fixing the issue of notations.