Binary matrix completion with nonconvex regularizers

TL;DR

This paper introduces a novel binary matrix completion model using nonconvex regularizers, providing theoretical recovery guarantees and an accelerated algorithm, outperforming existing convex regularization methods on synthetic and real data.

Contribution

The paper proposes a new nonconvex regularizer-based BMC model with recovery guarantees and an accelerated proximal algorithm for efficient solution.

Findings

The proposed method outperforms convex regularization approaches in experiments.

The accelerated algorithm achieves a convergence rate of 1/T.

The model demonstrates superior recovery performance on synthetic and real datasets.

Abstract

Many practical problems involve the recovery of a binary matrix from partial information, which makes the binary matrix completion (BMC) technique received increasing attention in machine learning. In particular, we consider a special case of BMC problem, in which only a subset of positive elements can be observed. In recent years, convex regularization based methods are the mainstream approaches for this task. However, the applications of nonconvex surrogates in standard matrix completion have demonstrated better empirical performance. Accordingly, we propose a novel BMC model with nonconvex regularizers and provide the recovery guarantee for the model. Furthermore, for solving the resultant nonconvex optimization problem, we improve the popular proximal algorithm with acceleration strategies. It can be guaranteed that the convergence rate of the algorithm is in the order of ${1/T}$, where $T$ is the number of iterations. Extensive experiments conducted on both synthetic and real-world data sets demonstrate the superiority of the proposed approach over other competing methods.

Tables

Table 1. Table 1: The functions r 𝑟 r and thresholds γ 𝛾 \gamma for nonconvex regularizers where μ > 0 𝜇 0 \mu>0 , η = λ ρ 𝜂 𝜆 𝜌 \eta={\textstyle{\lambda\over\rho}} . For the TNN regularizer, μ 𝜇 \mu is an integer denoting the number of leading singular values that are not penalized.

	$\eta r\left({\sigma }_i\left(𝐗\right)\right)$	$\gamma$
TNN	$\{\begin{array}{c}0,i\le \mu \\ \eta {\sigma }_i\left(𝐗\right),i>\mu \end{array}$	$\max ({\sigma }_{\mu +1}\left(𝐗\right),\eta )$
capped ${\ell }_1$	$\{\begin{array}{c}\eta {\sigma }_i\left(𝐗\right),{\sigma }_i\left(𝐗\right)\le \mu \\ \eta \mu ,{\sigma }_i\left(𝐗\right)>\mu \end{array}$	$\min (\eta ,\sqrt{2\mu \eta })$
LSP	$\eta \log \left({\sigma }_i\left(𝐗\right)/\mu +1\right)$	$\min (\eta /\mu ,\mu )$

Table 2. Table 2: Summary of Recommender system data sets.

	#users	#items	#ratings
MovieLens-100K	943	1,682	100,000
FlimTrust	1,508	2,071	35,497
Douban	3,000	3,000	136,891

Equations

$$\begin{array}[]{l}P\left({\left.{{A_{ij}}=0}\right|{M_{ij}}=1}\right)=1-\delta\\ P\left({\left.{{A_{ij}}=1}\right|{M_{ij}}=0}\right)=0\end{array}.$$

$$\underset{\mathbf{X}\in {\left\{0,1\right\}}^{m\times n}}{\min }F\left(\mathbf{X}\right)\equiv \ell \left(\mathbf{X}\right)+\lambda r_n\left(\mathbf{X}\right).$$

$${\ell }_{\omega }\left(x,a\right)=\omega I_{a=a_1}\ell \left(x,a_1\right)+\left(1-\omega \right)I_{a=a_2}\ell \left(x,a_2\right).$$

$$\underset{\mathbf{X},\mathbf{A}\in {\left\{0,1\right\}}^{m\times n}}{\min }F\left(\mathbf{X}\right)\equiv \lambda r_n\left(\mathbf{X}\right)+{\sum }_{i,j}{\ell }_{\omega }\left(X_{ij},A_{ij}\right).$$

$$R\left(\mathbf{X}\right)=\frac{1}{mn}{\Vert \mathbf{X}-\mathbf{M}\Vert }_F^2.$$

$$U_{\omega }\left(x,a\right)=\left(1-\omega \right)I_{x=1}{I}_{a=0}+\omega I_{x=0}{I}_{a=1}.$$

$$R_{\omega }\left(\mathbf{X}\right)=E\left[{\sum }_{i,j}{U}_{\omega }\left(X_{ij},A_{ij}\right)\right].$$

$$\begin{array}[]{l}P\left({\left.{{A_{ij}}=0}\right|{M_{ij}}=1}\right)=1-\delta={\rho_{+1}}\\ P\left({\left.{{A_{ij}}=1}\right|{M_{ij}}=0}\right)=0={\rho_{-1}}\end{array}.$$

$$\frac{1}{mn}{\Vert \hat{\mathbf{X}}-\mathbf{M}\Vert }_F^2\le \frac{C}{{\delta }^2}\left(\sqrt{\frac{\log \left(2/\alpha \right)}{mn}}+b\right).$$

$$\begin{array}[]{l}{L_{\omega}}\left({\mathbf{X},\mathbf{A}}\right)\\ =\left({1-\omega}\right)\sum\limits_{{A_{ij}}=0}{{{\left({{X_{ij}}-{A_{ij}}}\right)}^{2}}}+\omega\sum\limits_{{A_{ij}}=1}{{{\left({{X_{ij}}-{A_{ij}}}\right)}^{2}}}\\ =\left({1-\omega}\right)\left\|{\mathbf{X}-\mathbf{A}}\right\|_{F}^{2}+\left({2\omega-1}\right)\left\|{{{\cal P}_{\Omega}}\left({\mathbf{X}-\mathbf{A}}\right)}\right\|_{F}^{2}\end{array}.$$

$${\mathbf{X}}_{t+1}=pro{x}_{\frac{\lambda }{\rho }{f}_2}\left({\mathbf{X}}_t-\frac{1}{\rho }\nabla f_1\left({\mathbf{X}}_t\right)\right).$$

$$pro{x}_{\tau {\Vert \mathbf{X}\Vert }_{\ast }}\left(\mathbf{Z}\right)=\mathbf{U}{\left(\mathbf{\Sigma }-\tau \mathbf{I}\right)}_{+}{\mathbf{V}}^{\top }.$$

$${\mathbf{X}}_{t+1}=pro{x}_{\frac{\lambda }{\rho }{r}_n}\left({\mathbf{X}}_t-\frac{1}{\rho }\nabla L_{\omega }\left({\mathbf{X}}_t,\mathbf{A}\right)\right).$$

$$\frac{1}{2}\nabla L_{\omega }\left(\mathbf{X},\mathbf{A}\right)=\left(1-\omega \right)\left(\mathbf{X}-\mathbf{A}\right)+\left(2\omega -1\right){\mathcal{P}}_{\Omega }\left(\mathbf{X}-\mathbf{A}\right).$$

$$pro{x}_{r_n}\left(\mathbf{Z}\right)=\mathbf{U}diag\left(\hat{s}\right){\mathbf{V}}^{\top }.$$

$${\hat{s}}_i\in \arg \underset{s_i\ge 0}{\min }\frac{1}{2}{\left(s_i-{\sigma }_i\left(\mathbf{Z}\right)\right)}^2+\frac{\lambda }{\rho }r\left(s_i\right).$$

$$pro{x}_{r_n}\left(\mathbf{Z}\right)=\mathbf{W}pro{x}_{r_n}\left(\mathbf{W}{}^{\top }\mathbf{Z}\right).$$

$${\mathbf{X}}_{t+1}=\mathbf{W}{\mathbf{U}}_{a}diag(\hat{s}){\mathbf{V}}_a^{\top }.$$

$$\underset{t=1,\cdots ,T}{\min }{\Vert {\mathbf{X}}_{t+1}-{\mathbf{X}}_t\Vert }_F^2\le \frac{2}{\left(\rho -\beta \right)T}\left(F\left({\mathbf{X}}_1\right)-\inf F\right).$$

$${\mathbf{Z}}_{ig}=\left\{c_1{\mathbf{X}}_t+c_2{\mathbf{X}}_{t-1}+c_3\mathbf{A}\right\}+c_4{\mathcal{P}}_{\Omega }\left({\mathbf{Z}}_t-\mathbf{A}\right).$$

Full text

Binary matrix completion, Link prediction, Nonconvex regularizers, Topology inference

1 Introduction

The matrix completion problem attempts to recover a low-rank or an approximate low-rank matrix by observing only partial elements [1]. In recent years, many strong theoretical analyses have been developed on the matrix completion problem [2, 3, 4, 5, 6, 7], which has been applied to a wide variety of practical applications such as background modeling [8, 9], recommender systems [10], sensor localization [11, 12], image and video processing [13, 14], and link prediction [15]. In particular, these all results are based on a potential assumption that the observed entries are continuous-valued. However, in many practical applications, the observations are not only incomplete but also are often highly quantized to a single bit [16]. Therefore, there is a conspicuous gap between those existing approaches and practical situation, which promotes the rapid development of 1-bit matrix completion [16].

Instead of observing a subset of full entries, a more common situation in practice is that only the subset of positive elements can be observed. Thus, the observations are not only binary but also nonnegative. For instance, consider the link prediction problem in social networks, where only positive relationships, such as “friendships”, can be observed, while no “non-friendships” are observed. The goal here is to recover the whole social network from the observed friendships (positive entries). In the context of binary classification, the problems learned from positive and unlabeled examples are called positive and unlabeled learning (PU learning for short) [17]. Consequently, the unobserved entries were regarded as unlabeled samples, and then PU learning is applied to matrix completion [18].

The existing methods of PU matrix completion [18, 19] are all based on the convex regularizers such as nuclear norm and max-norm. However, many works [9, 13, 14, 20] stated that the (convex) nuclear norm might not be a good enough approximation of the rank function. In contrast, better recovery performance can be achieved by nonconvex surrogates [21, 22, 23]. Accordingly, we attempt to introduce the nonconvex regularizers into PU matrix completion.

In this paper, we propose a novel model of PU matrix completion with nonconvex regularizers and provide recovery guarantee for the model. To cope with the challenges of the resultant nonconvex optimization problem, we improve the proximal algorithm with two acceleration schemes: i) Instead of full singular value decomposition (SVD), only a few leading singular values are needed to generate the next iteration. ii) We replace a large matrix by its projection on leading subspace, and then the reduction of matrix size makes the calculation of proximal operator more efficient. Moreover, we show that further acceleration is available by taking advantage of the sparsity structure. Subsequently, the resultant algorithm, named “PU matrix completion with nonconvex regularizers (PUMC_N),” is analyzed in detail from the aspects of convergence and time complexity, respectively.

The primary contributions of our work can be summarized as follows:

• •

Employing the nonconvex regularization, we propose a novel PU matrix completion model and provide a strong guarantee for matrix recovery, i.e., the error in recovering an $m\times n$ 0-1 matrix is ${\cal O}\left({\frac{1}{{{\delta^{2}}\sqrt{mn}}}}\right)$, where $\delta$ denotes the sampling rate of positive entries.

• •

We develop an accelerated version of the proximal algorithm for solving the resultant nonconvex optimization model. It can be guaranteed that the proposed algorithm has a convergence rate of ${\cal O}\left({{1\mathord{\left/{\vphantom{1T}}\right.\kern-1.2pt}T}}\right)$, where $T$ denotes the number of iterations.

• •

We implement and analyze the proposed algorithm on both synthetic and real-world data sets. Experimental results demonstrate the superiority of the resultant algorithm to state-of-the-art methods.

The paper is organized as follows. The following section is a brief overview of the related work. In Section III, we propose the model and provide its recovery guarantee. A fast and efficient algorithm is proposed in Section IV, followed by the convergence and time complexity analysis. Experimental results on both synthetic and real-world data sets are presented in Section V. Finally, the conclusion is summarized in Section VI.

Notation: In this paper, vectors and matrices are denoted by lowercase and uppercase boldface, respectively. For a matrix ${\mathbf{X}}$, ${{\mathbf{X}}^{\top}}$ denotes its transpose, ${{\mathbf{X}}_{k}}={\mathbf{X}}\left({:,1:k}\right)$ is its leading $k$ columns, ${\left\|{\mathbf{X}}\right\|_{F}}=\sqrt{\sum\nolimits_{i,j}{X_{ij}^{2}}}$ is the Frobenius norm of ${\mathbf{X}}$, and ${\left\|{\mathbf{X}}\right\|_{*}}=\sum\nolimits_{i}{{\sigma_{i}}\left({\mathbf{X}}\right)}$ is the nuclear norm, where ${\sigma_{i}}\left({\mathbf{X}}\right)$ is the $i$-th largest singular value of ${\mathbf{X}}$. For a set $\Omega$, $\left|\Omega\right|$ is its cardinal number. In addition, we use $\nabla f$ for the gradient of a differentiable function $f$.

2 Related work

In the last decade, based on the remarkable result of low-rank matrix completion [1], a tremendous amount of work has focused on the problem, which enabled a burst of progress concerning the matrix completion theory. A strong theoretical basis of matrix completion [2, 3, 5], including the case of approximate low-rank matrices and noisy observations [4, 1, 9, 11], has been established.

However, these all results are based on the underlying assumption that the observed entries are continuous-valued. In practice, many applications, such as the popular Netflix111https://netflixprize.com/index.html and MovieLens [16] in recommender systems, have a rating matrix whose entries are discrete and quantized rather than continuous. Consequently, there is a conspicuous gap between standard matrix completion theory and practice, revealing the inadequacy of the corresponding methods in dealing with the above case.

Motivated by the above challenge, 1-bit matrix completion was advocated for the first time in [16] to deal with the binary (1-bit) observations. Theoretical guarantees were provided to show the efficiency of the method. In addition, a suite of experiments on both synthetic and real-world data sets illustrated some of the practice applications and demonstrated the superiority of 1-bit matrix completion. Then, [24] considered a general nonuniform sampling distribution concerning 1-bit matrix completion problem followed by corresponding theoretical guarantees. Moreover, the noisy version was studied in [25] under the same sampling scheme with [24]. Instead of nuclear norm, [25] used the max-norm as a convex relaxation for rank function. Similarly, [26] addressed the problem of social trust prediction with a 1-bit max-norm constrained formulation. Under constraints on infinity norm and exact rank, the noisy 1-bit matrix completion problem was explored in [27] and [28]. Furthermore, though the analysis on PAC-Bayesian bounds, [29] evaluated the performance of 1-bit matrix completion.

Relative to the settings of 1-bit matrix completion, there is a more common situation in practice. Consider the link prediction problem in social networks, instead of observing a subset of full entries, we can only observe a subset of the positive relationships, which is “one-sided” sampling in [18]. The similar situation also occurs in network topology inference problem [48]. In response to such case, [18] proposed PU matrix completion. This method introduced the idea of PU (positive and unlabeled) learning [17, 30], i.e., learning only in the presence of positive and unlabeled examples. Motivated by the development of semi-supervised classification [31] in recent years, [19] proposed the modified version of PU matrix completion.

In particular, the PU matrix completion was considered under the constraints on the nuclear norm. As the tightest convex lower bound of the matrix rank function, the nuclear norm is the most popular convex regularizer. Many algorithms based on the nuclear norm, such as accelerated inexact soft-impute algorithm (AIS-Impute) [32], singular value thresholding (SVT) [33], and inexact augmented Lagrange multipliers (IALM) [34], can solve the corresponding convex optimization problem effectively. Despite the nuclear norm is applied successfully and makes low-rank optimization easier, numerous attempts have recently been made to regard nonconvex regularizers as the better approximation of the matrix rank. For instance, nonconvex surrogates including truncated nuclear norm (TNN) [9, 35], log-sum penalty (LSP) [22, 36], and capped $\ell_{1}$ penalty [21] have been successfully applied in many fields and have better empirical performance than nuclear norm regularizers.

3 Problem formulation

Matrix completion is the problem of recovering the underlying target matrix given its partial information [1]. Following the “basic setting” of [18], let the target matrix ${\mathbf{M}}\in{\left\{{0,1}\right\}^{m\times n}}$ be a binary matrix that consists only of ones and zeros, and ${\Omega_{1}}=\left\{{\left.{\left({i,j}\right)}\right|{M_{ij}}=1}\right\}$ denotes the index set of all positive elements in ${\mathbf{M}}$. Equivalently, the observation matrix is denoted, herein, by ${\mathbf{A}}\in{\left\{{0,1}\right\}^{m\times n}}$, and $\Omega$ denotes the index set of observation elements. According to the “one-sided” sampling in [18], only a subset of positive entries of ${\mathbf{M}}$ can be observed, that is $\Omega\subseteq{\Omega_{1}}$. We suppose that the observation process follows the uniform sampling distribution, which is the popular choice for majority work, i.e., $\Omega$ is sampled randomly from ${\Omega_{1}}$. For the observation matrix ${\mathbf{A}}$, ${A_{ij}}=1$ if $\left({i,j}\right)\in\Omega$ and ${A_{ij}}=0$ otherwise. Here, our goal is to recover the underlying target matrix ${\mathbf{M}}$ from the observation matrix ${\mathbf{A}}$.

According to the above description, the relationship between ${\mathbf{M}}$ and ${\mathbf{A}}$ can be expressed in the following conditional probability.

[TABLE]

where $\delta={{\left|\Omega\right|}\mathord{\left/{\vphantom{{\left|\Omega\right|}{\left|{{\Omega_{1}}}\right|}}}\right.\kern-1.2pt}{\left|{{\Omega_{1}}}\right|}}$ denotes the sampling rate. We consider the problem of positive and unlabeled matrix completion (PU matrix completion) with nonconvex regularizers as the following form.

[TABLE]

where $\lambda$ is a regularization parameter, $\ell$ is a smooth loss function, ${r_{n}}$ is a nonconvex regularizer in Table 1. In addition, (2) has the following characteristics.

• •

$\ell$ is differentiable with $\beta$-Lipschitz continuous gradient, that is, it follows ${\left\|{\nabla\ell\left({{{\mathbf{X}}_{1}}}\right)-\nabla\ell\left({{{\mathbf{X}}_{2}}}\right)}\right\|_{F}}\leq\beta{\left\|{{{\mathbf{X}}_{1}}-{{\mathbf{X}}_{2}}}\right\|_{F}}$, $\beta>0$. Moreover, $\ell$ is bounded from below, i.e., $\inf\ell>-\infty$.

• •

${r_{n}}\left({\mathbf{X}}\right)=\sum\nolimits_{i=1}^{m}{r\left({{\sigma_{i}}\left({\mathbf{X}}\right)}\right)}$ is a nonconvex and nonsmooth function, where $r$ is a nondecreasing concave function and $r\left(0\right)=0$.

• •

${r_{n}}$ can be formulated as the difference of two convex functions [47], i.e., {r_{n}}\left({\mathbf{X}}\right)={\mathord{\buildrel{\hbox{\scriptscriptstyle\frown}}\over{r}}_{n}}\left({\mathbf{X}}\right)-{\mathord{\buildrel{\hbox{\scriptscriptstyle\smile}}\over{r}}_{n}}\left({\mathbf{X}}\right) , where {\mathord{\buildrel{\hbox{\scriptscriptstyle\frown}}\over{r}}_{n}} and {\mathord{\buildrel{\hbox{\scriptscriptstyle\smile}}\over{r}}_{n}} are convex. (The corresponding convex functions of the nonconvex regularizers mentioned in Section II are provided in Appendix A.)

For accurately quantifying the error in recovering the underlying binary matrix, we propose to adopt the $\omega$-weighted square loss [37, 38] as the loss function. The $\omega$-weighted square loss is defined as

[TABLE]

where $\ell\left({x,a}\right)={\left({x-a}\right)^{2}}$ is the square loss, ${I_{a={a_{1}}}}$ and ${I_{a={a_{2}}}}$ are indicator functions, i.e., ${I_{a={a_{1}}}}$ is 1 if $a={a_{1}}$ is true and 0 otherwise.

Consequently, (2) can be further formulated as follows.

[TABLE]

where ${\mathbf{X}}$ and ${\mathbf{A}}$ are the recovery matrix and observation matrix of the underlying target matrix $\mathbf{M}$, respectively, and ${\ell_{\omega}}\left({{X_{ij}},{A_{ij}}}\right)=\omega{I_{{A_{ij}}=1}}\ell\left({{X_{ij}},1}\right)+\left({1-\omega}\right){I_{{A_{ij}}=0}}\ell\left({{X_{ij}},0}\right)$.

Recovery error of (4). Following the definition of the recovery error in [16, 18], here the recovery error can be formulated as

[TABLE]

where ${\mathbf{M}},{\mathbf{X}}\in{\mathbb{R}^{m\times n}}$ is the underlying target matrix and its recovery matrix, respectively.

In [38], the label-dependent loss is defined as $U\left({x,a}\right)={I_{x=1}}{I_{a=0}}+{I_{x=0}}{I_{a=1}}$ . And similar to (3), we define the weighted version of the label-dependent loss as

[TABLE]

Therefore, the corresponding $\omega$-weighted expected error can be written as

[TABLE]

According to the class-conditional random noise model in [39], (1) can be written correspondingly as

[TABLE]

Theorem 1. For the choices $\hat{\omega}=\frac{{1-{\rho_{+1}}}}{2}$ and $\kappa=\frac{{1-{\rho_{+1}}}}{2}$, there exists a constant $c$ that is dependent of ${\mathbf{X}}$ such that, for any matrix ${\mathbf{X}}$, we have ${R_{\hat{\omega}}}\left({\mathbf{X}}\right)=\kappa R\left({\mathbf{X}}\right)+c$.

The above theorem (Theorem 1) is a special case of Theorem 9 in [39]. At this juncture, the linear mapping between the recovery error $R\left({\mathbf{X}}\right)$ and $\omega$-weighted expected error ${R_{\hat{\omega}}}\left({\mathbf{X}}\right)$ indicates that minimizing $R\left({\mathbf{X}}\right)$ is equivalent to minimizing ${R_{\hat{\omega}}}\left({\mathbf{X}}\right)$ on the partial information.

Theorem 2 (Main Result 1). Let ${\mathbf{\hat{X}}}\in{\mathbb{R}^{m\times n}}$ be the solution to (4), then with probability at least $1-\alpha$,

[TABLE]

where $C$ is absolute constant, $\delta$ denotes the sampling rate, $b=\frac{{\sqrt{m}+\sqrt{n}+\sqrt[4]{{\left|{{\Omega_{1}}}\right|}}}}{{mn}}$, $\Omega_{1}$ is the index set of all positive elements in ${\mathbf{M}}$. The proof can be found in Appendix C.

4 Algorithm

In this section, we will show the nonconvex model can be solved much faster. First, Subsection A shows the basic algorithm for nonconvex regularizers, followed by acceleration measures in Subsection B. Subsection C summarizes the whole algorithm and Subsection D analyses the resultant algorithm from the aspect of convergence and time complexity.

4.1 Basic algorithm

Let ${L_{\omega}}\left({{\mathbf{X}},{\mathbf{A}}}\right)=\sum\nolimits_{i,j}{{\ell_{\omega}}\left({{X_{ij}},{A_{ij}}}\right)}$, and in line with the definition of the Frobenius norm, we have the following derivation.

[TABLE]

where ${\left[{{{\cal P}_{\Omega}}\left({{\mathbf{X}}-{\mathbf{A}}}\right)}\right]_{ij}}={\left({{\mathbf{X}}-{\mathbf{A}}}\right)_{ij}}$ if $\left({i,j}\right)\in\Omega$ and 0 otherwise. In recent years, the proximal algorithm [32] has been regarded as an efficient method for solving the optimization problem $\mathop{\min}\limits_{\mathbf{X}}{f_{1}}\left(\mathbf{X}\right)+{f_{2}}\left(\mathbf{X}\right)$, when $f_{1}$ and $f_{2}$ are convex. The following theorem shows that the convergence of the proximal algorithm.

Theorem 3 [40]. Let $f_{1}$, $f_{2}$ be lower semicontinuous, and $f_{1}$ is differentiable with $\beta$-Lipschitz continuous gradient. If $f_{1}+f_{2}$ is coercive and strictly convex, the solution of the problem takes on uniqueness. For an arbitrary initial matrix $\mathbf{X}$, $\forall\rho\geq\beta$, the iterative sequence generated by the following statement can converge to the unique solution of the problem.

[TABLE]

where $pro{x_{{\textstyle{\lambda\over\rho}}{f_{2}}}}\left({\mathbf{Z}}\right)=\mathop{\arg\min}\limits_{{\mathbf{X}},{\mathbf{Z}}\in{\mathbb{R}^{m\times n}}}\left\{{\frac{1}{2}\left\|{{\mathbf{X}}-{\mathbf{Z}}}\right\|_{F}^{2}+{\textstyle{\lambda\over\rho}}{f_{2}}\left({\mathbf{X}}\right)}\right\}$ denotes the proximal operator.

If $f_{2}$ is the nuclear norm, the following theorem shows that the proximal operator of the nuclear norm has a closed form solution.

Theorem 4 [33]. For an arbitrary matrix ${\mathbf{Z}}\in{\mathbb{R}^{m\times n}}$, $\forall\tau>0$, the proximal operator of the nuclear norm of matrix ${\mathbf{X}}$ is

[TABLE]

where $\mathbf{I}$ denotes the identity matrix, $SVD\left({\mathbf{Z}}\right)={\mathbf{U\Sigma}}{{\mathbf{V}}^{\top}}$, and ${\left[{{\mathbf{M}_{+}}}\right]_{ij}}=max\left({{M_{ij}},0}\right)$.

For solving (4), we extend the proximal operator to nonconvex problem, similar to Theorem 3, at $t$-iteration, it products the iterative sequence as follows.

[TABLE]

where the learning rate, herein, denoted by $\rho$ is a fixed value, and $\nabla{L_{\omega}}\left({{\mathbf{X}},{\mathbf{A}}}\right)$ denotes the gradient of the $\omega$-weighted loss function, which can be computed efficiently as

[TABLE]

Recently, due to the successful application on convex optimization problem, the proximal algorithm has been extended to nonconvex situation [9, 13, 14, 20]. Similar to the nuclear norm (Theorem 4), the generalized singular value thresholding [20] was proposed to handle the nonconvex surrogates.

Theorem 5 Generalized singular value thresholding (GSVT) [20]. For an arbitrary matrix ${\mathbf{Z}}\in{\mathbb{R}^{m\times n}}$, let $r_{n}$ be a function that satisfies the characteristics in (2), then the proximal operator of $r_{n}$ has the following closed form solution.

[TABLE]

where ${\mathbf{U}}{\mathbf{\Sigma}}{{\mathbf{V}}^{\top}}$ is the SVD of $\mathbf{Z}$, and $\hat{s}=\left\{{{{\hat{s}}_{i}}}\right\}$ with

[TABLE]

Similar to Theorem 4, the above theorem indicates that the closed-form solutions of the nonconvex regularizer in Table 1 do exist. In addition, we generalize the above procedure as the Basic Algorithm shown in Algorithm 1.

4.2 Acceleration

However, the basic algorithm involves a full SVD (step 3) with the time complexity ${{\cal O}}\left({m{n^{2}}}\right)$. Next, we will take the following two schemes to make the basic algorithm much faster.

S1. The first thing that comes to our mind is to use partial SVD. However, $\hat{s}_{i}$ in (16) actually becomes zero when the singular value ${\sigma_{i}}\left({\mathbf{Z}}\right)$ is not larger than a threshold $\gamma$, that is automatic thresholding in [41]. This means only the leading few singular values, instead of all singular values, are needed to compute the proximal operator in Theorem 5. The thresholds $\gamma$ for the mentioned nonconvex regularizers are shown in Table 1.

S2. Another measure our carry out is reducing the size of SVD. The following theorem shows that the proximal operator on a large matrix can be replaced by the counterpart on a smaller one. The proof can be found in Appendix D.

Theorem 6. For an arbitrary matrix $\mathbf{Z}\in{\mathbb{R}^{m\times n}}$, let $k$ be the number of singular values of $\mathbf{Z}$ that are not less than $\gamma$, its rank-$k$ SVD is ${{\mathbf{U}}_{k}}{{\mathbf{\Sigma}}_{k}}{\mathbf{V}}_{k}^{\top}$, and $\mathbf{W}\in{\mathbb{R}^{m\times k}}$ be an orthogonal matrix, if $span\left({{{\mathbf{U}}_{k}}}\right)\subseteq span\left({\mathbf{W}}\right)$ ,then the following equation holds.

[TABLE]

Theorem 6 performs the process of replacing a large matrix by its projection on leading subspace. How to obtain such $\mathbf{W}$ in Theorem 6? Two approaches are available to find the $\mathbf{W}$ exactly in the same time complexity. The first method, PROPACK package [42], is widely applied to partial SVD. And the second one is the power method which has good approximation guarantee [43]. Compared with the former method, the latter one can benefit particularly from warm-start taking full advantage of the iterative nature of the proximal algorithm. Hence, we use the power method to get the $\mathbf{W}$, and the details are shown in Algorithm 2.

Let ${{\mathbf{X}}_{ig}}={{{\mathbf{X}}_{t}}-{\textstyle{1\over\rho}}\nabla{L_{\omega}}\left({{{\mathbf{X}}_{t}},{\mathbf{A}}}\right)}$, through the implementation of the above two acceleration measures, (13) can be rewritten as

[TABLE]

where ${{\mathbf{U}}_{a}}{{\mathbf{\Sigma}}_{a}}{\mathbf{V}}_{a}^{\top}$ is the rank-$a$ SVD of ${\mathbf{W}^{\top}{\mathbf{X}}_{ig}}$,$a$ is the number of singular values that are greater than the threshold $\gamma$, and $\hat{s}$ can be obtained from (16).

4.3 The whole algorithm

We summarize the whole procedure for solving (4) in Algorithm 3 and name it PUMC_N (Positive and Unlabeled Matrix Completion with Nonconvex regularizers). In step 3, similar to the nmAPG algorithm [44], a linear combination of ${{\mathbf{X}}_{t-1}}$ and ${{\mathbf{X}}_{t}}$ is used to accelerate the algorithm. The column spaces of the current iteration (${{\mathbf{V}}_{t}}$) and previous iteration (${{\mathbf{V}}_{t-1}}$) are used to accomplish the warm start in step 5 as in [32]. Step 6 and step 7 performs S2, and in step 8, a continuation strategy is introduced to speed up the algorithm further. Specifically, as the iteration proceeds, $\lambda$ is dynamic and gradually decreases from a large value. In addition, step 9-11 performs S1.

4.4 Algorithm analysis

Convergence analysis. Firstly, we present a lemma which provides the basic support for the convergence analysis of the proposed algorithm. The following lemma shows that the objective function $F$ is nonincreasing as iterations proceed.

Lemma 1 [45]. Let $\left\{{{{\mathbf{X}}_{t}}}\right\}$ be the iterative sequence produced by (13), for the optimization problem (2), we have $F\left({{{\mathbf{X}}_{t+1}}}\right)\!\leq\!F\left({{{\mathbf{X}}_{t}}}\right)\!-\!{\textstyle{{\rho-\beta}\over 2}}\left\|{{{\mathbf{X}}_{t+1}}\!-\!{{\mathbf{X}}_{t}}}\right\|_{F}^{2}$, where $\rho\!>\!\beta$.

The following theorem shows that the proposed algorithm generates a bounded iterative sequence. The proof can be found in Appendix E.

Theorem 7. Let $\left\{{{{\mathbf{X}}_{t}}}\right\}$ be the iterative sequence produced by (13), we say $\left\{{{{\mathbf{X}}_{t}}}\right\}$ is a bounded iterative sequence, i.e., $\sum\nolimits_{t=1}^{\infty}{\left\|{{{\mathbf{X}}_{t+1}}-{{\mathbf{X}}_{t}}}\right\|_{F}^{2}}<\infty$.

For the convex optimization problem in Theorem 3, the proximal mapping in [45] is denoted by ${G_{{\scriptstyle{1\over\rho}}{f_{2}}}}\left({{{\mathbf{X}}_{t}}}\right)=pro{x_{{\scriptstyle{1\over\rho}}{f_{2}}}}\left({{{\mathbf{X}}_{t}}-{\textstyle{1\over\rho}}\nabla{f_{1}}\left({{{\mathbf{X}}_{t}}}\right)}\right)-{{\mathbf{X}}_{t}}$. In particular, when $f_{2}$ is convex, ${\left\|{{G_{{\scriptstyle{1\over\rho}}{f_{2}}}}\left({{{\bf{X}}_{t}}}\right)}\right\|_{2}^{2}}$ can be used to conduct the convergence analysis. In contrast, if $f_{2}$ is nonconvex, it is no longer applicable. Hence, we use $\left\|{{G_{{\scriptstyle{1\over\rho}}{f_{2}}}}\left({{{\mathbf{X}}_{t}}}\right)}\right\|_{F}^{2}{\rm{=}}\left\|{{{\mathbf{X}}_{t+1}}-{{\mathbf{X}}_{t}}}\right\|_{F}^{2}$ instead to perform convergence analysis of the proposed algorithm. The convergence of Algorithm 3 is shown in the following theorem, and the proof can be found in Appendix F.

Theorem 8 (Main Result 2). Let $\left\{{{{\mathbf{X}}_{t}}}\right\}$ be the iterative sequence in Algorithm 3, for the consecutive elements ${\mathbf{X}_{t}}$ and ${\mathbf{X}_{t+1}}$, we have

[TABLE]

Time complexity. Assume that $\mathbf{Y}_{t}$ in step 5 of Algorithm 3 has $k_{t}$ columns at the current iteration. Consequently, step 5 takes ${{\cal O}}\left({nk_{t}^{2}}\right)$ time. Next, step 3 shows that $\mathbf{Z}_{t}$ is a linear combination of $\mathbf{X}_{t-1}$ and $\mathbf{X}_{t}$. Let ${c_{t}}={{\left({{\alpha_{t-1}}-1}\right)}\mathord{\left/{\vphantom{{\left({{\alpha_{t-1}}-1}\right)}{{\alpha_{t}}}}}\right.\kern-1.2pt}{{\alpha_{t}}}}$, combining step 4 and step 5, we have

[TABLE]

where ${c_{1}}=\left({1+{c_{t}}}\right)\left({1-{c_{3}}}\right)$, ${c_{2}}={c_{t}}\left({{c_{3}}-1}\right)$, ${c_{3}}={\textstyle{{2\left({1-\omega}\right)}\over\rho}}$, ${c_{4}}={\textstyle{{2\left({1-2\omega}\right)}\over\rho}}$. The first three terms involve low-rank matrices, whereas the last term involves a sparsity structure. The combined structure in (20) was studied specifically in [42]. Consider the multiplication of $\mathbf{Z}_{ig}$ and a vector $\mathbf{b}\in{\mathbb{R}^{n}}$. For the low-rank part the multiplication cost ${{\cal O}}\left({\left({m+n}\right){k_{t}}}\right)$ time, whereas the sparse part cost ${{\cal O}}\left({{{\left\|\Omega\right\|}_{1}}}\right)$ time. Hence, the cost is ${{\cal O}}\left({\left({m+n}\right){k_{t}}+{{\left\|\Omega\right\|}_{1}}}\right)$ per vector multiplication. Step 8 performs a rank-$k_{t}$ SVD of ${{\mathbf{W}}^{\top}}{{\mathbf{Z}}_{ig}}$, it takes ${{\cal O}}\left({\left({m+n}\right)k_{t}^{2}+{{\left\|\Omega\right\|}_{1}}{k_{t}}}\right)$ time. In summary, the order of the time complexity at the current iteration is ${{\cal O}}\left({\left({m+n}\right)k_{t}^{2}+{{\left\|\Omega\right\|}_{1}}{k_{t}}}\right)$.

5 Experiments

In this section, we perform experiments on synthetic and real-world data sets and demonstrate the effectiveness of the proposed algorithm in practical applications, including link prediction, topology inference, and recommender system. All experiments are implemented in Matlab on Windows 10 with Intel Xeon CPU (2.8GHz) and 128GB memory.

5.1 Synthetic Data

Data sets. As in [18, 41], we assume the matrix $\mathbf{Q}\in{\mathbb{R}^{m\times m}}$ is generated by $\mathbf{Q}=\mathbf{M}_{1}\mathbf{M}_{2}$, where the elements of $\mathbf{M}_{1}\in{\mathbb{R}^{m\times k}}$ and $\mathbf{M}_{2}\in{\mathbb{R}^{k\times m}}$ are obtained from the Gaussian distribution ${\cal N}\left({0,1}\right)$. Therefore, the underlying binary matrix $\mathbf{M}\in{\mathbb{R}^{m\times m}}$ can be generated by ${M_{ij}}={{I}_{{Q_{ij}}\geq q}}$, without loss of generality, we assume that $q=0.5$. We fix $k=5$ and very $m$ in $\left\{{50,100,500,1000,2000}\right\}$.

For each scenario, the mean square error ${\rm{MSE}}={{\left\|{{{\cal P}_{\bar{\Omega}}}\left({{\mathbf{X}}-{\mathbf{M}}}\right)}\right\|_{F}^{2}}\mathord{\left/{\vphantom{{\left\|{{{\cal P}_{\bar{\Omega}}}\left({{\bf{X}}-{\mathbf{M}}}\right)}\right\|_{F}^{2}}{\left\|{{{\cal P}_{\bar{\Omega}}}\left({\mathbf{M}}\right)}\right\|_{F}^{2}}}}\right.\kern-1.2pt}{\left\|{{{\cal P}_{\bar{\Omega}}}\left({\mathbf{M}}\right)}\right\|_{F}^{2}}}$ is used for performance evaluation, where $\mathbf{X}$ is the recovery matrix for underlying matrix $\mathbf{M}$ and $\bar{\Omega}$ is the index set of unobserved elements. Each experiment is repeated ten times with the sampling rate ($\delta$) varying from 0.3 to 0.9, and the average results are reported. From Theorem 1, if the sampling rate $\delta=0.3$ (only 30% 1’s in $\mathbf{M}$ are observed), $\omega=0.15$ is chosen.

Results are shown in Fig. 1. Only TNN regularizer (with $\mu$ in Table 1 set to 5) is used in synthetic experiments, similar results can be obtained from other nonconvex regularizers in Table 1. Fig. 1 shows the testing MSE at different sampling rates. Notice that for an increasing sampling rate we see a monotonous decrease in testing MSE until it is close to zero. This is reasonable since a larger number of observations give rise to more accurate information of the underlying matrix. In addition, there is also a negative correlation between MSE and the matrix size, which is particularly evident at a smaller sampling rate. In particular, if the underlying matrix $\mathbf{M}$ is large enough, it can be recovered accurately at a much small sampling rate (with small number of observations). In Fig. 1, we follow the settings in [41] and fix the sampling rate at 0.5 (from Theorem 1, $\omega=0.25$). It can be seen that the MSE drops sharply and precipitously at the beginning. Moreover, the larger the matrix, the more time the algorithm takes, and the smaller the MSE.

5.2 Recommender System

Data sets. In this practical setting, the popular data sets (Table 2), including MovieLens (100K)222http://vladowiki.fmf.uni-lj.si/doku.phpid=pajek:data:pajek:students [16], * FlimTrust333https://www.librec.net/datasets.html* [53], and Douban444https://github.com/fmonti/mgcnn [54], are used to evaluate the performance of our algorithm. We follow the setup in [16, 19] and convert these ratings in each data set to binary observations by comparing each rating to the average value (which is $\sim 3$, considering three data sets together) of whole data sets.

Methods. We compare with the nuclear norm based algorithm AIS-Impute [32], as well as 1-bit matrix completion in [19]. In particular, the AIS-Impute can be considered as an accelerated and inexact version of the proximal algorithm, and the 1-bit matrix completion is constrained by infinity norm and nuclear norm.

Results are shown in Fig. 2. Each point in the figure is the average across ten replicate experiments. Moreover, Table LABEL:table6 shows the performance of the mentioned methods on three data sets. From the above experiments, LSP regularizer usually has better or comparable performance than the other two regularizers, thus we only use LSP regularizer here. It can be seen that in the three recommended algorithms, as long as the sampling rate is not less than 0.4, PUMC_N will result in the lowest MSE in the least time. In addition, when the sampling rate is low, the performance of PUMC_N needs to be improved. Moreover, the MSE vs time on the above three data sets is provided in Appendix B.

6 Conclusion

In this paper, we addressed the problem of binary matrix completion with nonconvex regularizers, where the observations consist only of positive entries. We proposed a novel PU matrix completion model (4) for tacking the task based on the commonly-used nonconvex regularizers and the $\omega$-weighted loss. In particular, the error bound for the model is derived to show the underlying matrix ${\mathbf{M}}\in{\left\{{0,1}\right\}^{m\times n}}$ can be recovered accurately. Accordingly, we improved the proximal algorithm with two main acceleration strategies in nonconvex settings for solving (4), and the convergence can also be guaranteed. The experiments on both synthetic and real-world data sets have verified the effectiveness of the proposed approach and validated the superiority over the state-of-the-art methods.

There still remain several directions for further work. From the experimental results, it can be seen that there is still room for further performance improvements at a low sampling rate. Besides, as in [24, 25], a general nonuniform sampling distribution will be considered. In addition, to further speed up the proposed algorithm and apply it to massive data sets, we will focus on its distributed version.