Principled analytic classifier for positive-unlabeled learning via weighted integral probability metric

TL;DR

This paper introduces a scalable, theoretically justified positive-unlabeled learning algorithm that produces a closed-form classifier, with improved error bounds and demonstrated superior performance on synthetic and real datasets.

Contribution

The paper presents a novel PU learning algorithm with a closed-form solution and sharper error bounds, enhancing scalability and robustness over existing methods.

Findings

Achieves improved accuracy on synthetic and real datasets

Offers a closed-form classifier in RKHS hypothesis space

Provides sharper estimation error bounds and explicit excess risk bounds

Abstract

We consider the problem of learning a binary classifier from only positive and unlabeled observations (called PU learning). Recent studies in PU learning have shown superior performance theoretically and empirically. However, most existing algorithms may not be suitable for large-scale datasets because they face repeated computations of a large Gram matrix or require massive hyperparameter optimization. In this paper, we propose a computationally efficient and theoretically grounded PU learning algorithm. The proposed PU learning algorithm produces a closed-form classifier when the hypothesis space is a closed ball in reproducing kernel Hilbert space. In addition, we establish upper bounds of the estimation error and the excess risk. The obtained estimation error bound is sharper than existing results and the derived excess risk bound has an explicit form, which vanishes as sample sizes increase. Finally, we conduct extensive numerical experiments using both synthetic and real datasets, demonstrating improved accuracy, scalability, and robustness of the proposed algorithm.

Tables

Table 1. Table 1: A summary of elapsed training time and its ratio for the five PU learning algorithms based on 100 replications. We set n p = 100 , n u = 400 formulae-sequence subscript 𝑛 p 100 subscript 𝑛 u 400 n_{\rm p}=100,n_{\rm u}=400 , and π + = 0.5 subscript 𝜋 0.5 \pi_{+}=0.5 . Average and standard error are denoted by ‘average ± plus-or-minus \pm standard error’.

	LOG	DH	NNPU	tADJ	WMMD
in seconds $\times 10$	$90.0\pm 4.7$	$96.1\pm 6.1$	$6.0\pm 0.1$	$0.4\pm 0.0$	$0.2\pm 0.0$
in ratio	$347.9\pm 23.4$	$371.0\pm 28.5$	$23.2\pm 1.1$	$1.8\pm 0.0$	$1.0\pm 0.0$

Table 2. Table 2: A summary of the eleven binary classification datasets. ‘# of samples’ denotes the number of total samples after removing incomplete observations. We denote the number of positive, unlabeled, and test samples, by n p subscript 𝑛 p n_{\rm p} , n u subscript 𝑛 u n_{\rm u} , and n te subscript 𝑛 te n_{\rm te} after the random sampling, respectively. We categorize the eleven datasets into two groups: the first seven datasets as small-scale and the last four datasets as large-scale.

Dataset	$d$	# of samples	$n_p$	$n_u$	$n_{\mathrm{te}}$	${\pi }_{+}$	Scale
heart_scale	12	122	10	60	60	0.62	Small
sonar_scale	60	207	10	100	100	0.47	Small
australian_scale	12	449	20	220	220	0.51	Small
australian_scale2	12	449	10	130	130	0.15	Small
breast-cancer_scale	10	683	20	340	340	0.35	Small
breast-cancer_scale2	10	683	40	340	340	0.65	Small
diabetes_scale	8	759	50	380	370	0.65	Small
skin_nonskin	3	245,057	${10}^3$	${10}^5$	${10}^5$	0.79	Large
skin_nonskin2	3	245,057	${10}^3$	${10}^5$	${10}^5$	0.21	Large
epsilon_normalized	2,000	500,000	${10}^3$	$4\times {10}^5$	${10}^5$	0.50	Large
HIGGS	26	8,786,441	${10}^3$	${10}^6$	${10}^5$	0.50	Large

Table 3. Table 3: Accuracy and AUC comparison using the real datasets when the class-prior π + subscript 𝜋 \pi_{+} is known. We denote the memory error results for LOG and DH by the hyphen. Average and standard error are denoted by ‘average ± plus-or-minus \pm standard error’. Boldface numbers denote the best and equivalent algorithms with respect to a t-test with a significance level of 5%.

Dataset	LOG	DH	NNPU	tADJ	WMMD
Accuracy (in %)
heart_scale	$70.5\pm 0.8$	$68.4\pm 0.9$	$71.0\pm 0.8$	$65.1\pm 0.8$	$71.6\pm 0.8$
sonar_scale	$55.8\pm 0.6$	$52.9\pm 0.6$	$63.2\pm 0.6$	$60.7\pm 0.6$	$62.4\pm 0.6$
australian_scale	$85.4\pm 0.4$	$84.9\pm 0.6$	$79.2\pm 0.5$	$80.0\pm 0.7$	$84.2\pm 0.6$
australian_scale2	$85.7\pm 0.2$	$85.7\pm 0.2$	$86.7\pm 0.3$	$75.4\pm 2.0$	$86.2\pm 0.2$
breast-cancer_scale	$95.8\pm 0.1$	$96.0\pm 0.3$	$90.1\pm 0.3$	$91.0\pm 0.4$	$89.3\pm 0.5$
breast-cancer_scale2	$95.6\pm 0.3$	$94.4\pm 0.8$	$95.9\pm 0.1$	$92.5\pm 0.2$	$94.2\pm 0.3$
diabetes_scale	$66.7\pm 0.7$	$65.5\pm 0.9$	$69.4\pm 0.4$	$67.9\pm 0.3$	$66.4\pm 0.2$
skin_nonskin	-	-	$98.2\pm 0.1$	$78.0\pm 0.4$	$85.3\pm 0.7$
skin_nonskin2	-	-	$98.6\pm 0.0$	$93.9\pm 0.1$	$98.1\pm 0.2$
epsilon_normalized	-	-	$64.5\pm 0.3$	$63.1\pm 0.1$	$56.3\pm 1.3$
HIGGS	-	-	$56.3\pm 0.2$	$52.6\pm 0.1$	$54.0\pm 0.2$
AUC $\times 100$
heart_scale	$78.4\pm 1.0$	$78.3\pm 1.1$	$73.8\pm 0.9$	$72.5\pm 1.1$	$79.0\pm 0.9$
sonar_scale	$61.2\pm 0.8$	$60.6\pm 0.9$	$67.4\pm 0.7$	$66.2\pm 0.7$	$68.9\pm 0.8$
australian_scale	$91.1\pm 0.2$	$91.3\pm 0.3$	$87.8\pm 0.4$	$87.8\pm 0.5$	$90.4\pm 0.4$
australian_scale2	$89.2\pm 0.4$	$87.3\pm 0.4$	$84.3\pm 0.6$	$85.9\pm 0.7$	$88.6\pm 0.6$
breast-cancer_scale	$99.4\pm 0.0$	$99.3\pm 0.0$	$97.8\pm 0.1$	$95.6\pm 0.4$	$99.5\pm 0.0$
breast-cancer_scale2	$99.3\pm 0.0$	$99.2\pm 0.1$	$99.3\pm 0.0$	$97.2\pm 0.2$	$98.7\pm 0.2$
diabetes_scale	$74.0\pm 0.6$	$71.5\pm 1.1$	$73.5\pm 0.6$	$74.7\pm 0.5$	$74.5\pm 0.7$
skin_nonskin	-	-	$99.5\pm 0.1$	$94.8\pm 0.1$	$99.4\pm 0.1$
skin_nonskin2	-	-	$99.7\pm 0.0$	$94.6\pm 0.0$	$99.8\pm 0.0$
epsilon_normalized	-	-	$70.0\pm 0.4$	$69.3\pm 0.1$	$62.2\pm 2.3$
HIGGS	-	-	$59.6\pm 0.2$	$65.3\pm 0.1$	$55.7\pm 0.3$

Table 4. Table 4: Accuracy and AUC comparison using the real datasets when the class-prior π + subscript 𝜋 \pi_{+} is unknown. The ‘KM1’ method by Ramaswamy et al. ( 2016 ) is used for LOG, DH, and NNPU, and the density-based method is used for WMMD. We denote the infeasible cases due to ‘KM1’ method by the hyphen. Other details are given in Table 3 .

Dataset	LOG	DH	NNPU	tADJ	WMMD
Accuracy (in %)
heart_scale	$39.5\pm 0.8$	$39.1\pm 0.7$	$42.4\pm 0.9$	$65.1\pm 0.8$	$70.5\pm 0.7$
sonar_scale	$53.5\pm 0.6$	$52.1\pm 0.5$	$59.9\pm 0.7$	$60.7\pm 0.6$	$54.3\pm 0.8$
australian_scale	$50.0\pm 0.2$	$50.0\pm 0.2$	$50.0\pm 0.2$	$80.0\pm 0.7$	$79.4\pm 1.0$
australian_scale2	$84.9\pm 0.2$	$84.9\pm 0.2$	$85.6\pm 0.3$	$85.5\pm 1.0$	$80.2\pm 1.0$
breast-cancer_scale	$65.0\pm 0.2$	$65.2\pm 0.2$	$65.0\pm 0.2$	$91.0\pm 0.4$	$93.0\pm 1.0$
breast-cancer_scale2	$35.0\pm 0.2$	$35.0\pm 0.2$	$35.0\pm 0.2$	$92.5\pm 0.2$	$92.6\pm 0.4$
diabetes_scale	$36.0\pm 0.3$	$41.1\pm 1.1$	$37.9\pm 0.5$	$67.9\pm 0.3$	$65.1\pm 0.2$
skin_nonskin	-	-	-	$78.0\pm 0.4$	$82.2\pm 0.9$
skin_nonskin2	-	-	-	$93.9\pm 0.1$	$95.7\pm 0.4$
epsilon_normalized	-	-	-	$63.1\pm 0.1$	$49.9\pm 0.1$
HIGGS	-	-	-	$52.6\pm 0.1$	$50.9\pm 0.0$
AUC $\times 100$
heart_scale	$67.2\pm 1.8$	$67.2\pm 1.4$	$71.0\pm 0.9$	$72.5\pm 1.1$	$77.1\pm 0.9$
sonar_scale	$60.5\pm 0.9$	$62.2\pm 0.9$	$66.6\pm 0.8$	$66.2\pm 0.7$	$69.7\pm 0.8$
australian_scale	$78.4\pm 1.1$	$72.5\pm 1.4$	$80.3\pm 0.6$	$87.8\pm 0.5$	$90.3\pm 0.3$
australian_scale2	$92.9\pm 0.2$	$89.6\pm 0.6$	$85.9\pm 0.7$	$92.4\pm 0.3$	$93.3\pm 0.2$
breast-cancer_scale	$98.9\pm 0.1$	$93.5\pm 1.8$	$54.8\pm 1.1$	$95.6\pm 0.4$	$99.5\pm 0.0$
breast-cancer_scale2	$14.4\pm 2.0$	$19.4\pm 3.8$	$91.5\pm 0.3$	$97.2\pm 0.2$	$99.0\pm 0.1$
diabetes_scale	$64.0\pm 1.3$	$63.8\pm 1.4$	$72.6\pm 0.5$	$74.7\pm 0.5$	$75.9\pm 0.5$
skin_nonskin	-	-	-	$94.8\pm 0.1$	$99.5\pm 0.1$
skin_nonskin2	-	-	-	$94.6\pm 0.0$	$99.8\pm 0.0$
epsilon_normalized	-	-	-	$69.3\pm 0.1$	$59.7\pm 1.8$
HIGGS	-	-	-	$65.3\pm 0.1$	$55.4\pm 0.2$

Equations

$$R_L(f)\text{mathchar}58=$$

$$=$$

$$R_L(f)={\int }_{\mathcal{X}}L(-1,f(x))d{P}_X(x)+{\pi }_{+}{\int }_{\mathcal{X}}L(1,f(x))-L(-1,f(x))d{P}_{X\mid Y=1}(x).$$

$$\mathrm{IPM}(P,Q;\mathcal{F})\text{mathchar}58=\underset{f\in \mathcal{F}}{\sup }\left|{\int }_{\mathcal{X}}f(x)dP(x)-{\int }_{\mathcal{X}}f(x)dQ(x)\right|.$$

$$\mathrm{WIPM}(P,Q;\tilde{w},\mathcal{F})\text{mathchar}58=\underset{f\in \mathcal{F}}{\sup }\left|{\int }_{\mathcal{X}}f(x)dP(x)-{\int }_{\mathcal{X}}\tilde{w}(x)f(x)dQ(x)\right|.$$

$$\mathrm{WIPM}(P,Q;w,\mathcal{F})\text{mathchar}58=\underset{f\in \mathcal{F}}{\sup }\left|{\int }_{\mathcal{X}}f(x)dP(x)-w{\int }_{\mathcal{X}}f(x)dQ(x)\right|.$$

$$\underset{f\in \mathcal{F}}{\inf }{R}_{{\ell }_{\mathrm{h}}}(f)=1-\mathrm{WIPM}(P_X,P_{X\mid Y=1};2{\pi }_{+},\mathcal{F}).$$

$$\mathrm{WIPM}(P_X,P_{X\mid Y=1};2{\pi }_{+},\mathcal{F})={\int }_{\mathcal{X}}{g}_{\mathcal{F}}(x)d{P}_X(x)-2{\pi }_{+}{\int }_{\mathcal{X}}{g}_{\mathcal{F}}(x)d{P}_{X\mid Y=1}(x),$$

$$\mathrm{WIPM}(P_{X,n_{\mathrm{u}}},P_{X\mid Y=1,n_{\mathrm{p}}};w,\mathcal{F})=\underset{f\in \mathcal{F}}{\sup }\left|\frac{1}{n_{\mathrm{u}}}\sum _{i=1}^{n_{\mathrm{u}}}f(x_i^{\mathrm{u}})-\frac{w}{n_{\mathrm{p}}}\sum _{i=1}^{n_{\mathrm{p}}}f(x_i^{\mathrm{p}})\right|,$$

$$\mathrm{WIPM}(P_{X,n_{\mathrm{u}}},P_{X\mid Y=1,n_{\mathrm{p}}};w,\mathcal{F})=\frac{1}{n_{\mathrm{u}}}\sum _{i=1}^{n_{\mathrm{u}}}{\hat{g}}_{\mathcal{F}}(x_i^{\mathrm{u}})-\frac{w}{n_{\mathrm{p}}}\sum _{i=1}^{n_{\mathrm{p}}}{\hat{g}}_{\mathcal{F}}(x_i^{\mathrm{p}}).$$

$$R_{{\ell }_{\mathrm{h}}}({\hat{f}}_{\mathcal{F}})-\underset{f\in \mathcal{F}}{\inf }{R}_{{\ell }_{\mathrm{h}}}(f)\le$$

$$+C_{\tau ,{\rho }^2}^{(1)}{\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(1)}(2{\pi }_{+})+C_{\tau ,\nu ,\alpha }^{(2)}{\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(2)}(2{\pi }_{+}),$$

$$R_{{\ell }_{\mathrm{h}}}(\hat{f})-\underset{f\in \mathcal{F}}{\inf }{R}_{{\ell }_{\mathrm{h}}}(f)\le$$

$$+{\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(1)}(2{\pi }_{+})(1+\nu )\sqrt{2\tau }+\Delta .$$

$$\frac{1+\nu }{2}-\frac{5\sqrt{2\tau }{\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(2)}(2{\pi }_{+})\nu }{6{\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(1)}(2{\pi }_{+})}\ge \rho .$$

$$R_{{\ell }_{01}}({\hat{f}}_{\mathcal{F}})-\underset{f\in \mathcal{U}}{\inf }{R}_{{\ell }_{01}}(f)\le$$

$$+C_{\alpha }({\mathbb{E}}_{P_X^{n_{\mathrm{u}}}}({\mathfrak{R}}_{{\mathcal{X}}_{\mathrm{u}}}(\mathcal{F}))+2{\pi }_{+}{\mathbb{E}}_{P_{X\mid Y=1}^{n_{\mathrm{p}}}}({\mathfrak{R}}_{{\mathcal{X}}_{\mathrm{p}}}(\mathcal{F})))$$

$$+C_{\tau ,{\rho }^2}^{(1)}{\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(1)}(2{\pi }_{+})+C_{\tau ,\nu ,\alpha }^{(2)}{\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(2)}(2{\pi }_{+}).$$

$${\mathrm{WMMD}}_k(P_{X,n_{\mathrm{u}}},P_{X\mid Y=1,n_{\mathrm{p}}};w,r)=\frac{1}{n_{\mathrm{u}}}\sum _{i=1}^{n_{\mathrm{u}}}{\hat{g}}_{{\mathcal{H}}_{k,r}}(x_i^{\mathrm{u}})-\frac{w}{n_{\mathrm{p}}}\sum _{i=1}^{n_{\mathrm{p}}}{\hat{g}}_{{\mathcal{H}}_{k,r}}(x_i^{\mathrm{p}}).$$

$$\mathrm{sign}({\hat{f}}_{{\mathcal{H}}_{k,r}}(z))=\{\begin{array}{ll}+1& \mathrm{if}(2{\pi }_{+}{)}^{-1}<{\hat{\lambda }}_{n_{\mathrm{p}},n_{\mathrm{u}}}(z),\\ -1& otherwise,\end{array}$$

$${\hat{\lambda }}_{n_{\mathrm{p}},n_{\mathrm{u}}}(z)=\frac{n_{\mathrm{p}}^{-1}{\sum }_{i=1}^{n_{\mathrm{p}}}k(z,x_i^{\mathrm{p}})}{n_{\mathrm{u}}^{-1}{\sum }_{i=1}^{n_{\mathrm{u}}}k(z,x_i^{\mathrm{u}})}.$$

$$R_{{\ell }_{\mathrm{h}}}({\hat{f}}_{{\mathcal{H}}_{k,r_1}})-\underset{f\in {\mathcal{H}}_{k,r_1}}{\inf }{R}_{{\ell }_{\mathrm{h}}}(f)\le (C_{\alpha }+C_{\tau ,{\rho }^2}^{(1)}){\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(1)}(2{\pi }_{+})+C_{\tau ,\nu ,\alpha }^{(2)}{\chi }_{n_{\mathrm{p}},n_{\mathrm{u}}}^{(2)}(2{\pi }_{+}).$$

$$\underset{f\in {\mathcal{H}}_{k,r_1}}{\inf }{R}_{{\ell }_{\mathrm{h}}}(f)-\underset{f\in \beta \mathcal{M}}{\inf }{R}_{{\ell }_{\mathrm{h}}}(f)\le \beta \underset{g\in {\mathcal{H}}_{k,r_1/\beta }}{\inf }{\Vert g-f_1^{\ast }\Vert }_{L_2(P_X)},$$

$$R_{{\ell }_{01}}({\hat{f}}_{{\mathcal{H}}_{k,1}})-\underset{f\in \mathcal{U}}{\inf }{R}_{{\ell }_{01}}(f)\le C_{(\mathrm{p},\mathrm{u})}(n_{\mathrm{p}}\wedge n_{\mathrm{u}}{)}^{-\frac{{\alpha }_{\text{H}}(1+q)}{2{\alpha }_{\text{H}}+d}},$$

$$X|Y=y,U$$

$$X\mid Y=y\sim N\left(y\frac{1_2}{\sqrt{2}},I_2\right),Y\sim 2\times \mathrm{Bern}({\pi }_{+})-1,$$

$$R_{{\ell }_{\mathrm{h}}}(f)$$

$$=1+{\int }_{\mathcal{X}}f(x)d{P}_X(x)-2{\pi }_{+}{\int }_{\mathcal{X}}f(x)d{P}_{X\mid Y=1}(x).$$

$$\underset{f\in \mathcal{F}}{\inf }{R}_{{\ell }_{\mathrm{h}}}(f)$$

$$=1-\underset{f\in \mathcal{F}}{\sup }\left\{-{\int }_{\mathcal{X}}f(x)d{P}_X(x)+2{\pi }_{+}{\int }_{\mathcal{X}}f(x)d{P}_{X\mid Y=1}(x)\right\}$$

$$\stackrel{(\ast )}{=}1-\underset{f\in \mathcal{F}}{\sup }\left\{{\int }_{\mathcal{X}}f(x)d{P}_X(x)-2{\pi }_{+}{\int }_{\mathcal{X}}f(x)d{P}_{X\mid Y=1}(x)\right\}$$

$$=1-\mathrm{WIPM}(P_X,P_{X\mid Y=1};2{\pi }_{+},\mathcal{F}).$$

$$\underset{f\in \mathcal{F}}{\inf }{R}_{{\ell }_{\mathrm{h}}}(f)$$

$$=1-\left\{{\int }_{\mathcal{X}}{g}_{\mathcal{F}}(x)d{P}_X(x)-2{\pi }_{+}{\int }_{\mathcal{X}}{g}_{\mathcal{F}}(x)d{P}_{X\mid Y=1}(x)\right\}$$

$$=1-\left\{R_{{\ell }_{\mathrm{h}}}(g_{\mathcal{F}})-1\right\}$$

$$=2-R_{{\ell }_{\mathrm{h}}}(g_{\mathcal{F}})=R_{{\ell }_{\mathrm{h}}}(-g_{\mathcal{F}}).$$

Full text

1 Introduction

Supervised binary classification assumes that all the training data are labeled as either being positive or negative. However, in many practical scenarios, collecting a large number of labeled samples from the two categories is often costly, difficult, or not even possible. In contrast, unlabeled data are relatively cheap and abundant. As a consequence, semi-supervised learning is used for partially labeled data (Chapelle et al., 2006). In this paper, as a special case of semi-supervised learning, we consider Positive-Unlabeld (PU) learning, the problem of building a binary classifier from only positive and unlabeled samples (Denis et al., 2005; Li and Liu, 2005). PU learning provides a powerful framework when negative labels are impossible or very expensive to obtain, and thus has frequently appeared in many real-world applications. Examples include document classification (Elkan and Noto, 2008; Xiao et al., 2011), image classification (Zuluaga et al., 2011; Gong et al., 2018), gene identification (Yang et al., 2012, 2014), and novelty detection (Blanchard et al., 2010; Zhang et al., 2017).

Several PU learning algorithms have been developed over the last two decades. Liu et al. (2002) and Li and Liu (2003) considered a two-step learning scheme: in Step 1, assigning negative labels to some unlabeled observations believed to be negative, and in Step 2, learning a binary classifier with existing positive samples and the negatively labeled samples from Step 1. Liu et al. (2003) pointed out that the two-step learning scheme is based on heuristics, and suggested fitting a biased support vector machine by regarding all the unlabeled observations as being negative.

Scott and Blanchard (2009) and Blanchard et al. (2010) suggested a modification of supervised Neyman-Pearson classification, whose goal is to find a classifier minimizing the false positive rate keeping the false negative rate low. To circumvent the problem of lack of negative samples, they tried to build a classifier minimizing the marginal probability of being classified as positive while keeping the false negative rate low. Solving the empirical version of this constrained optimization problem is challenging, but the authors did not present an explicit algorithm.

Recently, many PU learning algorithms based on the empirical risk minimization principle have been studied. du Plessis et al. (2014) proposed the use of the ramp loss and provided an algorithm that requires solving a non-convex optimization problem. du Plessis et al. (2015) formulated a convex optimization problem by using the logistic loss or double hinge loss. However, all the aforementioned approaches involve solving a non-linear programming problem. This causes massive computational burdens for calculating the large Gram matrix when the sample size is large. Kiryo et al. (2017) suggested a stochastic algorithm for large-scale datasets with a non-negative risk estimator. However, to execute the algorithm, several hyperparameters are required, and choosing the optimal hyperparameter may demand substantial trials of running the algorithm (Oh et al., 2018), causing heavy computation costs.

In supervised binary classification, Sriperumbudur et al. (2012) proposed a computationally efficient algorithm building a closed-form binary discriminant function. The authors showed that their function estimator obtained by evaluating the negative of the empirical integral probability metric (IPM) is the minimizer of the empirical risk using the specific loss defined in Section 3.1. They further showed that a closed form can be derived as the result of restricting a hypothesis space to a closed unit ball in reproducing kernel Hilbert space (RKHS).

In this paper, capitalizing on the properties shown in the supervised learning method by Sriperumbudur et al. (2012), we extend it to PU learning settings. In addition, we derive new theoretical results on excess risk bounds. We first define a weighted version of IPM between two probability measures and call it the weighted integral probability metric (WIPM). We show that computing the negative of WIPM between the unlabeled data distribution and the positive data distribution is equivalent to minimizing the hinge risk. Based on this finding, we propose a binary discriminant function estimator that computes the negative of the empirical WIPM, and then derive associated upper bounds of the estimation error and the excess risk. Under a mild condition, our obtained upper bounds are shown to be sharper than the existing ones because of using Talagrand’s inequality over McDiarmid’s inequality (Kiryo et al., 2017). Moreover, we pay special attention to the case where the hypothesis space is a closed ball in RKHS and propose a closed-form classifier. We show that the associated excess risk bound has an explicit form that converges to zero as the sample sizes increase. To the best of our knowledge, this is the first result to explicitly show the excess risk bound in PU learning.

As a summary, our main contributions are:

• •

We formally define WIPM and establish a link with the infimum of the hinge risk (Theorem 3.1). We derive an estimation error bound and show that it is sharper than existing results (Theorem 3.2 and Proposition 1).

• •

The proposed algorithm produces a closed-form classifier when the underlying hypothesis space is a closed ball in RKHS (Proposition 2). Furthermore, we obtain a novel excess risk bound that converges to zero as sample sizes increase (Theorem 4.1).

• •

Numerical experiments using both synthetic and real datasets show that our method is comparable to or better than existing PU learning algorithms in terms of accuracy, scalability, and robustness in the case of unknown class-priors.

2 Preliminaries

In this section, we describe the $L$-risk for binary classification and present its PU representation. We briefly review several PU learning algorithms based on the $L$-risk minimization principle. We first introduce problem settings and notations.

2.1 Problem settings of PU learning

Let $X$ and $Y$ be random variables for input data and class labels, respectively, whose range is the product space $\mathcal{X}\times\{\pm 1\}\subseteq\mathbb{R}^{d}\times\{\pm 1\}$. The $d$ is a positive integer. We denote the joint distribution of $(X,Y)$ by $P_{X,Y}$ and the marginal distribution of $X$ by $P_{X}$. The distributions of positive and negative data are defined by conditional distributions, ${P}_{X\mid Y=1}$ and ${P}_{X\mid Y=-1}$, respectively. Let $\pi_{+}\mathrel{\mathop{\mathchar 58\relax}}=P_{X,Y}(Y=1)$ be the marginal probability of being positive and set $\pi_{-}=1-\pi_{+}$. We follow the two samples of data scheme (Ward et al., 2009; Niu et al., 2016). That is, let $\mathcal{X}_{\rm p}=\{x_{i}^{\rm p}\}_{i=1}^{n_{\rm p}}$ and $\mathcal{X}_{\rm u}=\{x_{i}^{\rm u}\}_{i=1}^{n_{\rm u}}$ be observed sets of independently identically distributed samples from the positive data distribution ${P}_{X\mid Y=1}$ and the marginal distribution ${P}_{X}$, respectively. Here, the $n_{\rm p}$ and $n_{\rm u}$ are the number of positive and unlabeled data points, respectively. Note that the unlabeled data distribution is the marginal distribution.

Let $\mathcal{U}$ be a class of real-valued measurable functions defined on $\mathcal{X}$. A function $f\in\mathcal{U}$, often called a hypothesis, can be understood as a binary discriminant function and we classify an input $x$ with the sign of a discriminant function, ${\rm sign}(f(x))$. Define $\mathcal{M}=\{f\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}\to\mathbb{R}\mid\mathinner{\!\left\lVert f\right\rVert}_{\infty}\leq 1\}\subseteq\mathcal{U}$, where $\mathinner{\!\left\lVert f\right\rVert}_{\infty}=\sup_{x\in\mathcal{X}}|f(x)|$ is the supremum norm. We restrict our attention to a class $\mathcal{F}\subseteq\mathcal{M}$ and call $\mathcal{F}$ a hypothesis space. Throughout this paper, we assume that the hypothesis space is symmetric, i.e., $f\in\mathcal{F}$ implies $-f\in\mathcal{F}$. In PU learning, the main goal is to construct a classifier ${\rm sign}(f(x))$ only from the positive dataset $\mathcal{X}_{\rm p}$ and the unlabeled dataset $\mathcal{X}_{\rm u}$ with $f\in\mathcal{F}$.

In this paper, the quantity $\pi_{+}$, often called the class-prior, is assumed to be known as in the literature (Kiryo et al., 2017; Kato et al., 2019) to focus on theoretical and practical benefits of our proposed algorithm. We examine the performance when $\pi_{+}$ is unknown in Experiment 3 of Section 6.1 and in Section 6.2.

2.2 $L$-risk minimization in PU learning

In supervised binary classification, the $L$-risk is defined by

[TABLE]

for a loss function $L\mathrel{\mathop{\mathchar 58\relax}}\{\pm 1\}\times\mathbb{R}\to\mathbb{R}$ (Steinwart and Christmann, 2008, Section 2.1). We denote the margin-based loss function by $\ell(yt)\mathrel{\mathop{\mathchar 58\relax}}=L(y,t)$ if a loss function $L(y,t)$ can be represented as a function of margin $yt$, the product of a label $y$ and a score $t$ for all possible $y\in\{\pm 1\}$ and $t\in\mathbb{R}$.

Under the PU learning framework, however, the right-hand side of Equation (1) cannot be directly estimated due to lack of negatively labeled observations. To circumvent this problem, many studies in the field of PU learning exploited the relationship $P_{X}=\pi_{+}P_{X\mid Y=1}+\pi_{-}P_{X\mid Y=-1}$ and replaced $P_{X\mid Y=-1}$ in Equation (1) with $(P_{X}-\pi_{+}P_{X\mid Y=1})/\pi_{-}$ (du Plessis et al., 2014; Sakai et al., 2017). That is, the $L$-risk can be alternatively expressed as:

[TABLE]

Now the right-hand side of Equation (2) can be empirically estimated by the positive dataset $\mathcal{X}_{\rm p}$ and the unlabeled dataset $\mathcal{X}_{\rm u}$. However, the $L$-risk $R_{L}(f)$ is not convex with respect to $f$ in general, and minimizing an empirical estimator for $R_{L}(f)$ is often formulated as a complicated non-convex optimization problem.

There have been several approaches to resolving the computational difficulty by modifying loss functions. du Plessis et al. (2014) proposed to use non-convex loss functions satisfying the symmetric condition, $L(1,f(x))+L(-1,f(x))=1$. They proposed to optimize the empirical risk based on the ramp loss $\ell_{\mathrm{ramp}}(yt)=0.5\times\max(0,$ $\min(2,1-yt))$ via the concave-convex procedure (Collobert et al., 2006). du Plessis et al. (2015) converted the problem to convex optimization through the linear-odd condition, $L(1,f(x))-L(-1,f(x))=-f(x)$. They showed that the logistic loss $\ell_{\mathrm{log}}(yt)=\log(1+\exp(-yt))$ and the double hinge loss $\ell_{\mathrm{dh}}(yt)=\max(0,\max(-yt,(1-yt)/2))$ satisfy the linear-odd condition. However, all the aforementioned methods utilized a weighted sum of $n_{\rm p}+n_{\rm u}$ predefined basis functions as a binary discriminant function, which triggered calculating the $(n_{\rm p}+n_{\rm u})$ $\times$ $(n_{\rm p}+n_{\rm u})$ Gram matrix. Hence, executing algorithms is not scalable and can be intractable when $n_{\rm p}$ and $n_{\rm u}$ are large (Sansone et al., 2018). Our first goal in this paper is to overcome this computational problem by providing a computationally efficient method.

3 Weighted integral probability metric and $L$-risk

In this section, we formally define WIPM, a key tool for constructing the proposed algorithm, and build a link with the $L$-risk in Theorem 3.1 below. Based on the link, we propose a new binary discriminant function estimator and present its theoretical properties in Theorem 3.2. We first introduce the earlier work by Sriperumbudur et al. (2012) that provided a closed-form classifier in supervised binary classification.

3.1 Relation between IPM and $L$-risk in supervised binary classification

Müller (1997) introduced an IPM for any two probability measures ${P}$ and ${Q}$ defined on $\mathcal{X}$ and a class $\mathcal{F}$ of bounded measurable functions, given by

[TABLE]

IPM has been studied as either a metric between two probability measures (Sriperumbudur et al., 2010a; Arjovsky et al., 2017; Tolstikhin et al., 2018) or a hypothesis testing tool (Gretton et al., 2012).

Under the supervised binary classification setting, Sriperumbudur et al. (2012) showed that calculating IPM between $P_{X\mid Y=1}$ and $P_{X\mid Y=-1}$ is negatively related to minimizing the risk with a loss function, i.e., ${\rm IPM}(P_{X\mid Y=1},P_{X\mid Y=-1};\mathcal{F})=-\inf_{f\in\mathcal{F}}R_{L_{\mathrm{c}}}(f)$, where $L_{\mathrm{c}}(1,t)=-t/\pi_{+}$ and $L_{\mathrm{c}}(-1,t)=t/\pi_{-}$ for all $t\in\mathbb{R}$. They further showed that a discriminant function minimizing the $L_{\mathrm{c}}$-risk can be obtained analytically when $\mathcal{F}$ is a closed unit ball in RKHS. This result cannot be directly extended to PU learning due to absence of negatively labeled observations. In the next subsection, we define a generalized version of IPM and extend the previous results for supervised binary classification to PU learning.

3.2 Extension to WIPM and $L$-risk in PU learning

Let $\mathcal{F}$ be a given class of bounded measurable functions and let $\tilde{w}\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}\to\mathbb{R}$ be a weight function such that $\mathinner{\!\left\lVert\tilde{w}\right\rVert}_{\infty}<\infty$. We define WIPM111Although WIPM is not a metric in general, we keep saying the name WIPM to emphasize that it is a weighted version of IPM. between two probability measures ${P}$ and ${Q}$ with a function class $\mathcal{F}$ and a weight function $\tilde{w}$ by

[TABLE]

Note that WIPM reduces to IPM if $\tilde{w}(x)=1$ for all $x\in\mathcal{X}$. Other special cases of Equation (3) have been discussed in many applications. In the covariate shift problem, Huang et al. (2007) and Gretton et al. (2009) proposed to minimize WIPM with respect to $\tilde{w}$ when $\mathcal{F}$ is the unit ball in RKHS and $P,Q$ are empirical distributions of test and training data, respectively. In unsupervised domain adaptation, Yan et al. (2017) regarded $P,Q$ as empirical distributions of target and source data, respectively, where in this case, $\tilde{w}$ is a ratio of two class-prior distributions.

We pay special attention to the case where $\tilde{w}(x)$ is constant, $w\in\mathbb{R}$, for every input value and denote WIPM by ${\rm WIPM}({P},{Q};w,\mathcal{F})$,

[TABLE]

In the following theorem, we establish a link between ${\rm WIPM}(P_{X},P_{X\mid Y=1};2\pi_{+},\mathcal{F})$ and the infimum of the $\ell_{\mathrm{h}}$-risk over $\mathcal{F}$ for the hinge loss $\ell_{\mathrm{h}}(yt)=\max(0,1-yt)$. A proof is in Appendix A.

Theorem 3.1 (Relationship between $\ell_{\mathrm{h}}$-risk and WIPM).

Let $\mathcal{F}$ be a symmetric hypothesis space in $\mathcal{M}$ and $\ell_{\mathrm{h}}(yt)=\max(0,1-yt)$ be the hinge loss. Then, we have

[TABLE]

Moreover, if $g_{\mathcal{F}}$ satisfies

[TABLE]

then $\inf_{f\in\mathcal{F}}R_{\ell_{\mathrm{h}}}(f)$ $=R_{\ell_{\mathrm{h}}}(-g_{\mathcal{F}})$.

Theorem 3.1 shows that the infimum of the $\ell_{\mathrm{h}}$-risk over a hypothesis space $\mathcal{F}$ equals the negative WIPM between the unlabeled data distribution $P_{X}$ and the positive data distribution $P_{X\mid Y=1}$ with the same hypothesis space $\mathcal{F}$ and the weight $2\pi_{+}$ up to addition by constant. Furthermore, by negating the WIPM optimizer $g_{\mathcal{F}}$, we obtain the minimizer of the $\ell_{\mathrm{h}}$-risk over the hypothesis space $\mathcal{F}$. Here, we define a WIPM optimizer $g_{\mathcal{F}}$ as a function that attains the supremum, i.e., ${\rm WIPM}(P_{X},P_{X\mid Y=1};2\pi_{+},\mathcal{F})=\int_{\mathcal{X}}g_{\mathcal{F}}(x)$ $dP_{X}(x)$ $-2\pi_{+}\int_{\mathcal{X}}g_{\mathcal{F}}(x)dP_{X\mid Y=1}(x)$ and we set $f_{\mathcal{F}}=-g_{\mathcal{F}}$ for later notational convenience. Sriperumbudur et al. (2012) derived a similar result to Theorem 3.1 by showing ${\rm IPM}(P_{X\mid Y=1},P_{X\mid Y=-1};\mathcal{F})$ $=$ $-\inf_{f\in\mathcal{F}}R_{L_{\mathrm{c}}}(f)$ in supervised binary classification. However, as we mentioned in Section 3.1, their method is only applicable to supervised binary classification settings.

3.3 Theoretical properties of empirical WIPM optimizer

We denote the empirical distributions of ${P}_{X\mid Y=1}$ and ${P}_{X}$ by ${P}_{X\mid Y=1,n_{\rm p}}$ and ${P}_{X,n_{\rm u}}$, respectively. Let ${P}_{X\mid Y=1,n_{\rm p}}=n_{\rm p}^{-1}\sum_{i=1}^{n_{\rm p}}\delta_{x_{i}^{\rm p}}$ and ${P}_{X,n_{\rm u}}=n_{\rm u}^{-1}\sum_{i=1}^{n_{\rm u}}\delta_{x_{i}^{\rm u}}$, where $\delta(\cdot)$ defined on $\mathcal{X}$ is the Dirac delta function and $\delta_{x}(\cdot)\mathrel{\mathop{\mathchar 58\relax}}=\delta(\cdot-x)$ for $x\in\mathcal{X}$. The empirical Rademacher complexity of $\mathcal{F}$ given a set $S=\{z_{1},\dots,z_{m}\}$ is defined by $\mathfrak{R}_{S}(\mathcal{F})\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}_{\sigma}\left(\frac{1}{m}\sup_{f\in\mathcal{F}}\left|\sum_{i=1}^{m}\sigma_{i}f(z_{i})\right|\right)$. Here, $\{\sigma_{i}\}_{i=1}^{m}$ is a set of independent Rademacher random variables taking $1$ or $-1$ with probability $0.5$ each and $\mathbb{E}_{\sigma}(\cdot)$ is the expectation operator over the Rademacher random variables (Bartlett and Mendelson, 2002). Denote a maximum by $a\vee b\mathrel{\mathop{\mathchar 58\relax}}=\max(a,b)$, a minimum by $a\wedge b\mathrel{\mathop{\mathchar 58\relax}}=\min(a,b)$. For a probability measure ${Q}$ defined on $\mathcal{X}$, denote the expectation of a discriminant function $f$ by $\mathbb{E}_{Q}(f)\mathrel{\mathop{\mathchar 58\relax}}=\int_{\mathcal{X}}f(x)d{Q}(x)$ and the variance by $\mathrm{Var}_{{Q}}(f)\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}_{Q}(f^{2})-(\mathbb{E}_{Q}(f))^{2}$.

The empirical estimator for ${\rm WIPM}({P}_{X},{P}_{X\mid Y=1};w,\mathcal{F})$ is given by plugging in the empirical distributions,

[TABLE]

and we define an empirical WIPM optimizer $\hat{g}_{\mathcal{F}}\in\mathcal{F}$ that satisfies the following equation,

[TABLE]

We set $\hat{f}_{\mathcal{F}}=-\hat{g}_{\mathcal{F}}$ for notational convenience as in Section 3.2.

We analyze the estimation error $R_{\ell_{\mathrm{h}}}(\hat{f}_{\mathcal{F}})-\inf_{f\in\mathcal{F}}R_{\ell_{\mathrm{h}}}(f)$ in the following theorem. A proof is provided in Appendix B.2. To begin, let $\chi_{n_{\rm p},n_{\rm u}}^{(1)}(w)=w/\sqrt{n_{\rm p}}+1/\sqrt{n_{\rm u}}$ and $\chi_{n_{\rm p},n_{\rm u}}^{(2)}(w)=2(w/n_{\rm p}+1/n_{\rm u})$.

Theorem 3.2 (Estimation error bound for general function space).

Let $\hat{g}_{\mathcal{F}}$ be an empirical WIPM optimizer defined in Equation (4) and set $\hat{f}_{\mathcal{F}}=-\hat{g}_{\mathcal{F}}$. Let $\mathcal{F}$ be a symmetric hypothesis space such that $\mathinner{\!\left\lVert f\right\rVert}_{\infty}\leq\nu\leq 1$, $\mathrm{Var}_{{P}_{X\mid Y=1}}(f)\leq\sigma_{X\mid Y=1}^{2}$, and $\mathrm{Var}_{{P}_{X}}(f)\leq\sigma_{X}^{2}$. Denote $\rho^{2}=\sigma_{X\mid Y=1}^{2}\vee\sigma_{X}^{2}$. Then, for all $\alpha,\tau>0$, the following holds with probability at least $1-e^{-\tau}$,

[TABLE]

where $C_{\alpha}=4(1+\alpha)$, $C_{\tau,\rho^{2}}^{(1)}=2\sqrt{2\tau\rho^{2}}$, $C_{\tau,\nu,\alpha}^{(2)}=2\tau\nu\left(\frac{2}{3}+\frac{1}{\alpha}\right)$.

Due to Talagrand’s inequality, Theorem 3.2 provides a sharper bound than the existing result based on McDiarmid’s inequality. Specifically, Kiryo et al. (2017, Theorem 4) utilized McDiarmid’s inequality and showed that for $\tau>0$ and some $\Delta>0$ the following holds with probability at least $1-e^{-\tau}$,

[TABLE]

The following proposition shows that the proposed upper bound (5) is sharper than the upper bound (6) under a certain condition. A proof is provided in Appendix B.3.

Proposition 1.

With the notations defined in Theorem 3.2, suppose that the following holds,

[TABLE]

Then, the proposed upper bound (5) is sharper than the previous result (6) proposed by Kiryo et al. (2017).

It is noteworthy that the second term in the left-hand side of (7) converges to zero as $n_{\rm p}$ and $n_{\rm u}$ increase because $\chi_{n_{\rm p},n_{\rm u}}^{(1)}(2\pi_{+})$ $=$ $O_{{P}_{X\mid Y=1},{P}_{X}}((n_{\rm p}\wedge n_{\rm u})^{-1/2})$ and $\chi_{n_{\rm p},n_{\rm u}}^{(2)}(2\pi_{+})$ $=$ $O_{{P}_{X\mid Y=1},{P}_{X}}((n_{\rm p}\wedge n_{\rm u})^{-1})$. Due to $(1+\nu)/2\geq\nu\geq\rho$, the condition (7) is quite reasonable if the upper bounds of the variances, $\sigma_{X}^{2}$ and $\sigma_{X\mid Y=1}^{2}$, are sufficiently small.

In binary classification, one ultimate goal is to find a classifier minimizing the misclassification error, or equivalently, minimizing the excess risk. Bartlett et al. (2006) showed that there is an invertible function $\psi\mathrel{\mathop{\mathchar 58\relax}}[-1,1]\to[0,\infty)$ such that the excess risk $R_{\ell_{01}}(\hat{f}_{\mathcal{F}})-\inf_{f\in\mathcal{U}}R_{\ell_{01}}(f)$ is bounded above by $\psi^{-1}(R_{\ell}(\hat{f}_{\mathcal{F}})-\inf_{f\in\mathcal{U}}R_{\ell}(f))$ if the margin-based loss $\ell$ is classification-calibrated. In particular, Zhang (2004) showed that the excess risk is bounded above by the excess $\ell_{\mathrm{h}}$-risk, i.e., $R_{\ell_{01}}(\hat{f}_{\mathcal{F}})-\inf_{f\in\mathcal{U}}R_{\ell_{01}}(f)$ $\leq$ $R_{\ell_{\mathrm{h}}}(\hat{f}_{\mathcal{F}})-\inf_{f\in\mathcal{U}}R_{\ell_{\mathrm{h}}}(f)$. This implies that an excess risk bound can be obtained by analyzing the excess $\ell_{\mathrm{h}}$-risk bound with Theorem 3.2. The following corollary provides the excess risk bound.

Corollary 1 (Excess risk bound for general function space).

With the notations defined in Theorem 3.2, for all $\alpha,\tau>0$, the following holds with probability at least $1-e^{-\tau}$,

[TABLE]

4 WIPM optimizer with reproducing kernel Hilbert space

In this section, we provide a computationally efficient PU learning algorithm which builds an analytic classifier when a hypothesis space is a closed ball in RKHS. In addition, unlike the excess risk bound in Corollary 1, we explicitly derive the bound that converges to zero when the sample sizes $n_{\rm p}$ and $n_{\rm u}$ increase.

4.1 An analytic classifier via WMMD optimizer

To this end, we assume that $\mathcal{X}\subseteq[0,1]^{d}$ is compact. Let $k\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ be a reproducing kernel defined on $\mathcal{X}$ and $\mathcal{H}_{k}$ be the associated RKHS with the inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}_{k}}\mathrel{\mathop{\mathchar 58\relax}}\mathcal{H}_{k}\times\mathcal{H}_{k}\to\mathbb{R}$. We denote the induced norm by $\mathinner{\!\left\lVert\cdot\right\rVert}_{\mathcal{H}_{k}}$. Denote a closed ball in RKHS $\mathcal{H}_{k}$ with a radius $r>0$, by $\mathcal{H}_{k,r}=\{f\mathrel{\mathop{\mathchar 58\relax}}\mathinner{\!\left\lVert f\right\rVert}_{\mathcal{H}_{k}}\leq r\}$. We define the weighted maximum mean discrepancy (WMMD) between two probability measures ${P}$ and ${Q}$ with a weight $w$ and a closed ball $\mathcal{H}_{k,r}$ by ${\rm WMMD}_{k}({P},{Q};w,r)$ $\mathrel{\mathop{\mathchar 58\relax}}=$ ${\rm WIPM}({P},{Q};w,\mathcal{H}_{k,r})$. The name of WMMD comes from the maximum mean discrepancy (MMD), a popular example of the IPM whose function space is the unit ball ${\mathcal{H}}_{k,1}$, i.e., ${\rm MMD}_{k}({P},{Q})\mathrel{\mathop{\mathchar 58\relax}}={\rm IPM}({P},{Q};\mathcal{H}_{k,1})$ (Sriperumbudur et al., 2010a, b). As defined in Equation (4), let $\hat{g}_{\mathcal{H}_{k,r}}\in\mathcal{H}_{k,r}$ be the empirical WMMD optimizer such that

[TABLE]

In addition, we set $\hat{f}_{\mathcal{H}_{k,r}}=-\hat{g}_{\mathcal{H}_{k,r}}$, which leads the corresponding classification rule to ${\rm sign}(\hat{f}_{\mathcal{H}_{k,r}}(z))$. In the following proposition, we show that this classification rule has an analytic expression by exploiting the reproducing property $f(x)=\langle f,k(\cdot,x)\rangle_{\mathcal{H}_{k}}$ and the Cauchy-Schwarz inequality. A proof is provided in Appendix C.1.

Proposition 2.

Let $k\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ be a bounded reproducing kernel. Then, the classification rule has a closed-form expression given by

[TABLE]

where

[TABLE]

We call the classifier defined in Equation (8) the WMMD classifier and the score $\hat{\lambda}_{n_{\rm p},n_{\rm u}}(z)$ the WMMD score for $z$. One strength of the WMMD classifier is that the classification rule has a closed-form expression, resulting in computational efficiency. Furthermore, the WMMD score $\hat{\lambda}_{n_{\rm p},n_{\rm u}}$ is independent of the class-prior $\pi_{+}$, and thus we can obtain the score function without prior knowledge of the class-prior.

4.2 Explicit excess risk bound of WMMD classifier

Since the empirical WMMD optimizer $\hat{g}_{\mathcal{H}_{k,r}}$ is a special case of the empirical WIPM optimizer, we have an excess risk bound from the result of Corollary 1. However, without knowing convergence rates of the Rademacher complexities, $\mathbb{E}_{P_{X}^{n_{\rm u}}}(\mathfrak{R}_{\mathcal{X}_{\rm u}}(\mathcal{F}))$ and $\mathbb{E}_{P_{X\mid Y=1}^{n_{\rm p}}}(\mathfrak{R}_{\mathcal{X}_{\rm p}}(\mathcal{F}))$, and the approximation error, the consistency of the classifier remains unclear. In this subsection, we establish an explicit excess risk bound that vanishes. We first derive an explicit estimation error bound in the following lemma.

Lemma 1 (Explicit estimation error bound).

With the notations defined in Theorem 3.2, assume that a reproducing kernel $k$ defined on a compact space $\mathcal{X}$ is bounded. Let $r_{1}^{-1}=\sup_{x\in\mathcal{X}}\sqrt{k(x,x)}$. Then, we have $\mathcal{H}_{k,r_{1}}\subseteq\mathcal{M}$. Moreover, for all $\alpha,\tau>0$, the following holds with probability at least $1-e^{-\tau}$,

[TABLE]

While the bound in Theorem 3.2 is expressed in terms $\mathbb{E}_{P_{X}^{n_{\rm u}}}(\mathfrak{R}_{\mathcal{X}_{\rm u}}(\mathcal{F}))$ and $\mathbb{E}_{P_{X\mid Y=1}^{n_{\rm p}}}(\mathfrak{R}_{\mathcal{X}_{\rm p}}(\mathcal{F}))$, these are evaluated in terms of $n_{p}$ and $n_{u}$ in the upper bound in Lemma 1, giving an explicit estimation error bound with $O((n_{\rm p}\wedge n_{\rm u})^{-1/2})$ convergence rate. The key idea is to use reproducing property $f(x)=\langle f,k(\cdot,x)\rangle_{\mathcal{H}_{k}}$ and the Cauchy-Schwarz inequality to obtain an upper bound for the Rademacher complexity. Detailed proofs are given in Appendix C.2.

In the following lemma, we elaborate on the approximation error bound. To begin, for any $0<\beta\leq 1$, let $\beta\mathcal{M}\mathrel{\mathop{\mathchar 58\relax}}=\{\beta f\mathrel{\mathop{\mathchar 58\relax}}f\in\mathcal{M}\}$. Set $f_{1}^{*}(x)={\rm sign}(P(Y=1\mid X=x)-\frac{1}{2})$.

Lemma 2 (Approximation error bound over uniformly bounded hypothesis space).

With the notations defined in Lemma 1, we have

[TABLE]

for any $0<\beta\leq 1$.

When $\beta=1$, Lemma 2 implies that the approximation error $\inf_{f\in\mathcal{H}_{k,r_{1}}}R_{\ell_{\mathrm{h}}}(f)$ $-\inf_{f\in\mathcal{U}}R_{\ell_{\mathrm{h}}}(f)$ is bounded above by $\inf_{g\in\mathcal{H}_{k,{r_{1}}}}\mathinner{\!\left\lVert g-f_{1}^{*}\right\rVert}_{L_{2}(P_{X})}$ due to $\inf_{f\in\mathcal{U}}R_{\ell_{\mathrm{h}}}(f)$ $=\inf_{f\in\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)$ (Lin, 2002). Hence, a naive substitution to Corollary 1 will give a sub-optimal bound because $\inf_{g\in\mathcal{H}_{k,{r_{1}}}}\mathinner{\!\left\lVert g-f_{1}^{*}\right\rVert}_{L_{2}(P_{X})}$ is non-zero in general.

In the following theorem, we rigorously establish the explicit excess risk bound which vanishes as $n_{\rm p}$ and $n_{\rm u}$ increase.

Theorem 4.1.

Assume that the Gaussian kernel $k(x,y)=\exp$ $(-\frac{\mathinner{\!\left\lVert x-y\right\rVert}^{2}}{2h^{2}})$ is used. Under the assumptions (A1)-(A4) in Appendix C.3, for $h=(n_{\rm p}\wedge n_{\rm u})^{-\frac{1}{2\alpha_{\text{H}}+d}}$, we have the following holds with probability at least $1-1/n_{\rm p}-1/n_{\rm u}$:

[TABLE]

where $\alpha_{\text{H}}$ and $q$ are defined in Appendix C.3 and $C_{(\rm p,\rm u)}>0$ is some constant.

In supervised binary classification settings, a similar result is obtained by Audibert et al. (2007, Theorem 3.3) and the convergence rate is called a super-fast rate when $\alpha_{\text{H}}q>d$. However, $\alpha_{\text{H}}$ and $q$ cannot be simultaneously very large (Audibert et al., 2007).

Niu et al. (2016) provided the excess risk bound expressed as a function of $n_{\rm p},n_{\rm u}$. However, their bound included combined terms of the approximation error and the Rademacher complexity, as in Corollary 1. To the best of our knowledge, we are the first to explicitly derive the excess risk bound with convergence rate in terms of a function of $n_{\rm p},n_{\rm u}$ in PU learning.

5 Related work

Excess risk bound in noisy label literature: PU learning can be considered as a special case of classification with asymmetric label noise, and many studies in this literature have shown consistency results similar to Theorem 4.1 (Natarajan et al., 2013). Patrini et al. (2016) derived an explicit estimation error when $\mathcal{F}$ is a set of linear hypotheses and Blanchard et al. (2016) showed a consistency result of the excess risk bound when the hypothesis space is RKHS with universal kernels. While the two studies assumed the one sample of data scheme, the proposed bound is based on the two samples of data scheme. Therefore, our proposed excess risk bound is expressed in $n_{\rm p}$ and $n_{\rm u}$, giving a new consistency theory.

Closed-form classifier: Blanchard et al. (2010) suggested a score function similar to the WMMD score by using different bandwidth hyperparameters for the denominator and the numerator. However, with these differences, our method gains theoretical justification while their score function does not. du Plessis et al. (2015) derived a closed-form classifier based on the squared loss. They estimated $P(Y=1\mid X)-P(Y=-1\mid X)$ and showed the consistency of the estimation error bound in the two samples of data scheme. However, the classifier is not scalable because it requires to compute the inverse of a $(n_{\rm p}+n_{\rm u})$ $\times$ $(n_{\rm p}+n_{\rm u})$ matrix.

6 Numerical experiments

In this section, we empirically analyze the proposed algorithm to demonstrate its practical efficacy using synthetic and real datasets. Optimization procedures and the selection of hyperparameters are detailed in Appendix D. Pytorch implementation for the experiments is available at https://github.com/eraser347/WMMD_PU.

6.1 Synthetic data analysis

We first visualize the effect of increasing the sample sizes $n_{\rm p}$ and $n_{\rm u}$ on the discriminant ability of the proposed algorithm (Experiment 1). Then we compare performance with (i) the logistic loss $\ell_{\mathrm{log}}$, denoted by LOG, (ii) the double hinge loss $\ell_{\mathrm{dh}}$, denoted by DH, both proposed by du Plessis et al. (2015), (iii) the non-negative risk estimator method, denoted by NNPU, proposed by Kiryo et al. (2017), (iv) the threshold adjustment method, denoted by tADJ, proposed by Elkan and Noto (2008), and (v) the proposed algorithm, denoted by WMMD (Experiments 2, 3, and 4).

Experiment 1: In this case, we used the two_moons dataset whose underlying distributions are

[TABLE]

where $U$ refers to a uniform random variable ranges from 0 to 1 and $N(\mu,\Sigma)$ is the normal distribution with mean $\mu$ and covariance $\Sigma$. We used the ‘make_moons’ function in the Python module ‘sklearn.datasets’ (Pedregosa et al., 2011) to generate the datasets.

Figure 1 illustrates the decision boundaries of WMMD using the two_moons dataset. The first row displays the case where the unlabeled sample size is small, $n_{\rm u}=50$, and the second row displays the case where the unlabeled sample size is large, $n_{\rm u}=400$. The first and second columns display the case where the positive sample sizes are $n_{\rm p}=5$ and $n_{\rm p}=10$, respectively. The class-prior is fixed to $\pi_{+}=0.5$, and we assumed that the class-prior is known. We visualize the true mean function of the positive and negative data distributions with blue and red lines, respectively. The positive data are represented by blue diamond points, and the unlabeled data are represented by gray points. The decision boundaries of the WMMD classifier tend to correctly separate the two clusters as $n_{\rm p}$ and $n_{\rm u}$ increase.

In Experiments 2, 3, and 4, we evaluate: (i) the accuracy and area under the receiver operating characteristic curve (AUC) as $n_{\rm u}$ and $\pi_{+}$ change when the class-prior is known (Experiment 2) and unknown (Experiment 3); (ii) the elapsed training time (Experiment 4). In these experiments, we set up the underlying joint distribution as follows:

[TABLE]

where $\mathrm{Bern}(p)$ is the Bernoulli distribution with mean $p$, $\boldsymbol{1}_{2}=(1,1)^{T}$ is the 2 dimensional vector of all ones and $I_{2}$ is the identity matrix of size 2.

Experiment 2: In this experiment, we compare the accuracy and AUC of the five PU learning algorithms when the true class-prior $\pi_{+}$ is known. Figures 2(a) and 2(c) show the accuracy and AUC on various $n_{\rm u}$. The training sample size for the positive data is $n_{\rm p}=100$ and the class prior is $\pi_{+}=0.5$. The unlabeled sample size changes from 40 to 500 by 20. We repeat a random generation of training and test data 100 times. For comparison purposes, we add the $1-$Bayes risk for each unlabeled sample size. In terms of accuracy, the proposed WMMD tends to be closer to the $1-$Bayes risk as the $n_{\rm u}$ increases. Compared with other PU learning algorithms, WMMD achieves higher accuracy in every $n_{\rm u}$ and achieves comparable to or better AUC.

Figures 2(b) and 2(d) show a comparison of accuracy and AUC as $\pi_{+}$ changes. The training sample size for the positive and unlabeled data are $n_{\rm p}=100$ and $n_{\rm u}=400$, respectively. The class-prior $\pi_{+}$ changes from $0.05$ to $0.95$ by $0.05$. The test sample size is $10^{3}$. Training and test data are repeatedly generated 100 times with different random seeds. In terms of accuracy, the proposed WMMD performs comparably with LOG and NNPU, showing advantages over DH and tADJ. When the true class-prior is less than equal to $0.8$, WMMD performs better in terms of AUC, except for tADJ. The tADJ achieves the highest AUC because $P(Y=1\mid X=x)$ is proportional to $P(\{x\text{ is from the positive dataset}\}\mid X=x)$. This empirically shows that WMMD has a comparable discriminant ability to the other algorithms for a wide range of class-priors.

Experiment 3: The main goal of this subsection is to show the robustness of the proposed classifier in the case of unknown class-prior ${\pi}_{+}$. In PU learning literature, $\pi_{+}$ has been frequently assumed to be known (du Plessis et al., 2015; Niu et al., 2016; Kiryo et al., 2017; Kato et al., 2019). However, this assumption can be considered to be strong in real-world applications, and to correctly execute existing PU learning algorithms, an accurate estimate of $\pi_{+}$ is necessary. In this experiment, we compare the accuracy and AUC when the class-prior ${\pi}_{+}$ is unknown. For the WMMD classifier, we used a density-based method for the class-prior estimation described in Appendix D.1, which can be obtained as a byproduct of the proposed algorithm. The results of LOG, DH, and NNPU are given for completeness sake using the ‘KM1’ method222While the ‘KM2’ method by Ramaswamy et al. (2016) is often considered to be a state-of-the-art method for estimating $\pi_{+}$, in our experiments, estimates based on the ‘KM2’ method have a larger estimation error than that of the ‘KM1’ method and thus we omitted it. by Ramaswamy et al. (2016). We take these estimates as true values and repeat the same comparative numerical experiments in Experiment 2.

Since the objective functions of the LOG, DH, and NNPU algorithms depend on the estimate $\hat{\pi}_{+}$, we anticipate that both the accuracy and AUC rely on the quality of the estimation. On the other hand, the tADJ algorithm does not depend on the class-prior, so the performance is not affected. Also, as the proposed score function does not depend on the class-prior $\pi_{+}$, and since $\pi_{+}$ is used only to determine a cutoff, the AUC of the proposed algorithm is less affected by the estimation of $\pi_{+}$.

Figures 3(a) and 3(c) compare the accuracy and AUC as a function of $n_{\rm u}$. WMMD performs worse than LOG, DH, and NNPU, while AUC is higher. Though tADJ shows poor accuracy in a wide range, it achieves high AUC comparable to WMMD. As we anticipated, WMMD is more robust than LOG, DH, and NNPU in AUC. This is possibly because our score function $\hat{\lambda}_{n_{\rm p},n_{\rm u}}$ does not depend on $\pi_{+}$. A similar trend can be found in Figures 3(b) and 3(d). We note that the ‘KM1’ method is not scalable and thus may not be used for large-scale datasets.

Experiment 4: In this experiment, we compare the elapsed training time, including hyperparameter optimization, of the five PU learning algorithms. The data are generated from the distributions described in Equation (9), and we set $n_{\rm p}=100,n_{\rm u}=400$, and $\pi_{+}=0.5$. The elapsed time is measured with 20 Intel® Xeon® E5-2630 v4@2.20GHz CPU processors.

Table 1 compares the elapsed training time and its ratio relative to that of WMMD. WMMD takes the shortest time among the five baseline methods. In particular, the training time for WMMD is at least about 300 times shorter than that of the LOG and DH methods. This is because the WMMD classifier has an analytic form while the LOG and DH methods require solving a non-linear programming problem.

6.2 Real data analysis

We demonstrate the practical utility of the proposed algorithm using the eight real binary classification datasets from the LIBSVM333https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ (Chang and Lin, 2011). Since some observations from the raw datasets are not completely recorded, we removed such observations and construct the dataset with fully recorded data. Next, to investigate the effect of varying $\pi_{+}$, we artificially reconstructed $\mathcal{X}_{\rm p}$ and $\mathcal{X}_{\rm u}$ through a random sampling from the fully recorded datasets. For the three datasets australian_scale, breast-cancer_scale and skin_nonskin, we reconstructed the data so that the resulting class-prior $\pi_{+}$ ranges from $0.15$ to $0.79$. We add the suffix 2 for those datasets. We randomly resampled data $100$ times for the seven small datasets and $10$ times for the four big datasets: skin_nonskin, skin_nonskin2, epsilon_normalized, and HIGGS. Table 2 summarizes statistics for the eleven real datasets. We conduct two comparative numerical experiments when $\pi_{+}$ is known and unknown.

Table 3 shows the average and the standard error of the accuracy and AUC when the class-prior $\pi_{+}$ is known. LOG and DH fail to compute the $(n_{\rm p}+n_{\rm u})$ $\times$ $(n_{\rm p}+n_{\rm u})$ Gram matrix due to out of memory in the 12 GB GPU memory limit. WMMD achieves comparable to or better accuracy and AUC than LOG, DH, and tADJ on most datasets. Compared to NNPU, WMMD performs comparably on the small datasets. However, NNPU achieves higher accuracy on skin_nonskin, epsilon_normalized, and HIGGS. The neural network used in NNPU fits well to the complicated and high-dimensional structure of data and shows high accuracy.

Table 4 compares the average and the standard error of the accuracy and AUC when the class-prior $\pi_{+}$ is unknown. As in Experiment 3 in Section 6.1, we estimate $\pi_{+}$ using the ‘KM1’ method for LOG, DH, and NNPU, and using the density-based method for WMMD. The LOG, DH, and NNPU algorithms are implemented on the seven small-scale datasets alone because the method by Ramaswamy et al. (2016) is not feasible with the large-scale datasets (Bekker and Davis, 2018). Overall, WMMD shows comparable to or better performances than other PU learning algorithms on most datasets. Compared to Table 3, WMMD and tADJ show robustness to unknown $\pi_{+}$ in terms of AUC. This is because WMMD and tADJ do not require estimation of $\pi_{+}$ to construct score functions. In contrast, the other methods require an estimate $\hat{\pi}_{+}$, and we observe a substantial drop in accuracy and AUC when the ‘KM1’ method estimate is used.

7 Concluding remarks

Existing methods use different objective functions and hypothesis spaces, and as a consequence, different optimization algorithms. Hence, there is no reason that one method outperforms uniformly for all scenarios. It is possible that one particular method may outperform in one scenario, for example, NNPU proposed by Kiryo et al. (2017) would perform better in complicated data settings because of the expressive power of neural networks. However, the proposed method has a clear computational advantage due to the closed-form as well as theoretical strength in terms of the explicit excess risk bound. Further, the proposed method works reasonably well in both cases in which $\pi_{+}$ is known or unknown. In this regard, we believe the proposed method can be used as a principled and easy-to-compute baseline algorithm in PU learning.

8 Acknowledgement

YK, WK, and MCP were supported by the National Research Foundation of Korea under grant NRF-2017R1A2B4008956. MS was supported by JST CREST JPMJCR1403.

Appendix A Proof of Theorem 3.1

Proof (Proof of Theorem 3.1).

Since a function $f\in\mathcal{F}$ is bounded by 1, we have $\ell_{\mathrm{h}}(yf(x))=\max(0,1-yf(x))=1-yf(x)$. Then, from Equation (2),

[TABLE]

Thus, we have

[TABLE]

Equation ($\ast$) holds because $\mathcal{F}$ is symmetric.

For the second result, note that a WIPM optimizer $g_{\mathcal{F}}$ satisfies ${\rm WIPM}(P_{X},P_{X\mid Y=1};2\pi_{+},\mathcal{F})$ $=\int_{\mathcal{X}}g_{\mathcal{F}}(x)dP_{X}(x)$ $-2\pi_{+}\int_{\mathcal{X}}g_{\mathcal{F}}(x)dP_{X\mid Y=1}(x)$. Thus, we have

[TABLE]

The last equality is from $R_{\ell_{\mathrm{h}}}(g_{\mathcal{F}})+R_{\ell_{\mathrm{h}}}(-g_{\mathcal{F}})=2$ due to $g_{\mathcal{F}}\in\mathcal{F}\subseteq\mathcal{M}$.

Appendix B Proofs for Section 3.3: Theoretical properties of empirical WIPM optimizer

In this section, we present a proof of Theorem 3.2 in Appendix B.2. We also provide a proof for Proposition 1 in Appendix B.3. Before presenting the proof for Theorem 3.2, we begin with necessary technical proposition and lemma in Appendix B.1.

B.1 Preliminaries for Theorem 3.2

In supervised binary classification settings, Sriperumbudur et al. (2012) introduced an empirical estimator for IPM and developed its consistency result. In Proposition 3, we recreate theoretical results for PU learning settings, giving a consistency result of empirical WIPM estimator.

Proposition 3 (Consistency result of WIPM estimator).

Let $\mathcal{F}$ be the symmetric function space such that $\mathinner{\!\left\lVert f\right\rVert}_{\infty}\leq\nu$, $\mathrm{Var}_{{P}_{X\mid Y=1}}(f)\leq\sigma_{X\mid Y=1}^{2}$, and $\mathrm{Var}_{{P}_{X}}(f)\leq\sigma_{X}^{2}$. Denote $\rho^{2}=\sigma_{X\mid Y=1}^{2}\vee\sigma_{X}^{2}$. Then for all $w,\alpha,\tau>0$, the following holds with probability at least $1-e^{-\tau}$ over the choice of $\mathcal{X}_{\rm pu}\mathrel{\mathop{\mathchar 58\relax}}=\{x_{1}^{\rm p},\dots,x_{n_{\rm p}}^{\rm p},x_{1}^{\rm u},\dots,x_{n_{\rm u}}^{\rm u}\}\sim P_{\rm pu}\mathrel{\mathop{\mathchar 58\relax}}=P_{X\mid Y=1}^{n_{\rm p}}\times P_{X}^{n_{\rm u}}$,

[TABLE]

Proof (Proof of Proposition 3).

The following proof is a slight modification of the proof of Theorem 3.3 in Sriperumbudur et al. (2012). Without loss of generality, by changing an order, we define a set of observations and a set of weights as follows,

[TABLE]

and

[TABLE]

respectively. For independent Rademacher random variables $\{\sigma_{i}\}_{i=1}^{n_{\rm p}+n_{\rm u}}$, we define the empirical Rademacher complexity-like term given by

[TABLE]

Note that $\tilde{\mathfrak{R}}_{\mathcal{X}_{\rm pu}}(\mathcal{F})\leq\mathfrak{R}_{\mathcal{X}_{\rm u}}(\mathcal{F})+w\mathfrak{R}_{\mathcal{X}_{\rm p}}(\mathcal{F})$. Define $\mu_{i}={P}_{X\mid Y=1}$ for $i\in\{1,\dots,n_{\rm p}\}$ and $\mu_{i}={P}_{X}$ for $i\in\{n_{\rm p}+1,\dots,n_{\rm p}+n_{\rm u}\}$, respectively. That is, $P_{\rm pu}=\times_{i=1}^{n_{\rm p}+n_{\rm u}}\mu_{i}$. Let $(X_{1},\dots,X_{n_{\rm p}+n_{\rm u}})\sim P_{\rm pu}$ and define random variables $\theta_{i}(f,X_{i})=w\{f(X_{i})-{P}_{X\mid Y=1}(f)\}/n_{\rm p}$ for $i\in\{1,\dots,n_{\rm p}\}$ and $\theta_{i}(f,X_{i})=\{f(X_{i})-{P}_{X}(f)\}/n_{\rm u}$ for $i\in\{n_{\rm p}+1,\dots,n_{\rm p}+n_{\rm u}\}$, respectively.

Then, using the fact that $\left|\sup|C|-\sup|D|\right|\leq\sup|C-D|$, we have

[TABLE]

Further, it is easy to verify that (i) $\int_{\mathcal{X}}\theta_{i}(f,z)d\mu_{i}(z)=0$ for all $i$ and $f\in\mathcal{F}$ and (ii) $\int_{\mathcal{X}}\theta_{i}^{2}(f,z)d\mu_{i}(z)$ is bounded by $w^{2}{\sigma^{2}_{X\mid Y=1}}/{n_{\rm p}^{2}}$ for $i\in\{1,\dots,n_{\rm p}\}$ and ${\sigma^{2}_{X}}/{n_{\rm u}^{2}}$ for $i\in\{n_{\rm p}+1,\dots,n_{\rm p}+n_{\rm u}\}$, respectively. Finally, (iii) $\mathinner{\!\left\lVert\theta_{i}(f,\cdot)\right\rVert}_{\infty}\leq 2(w/n_{\rm p}+1/n_{\rm u})\nu=\chi_{n_{\rm p},n_{\rm u}}^{(2)}(w)\nu$ for all $i\in\{1,\dots,n_{\rm p}+n_{\rm u}\}$.

Then, for all $\alpha>0$, the following holds with probability at least $1-e^{-\tau}$,

[TABLE]

The first inequality is derived by the second inequality of Lemma 3 in Appendix B.1, a variant of the Talagrand’s inequality. The second inequality is from using a symmetrization lemma: with corresponding independent ghost empirical distributions $\tilde{P}_{X,n_{\rm u}}$ and $\tilde{P}_{X\mid Y=1,n_{\rm p}}$,

[TABLE]

Next, simply using the fact $\tilde{\mathfrak{R}}_{\mathcal{X}_{\rm pu}}(\mathcal{F})\leq\mathfrak{R}_{\mathcal{X}_{\rm u}}(\mathcal{F})+w\mathfrak{R}_{\mathcal{X}_{\rm p}}(\mathcal{F})$, we have for all $w,\alpha,\tau>0$, the following holds with probability at least $1-e^{-\tau}$,

[TABLE]

It concludes the proof.

For Lemmas 3, we quote the Proposition B.1 of Sriperumbudur et al. (2012) without proofs.

Lemma 3 (Proposition B.1 of Sriperumbudur et al. (2012): A variant of Talagrand’s inequality).

Let $B\geq 0,n\geq 1,(\Omega_{i},\mathcal{A}_{i},\mu_{i}),i=1,\dots,n$ be a probability space and $\theta_{i}\mathrel{\mathop{\mathchar 58\relax}}\mathcal{F}\times\Omega_{i}\to\mathbb{R}$ be bounded measurable functions, where $\mathcal{F}$ is the space of real-valued $\mathcal{A}_{i}$-measurable functions for all $i$. Suppose

• (a)

$\int_{\Omega_{i}}\theta_{i}(f,\omega)d\mu_{i}(\omega)=0$* for all $i$ and $f\in\mathcal{F}$*

• (b)

$\int_{\Omega_{i}}\theta_{i}^{2}(f,\omega)d\mu_{i}(\omega)\leq\rho_{i}^{2}$* for all $i$ and $f\in\mathcal{F}$*

• (c)

$\mathinner{\!\left\lVert\theta_{i}(f,\cdot)\right\rVert}_{\infty}\leq B$* for all $i$ and $f\in\mathcal{F}$.*

Define $Z\mathrel{\mathop{\mathchar 58\relax}}=\times_{i=1}^{n}\Omega_{i}$ and $P\mathrel{\mathop{\mathchar 58\relax}}=\times_{i=1}^{n}\mu_{i}$. Furthermore, define $g\mathrel{\mathop{\mathchar 58\relax}}Z\to\mathbb{R}$ by

[TABLE]

Then, for all $\tau>0$, we have

[TABLE]

In addition, for all $\tau>0$ and $\alpha>0$,

[TABLE]

B.2 Proof of Theorem 3.2: estimation error bound of WIPM optimizer

Proof (Proof of Theorem 3.2).

We first define the empirical risk estimator $\hat{R}_{\ell_{\mathrm{h}}}(f)$ by replacing data distributions in Equation (2) with the empirical distributions. To be more specific, we define

[TABLE]

Since $\mathinner{\!\left\lVert f\right\rVert}_{\infty}\leq 1$, using the similar derivations in Appendix A, we have

[TABLE]

By the result of Theorem 3.1, the negative of an WIPM optimizer ${f}_{\mathcal{F}}$ is minimizer of $R_{\ell_{\mathrm{h}}}(f)$, i.e., $R_{\ell_{\mathrm{h}}}({f}_{\mathcal{F}})=\inf_{f\in\mathcal{F}}R_{\ell_{\mathrm{h}}}(f)$. Thus, we have

[TABLE]

The first inequality holds since $\hat{R}_{\ell_{\mathrm{h}}}(\hat{f}_{\mathcal{F}})=1-{\rm WIPM}({P}_{X,n_{\rm u}},{P}_{X\mid Y=1,n_{\rm p}};2\pi_{+},\mathcal{F})\leq\hat{R}_{\ell_{\mathrm{h}}}(f)$ for any $f\in\mathcal{F}$. Thus it is enough to bound $\sup_{f\in\mathcal{F}}|\hat{R}_{\ell_{\mathrm{h}}}({f})-R_{\ell_{\mathrm{h}}}({f})|$, and

[TABLE]

Note that this is a special case of Equation (B.1). Therefore, applying Equation (10) in Proposition 3 with $w=2\pi_{+}$, we have for all $\alpha,\tau>0$, the following holds with probability at least $1-e^{-\tau}$,

[TABLE]

This concludes a proof.

B.3 Proof of Proposition 1

Proof (Proof of Proposition 1).

After omitting the positive term $\Delta$ from the upper bound (6) and plugging $\alpha=1$ in the upper bound (5), it is enough to show that, under the condition (7), the following is satisfied.

[TABLE]

where $R=\mathbb{E}_{P_{X}^{n_{\rm u}}}\{\mathfrak{R}_{\mathcal{X}_{\rm u}}(\mathcal{F})\}+2\pi_{+}\mathbb{E}_{P_{X\mid Y=1}^{n_{\rm p}}}\{\mathfrak{R}_{\mathcal{X}_{\rm p}}(\mathcal{F})\}$.

Using simple algebras, we have

[TABLE]

Thus, the proposed upper bound is sharper than that of Kiryo et al. (2017) if the condition (7) holds.

Appendix C Proofs for Section 4: The empirical WMMD optimizer and the WMMD classifier

In this section, we first show that the empirical WMMD optimizer has a closed-form expression in Appendix C.1. We provide proofs of Lemmas 1 and 2 in Appendix C.2 and Theorem 4.1 in Appendix C.3.

C.1 Proof of Proposition 2: WMMD optimizer has a closed-form expression

We first state and prove the following Proposition 4, which is an extended version of Proposition 2. Please note that Proposition 2 can be directly obtained by plugging the two empirical distributions ${P}_{X\mid Y=1,n_{\rm p}}$ and ${P}_{X,n_{\rm u}}$ and $w=2\pi_{+}$.

Proposition 4 (Weighted maximum mean discrepancy).

Let ${P}$ and ${Q}$ be two probability measures defined on $\mathcal{X}$ and let $k\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ a bounded reproducing kernel.

(a) WMMD between two probability measures ${P}$ and ${Q}$ with a weight $w$ and a closed ball $\mathcal{H}_{k,r}$ with the radius $r$ can be represented in a closed form,

[TABLE]

where $x$ and $x^{\prime}$ independently follow ${P}$ and $y$, and $y^{\prime}$ independently follow ${Q}$.

(b) we also have a closed-form expression for the unique optimizer $g_{\mathcal{H}_{k,r}}\in\mathcal{H}_{k,r}$ given by

[TABLE]

where $T$ is the normalizing operator defined by $T(g)=g/{\mathinner{\!\left\lVert g\right\rVert}_{\mathcal{H}_{k}}}$.

(c) The associated classifier is given by

[TABLE]

where

[TABLE]

Proof (Proof of Proposition 4).

The main idea of this proof is to use $f(x)={\langle f,k(\cdot,x)\rangle}_{\mathcal{H}_{k}}$ as Gretton et al. (2012) showed. From the definition of the WMMD, we have

[TABLE]

The last equation is obtained by Cauchy-Schwarz inequality with the WMMD optimizer $g_{\mathcal{H}_{k,r}}\in\mathcal{H}_{k,r}$ given by

[TABLE]

It concludes a proof for (b).

Furthermore,

[TABLE]

It concludes a proof for (a).

Lastly, we prove the statement (c). Note that the associated classifier is determined by the sign of $f_{\mathcal{H}_{k,r}}=-g_{\mathcal{H}_{k,r}}$. From the statement (b), we have

[TABLE]

Hence, we have

[TABLE]

C.2 Lemmas for Section 4

In this subsection, we provide the proof for the two useful Lemmas: (i) explicitly showing the estimation error bound in Lemma 1 and (ii) deriving an approximation error bound in Lemma 2.

Proof (Proof of Lemma 1).

We first prove $\mathcal{H}_{k,r_{1}}\subseteq\mathcal{M}$. By the reproducing property of $\mathcal{H}_{k}$ and Cauchy-Schwarz inequality, for any $f\in\mathcal{H}_{k,r_{1}},x\in\mathcal{X}$

[TABLE]

Thus, $\mathinner{\!\left\lVert f\right\rVert}_{\infty}\leq 1$ for all $f\in\mathcal{H}_{k,r_{1}}$. This proves $\mathcal{H}_{k,r_{1}}\subseteq\mathcal{M}$.

Now, we prove the inequality. First, we apply Theorem 3.2 with $\mathcal{H}_{k,r_{1}}$.

[Step 1] From the result of Theorem 3.2, for all $\alpha,\tau>0$, the estimation error term is bounded above with probability at least $1-e^{-\tau}$,

[TABLE]

Using the notations $(\tilde{y}_{1},\dots,\tilde{y}_{n_{\rm p}},\tilde{y}_{n_{\rm p}+1},\dots,\tilde{y}_{n_{\rm p}+n_{\rm u}})$ $=(2\pi_{+}/n_{\rm p},$ $\dots,$ $2\pi_{+}/n_{\rm p},$ $-1/n_{\rm u},\dots,-1/n_{\rm u})$, we obtain upper bound of the empirical Rademacher complexity of $\mathcal{H}_{k,r}$ given the positive samples $\mathfrak{R}_{\mathcal{X}_{\rm p}}(\mathcal{H}_{k,r})$ as follows.

[TABLE]

We continue the similar method to the unlabeled dataset, and applying expectation operator gives

[TABLE]

[Step 2] Using Equations (12) and (13), we then conclude that for all $\alpha,\tau>0$, the following holds with probability at least $1-e^{-\tau}$,

[TABLE]

These equations conclude the proof.

Proof (Proof of Lemma 2).

[Step 1] In this step, we first claim that $\inf_{f\in\beta\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)=R_{\ell_{\mathrm{h}}}(\beta f^{*}_{1})$. By Lin (2002, Lemma 3.1), the $f^{*}_{1}$ satisfies that $\inf_{f\in\mathcal{U}}R_{\ell_{\mathrm{h}}}(f)=\inf_{f\in\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)$ $=R_{\ell_{\mathrm{h}}}(f^{*}_{1})$. Note that $\ell_{\mathrm{h}}(yf(x))=\max(0,1-yf(x))=1-yf(x)$ for all $\mathinner{\!\left\lVert f\right\rVert}_{\infty}\leq 1$. It is obvious that $\beta f_{1}^{*}\in\beta\mathcal{M}$. By definition of the infimum, we have $R_{\ell_{\mathrm{h}}}(\beta f^{*}_{1})\geq\inf_{f\in\beta\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)$. Suppose $R_{\ell_{\mathrm{h}}}(\beta f^{*}_{1})>\inf_{f\in\beta\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)$. Let $f_{\beta}^{*}$ be a function in $\beta\mathcal{M}$ such that $R_{\ell_{\mathrm{h}}}(f^{*}_{\beta})=\inf_{f\in\beta\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)$ Then,

[TABLE]

Since $\beta^{-1}f_{\beta}^{*}\in\mathcal{M}$,

[TABLE]

Note that $\inf_{f\in\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)=R_{\ell_{\mathrm{h}}}(f^{*}_{1})$ This contradicts with the assumption $R_{\ell_{\mathrm{h}}}(\beta f^{*}_{1})>\inf_{f\in\beta\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)$, and we have $\inf_{f\in\beta\mathcal{M}}R_{\ell_{\mathrm{h}}}(f)=R_{\ell_{\mathrm{h}}}(\beta f^{*}_{1})$.

[Step 2] By Lemma 1, we have $\mathcal{H}_{k,r_{1}}\subseteq\mathcal{M}$. Thus, for all $g\in\mathcal{H}_{k,r_{1}}\subseteq\mathcal{M}$ and $0<\beta\leq 1$, we have

[TABLE]

The first equality holds because $\ell_{\mathrm{h}}(yg(x))=\max(0,1-yg(x))=1-yg(x)$ for all $g\in\mathcal{H}_{k,r_{1}}\subseteq\mathcal{M}$. Hence,

[TABLE]

Therefore, by [Step 1] and [Step 2],

[TABLE]

for any $0<\beta\leq 1$.

C.3 Proof of Theorem 4.1

We first state the following assumptions for Theorem 4.1.

• (A1)

The distribution functions $P_{X}$ and $P_{X\mid Y=1}$ have probability density functions $p_{X}(x)$ and $p_{X\mid Y=1}(x)$, respectively.

• (A2)

For $\alpha_{\text{H}}\in(0,1)$, the density functions $p_{X}(x)$ and $p_{X\mid Y=1}(x)$ are $\alpha_{\text{H}}$-Hölder continuous, i.e., $|g(x)-g(x^{\prime})|\leq C_{\text{H}}|x-x^{\prime}|^{\alpha_{\text{H}}}$ for some constant $C_{\text{H}}$ and for all $x,x^{\prime}\in\mathbb{R}^{d}$ and $g\in\{p_{X},p_{X\mid Y=1}\}$.

• (A3)

The marginal density function is bounded away from zero, i.e., $p_{X}(x)\geq p_{\rm min}>0$ for all $x$ on its support.

• (A4)

The marginal distribution $P_{X}$ has Tsybakov’s noise exponent $q\in[0,\infty)$, i.e., there exists a constant $C_{\text{noise}}>0$ such that for all sufficiently small $t>0$, we have

[TABLE]

where $\eta(x)\mathrel{\mathop{\mathchar 58\relax}}=P(Y=1\mid X=x)$.

To begin, we quote a useful Lemma from Jiang (2017, Theorem 2) into our contexts. We let $\tilde{f}_{\mathcal{H}_{k}}(z)=\frac{2\pi_{+}}{\sqrt{2\pi}n_{\rm p}h^{d}}\sum_{i=1}^{n_{\rm p}}k(z,x_{i}^{{\rm p}})-\frac{1}{\sqrt{2\pi}n_{\rm u}h^{d}}$ $\sum_{i=1}^{n_{\rm u}}k(z,x_{i}^{{\rm u}})$ and $\tilde{f}_{{\rm Bayes}}(z)=2\pi_{+}p_{X\mid Y=1}(z)-p_{X}(z)$.

Lemma 4 (Theorem 2 of Jiang (2017)).

Suppose $p_{X}(x)$ and $p_{X\mid Y=1}(x)$ are $\alpha_{\text{H}}$-Hölder continuous. Then there exist a constant $\tilde{C}_{(\rm p,\rm u)}$ such that the following holds with probability at least $1-1/n_{\rm p}-1/n_{\rm u}$.

[TABLE]

where the bandwidth $h=(n_{\rm p}\wedge n_{\rm u})^{-\frac{1}{2\alpha_{\text{H}}+d}}$.

Proof (Proof of Theorem 4.1).

Due to the $\ell_{01}(z)=\ell_{01}(cz)$ for any $c>0$ and $z\in\mathbb{R}$, we have

[TABLE]

where $\mathbbm{1}_{\{\cdot\}}$ denotes the indicator function. The second equality is due to

[TABLE]

Denote $\varepsilon(n_{\rm p},n_{\rm u})\mathrel{\mathop{\mathchar 58\relax}}=4\pi_{+}\tilde{C}_{(\rm p,\rm u)}(n_{\rm p}\wedge n_{\rm u})^{-\frac{\alpha_{\text{H}}}{2\alpha_{\text{H}}+d}}$ and let $\mathcal{E}$ be the event that Equation (14) holds. By Lemma 4, under the event $\mathcal{E}$, $\mathbbm{1}_{\{\tilde{f}_{\mathcal{H}_{k}}(X)\tilde{f}_{{\rm Bayes}}(X)<0\}}\leq\mathbbm{1}_{\{|{\tilde{f}_{{\rm Bayes}}(X)}|\leq\varepsilon(n_{\rm p},n_{\rm u})\}}$. Thus, we have

[TABLE]

for some constant $C_{(\rm p,\rm u)}>0$. Thus, under the assumptions (A1)-(A4), we have the following with probability at least $1-1/n_{\rm p}-1/n_{\rm u}$:

[TABLE]

Appendix D Implementation details

In this section, we provide implementation details for WMMD and other baseline PU learning algorithms in Appendix D.1 and D.2, respectively.

D.1 Proposed PU learning algorithm: WMMD algorithm

We divided the original training data into training and validation sets, with an 80-20 random split. Let $\tilde{x}_{j}^{\rm p}$ and $\tilde{x}_{j}^{\rm u}$ be the positive and the unlabeled samples in validation set, respectively. Similarly, $\tilde{n}_{\rm p}$ and $\tilde{n}_{\rm u}$ be the number of samples in the positive and the unlabeled validation set, respectively. With the validation set, we conducted a grid search method for the hyperparameter selection with a grid $\gamma\in\{1,0.4,0.2,0.1,0.05\}$ for all the numerical experiments. With the grids, we selected the optimal hyperparameters $\gamma^{*}$ which minimized

[TABLE]

where

[TABLE]

[TABLE]

and $k_{\gamma}(z_{1},z_{2})=\exp(-\gamma\mathinner{\!\left\lVert z_{1}-z_{2}\right\rVert}_{2}^{2})$. Note that Equation (15) is an empirical estimation of the misclassification error since

[TABLE]

Final classification for a test datum $z$ was determined by ${\rm sign}(\hat{f}_{\mathcal{H}_{k,r}}^{\gamma^{*}}(z))$ and AUC was computed by using $\hat{\lambda}_{n_{\rm p},n_{\rm u}}^{\gamma^{*}}(z)$.

When the class-prior $\pi_{+}$ is unknown, we suggest a simple $\pi_{+}$ estimation method, called the density-based method, to find $\hat{\pi}_{+}^{\rm WMMD}$ such that for some $\eta\in(0,1)$,

[TABLE]

This estimator is sensible because $\pi_{+}\leq p_{X}(x)/p_{X\mid Y=1}(x)$ and $\{\hat{\lambda}_{n_{\rm p},n_{\rm u}}^{\gamma}(x)\}^{-1}$ can be considered as a kernel density estimation of $p_{X}(x)/p_{X\mid Y=1}(x)$ for $x\in supp(P_{X\mid Y=1})$. Here, we denote the density functions by $p_{X}$ and $p_{X\mid Y=1}$. In our experiments, we fix $\eta=0.1$ and using $\hat{\pi}_{+}^{\rm WMMD}$ leaded better performance than using ‘KM1’ method in terms of accuracy.

D.2 The baseline PU learning algorithms

We compared the following 4 PU learning algorithms: (i) the logistic loss $\ell_{\mathrm{log}}$, denoted by LOG, (ii) the double hinge loss $\ell_{\mathrm{dh}}$, denoted by DH, both proposed by du Plessis et al. (2015), (iii) the non-negative risk estimator method, denoted by NNPU, proposed by Kiryo et al. (2017), and (iv) the threshold adjustment method, denoted by tADJ, proposed by Elkan and Noto (2008).

**General: ** Similar to the procedure in Appendix D.1, we set training and validation sets with the 80-20 random split of the original training dataset and we conducted a grid search method for hyperparameter selection.

**LOG and DH: ** As du Plessis et al. (2015) proposed, we followed a binary discriminant function as

[TABLE]

where $N=n_{\rm p}+n_{\rm u}$ and $\varphi_{i}(x)=\exp(-\gamma\mathinner{\!\left\lVert x-c_{i}\right\rVert}_{2}^{2})$ for $\{c_{1},\dots,c_{N}\}=\{x_{1}^{\rm p},\dots,x_{n_{\rm p}}^{\rm p},$ $x_{1}^{\rm u},\dots,x_{n_{\rm u}}^{\rm u}\}$. For a loss function $\ell\in\{\ell_{\mathrm{log}},\ell_{\mathrm{dh}}\}$, the empirical risk function is given by

[TABLE]

Here, the hyperparameter grids are $\lambda\in\{1,0.4,0.2,0.1,0.05\}$ and $\gamma\in\{1,0.4,0.2,$ $0.1,0.05\}$ for all the numerical experiments. With the grids, we selected the optimal hyperparameter $(\lambda^{*},\gamma^{*})$ which minimized the empirical risk on the validation set $\tilde{J}_{\lambda,\gamma}(\alpha,b)$ defined by

[TABLE]

After selecting the optimal hyperparameter $(\lambda^{*},\gamma^{*})$, we minimized $\hat{J}_{(\lambda^{*},\gamma^{*})}(\alpha,b)$ with the gradient descent algorithm. Learning rate was fixed by $0.1$ and the number of epochs was $100$. During the training, we applied the early stopping rule: we stopped training if the validation error is not minimized in 10 successive epochs. After the training phase, with the trained $\hat{\alpha}$ and $\hat{b}$, we classified a test datum $z$ as ${\rm sign}(g_{\hat{\alpha},\hat{b}}(z))$. AUC was computed by using $g_{\hat{\alpha},\hat{b}}(z)$.

**NNPU: ** We followed the method by Kiryo et al. (2017). The model for NNPU was a 5-layer multilayer perceptron with ReLU nonlinearity ($d$-300-300-300-1). We applied the batch normalization before each ReLU nonlinearity. Please note that this network architecture is quite similar to the model in Kiryo et al. (2017). We used a stochastic gradient descent algorithm with a learning rate $0.01$. Loss function was the sigmoid function. The number of epochs was $100$ and the optimal weights were selected at the best validation error during the training.

**tADJ: ** We followed the method by Elkan and Noto (2008). We used ‘LogisticRegressionCV’ function in the Python module ‘sklearn.linear_model’ (Pedregosa et al., 2011) to estimate $P(\{x\text{ is from the positive dataset}\}\mid X=x)$. The hyperparameter grid for $L_{2}$-regularizer was $\{0.01,0.1,1,10,100\}$ and the optimal hyperparameter was chosen based on 5-fold cross validation on the split training dataset, i.e., 80% of the original training dataset. Then, $P(\{x\text{ is from the positive dataset}\}\mid Y=1)$ was estimated by the split validation set, i.e., 20% of the original training dataset.

Appendix E Comparison between Gaussian and inverse kernels

We compared LOG, DH, and WMMD using two kernels: (i) the Gaussian kernel $k(x,y)=\exp(-\gamma\mathinner{\!\left\lVert x-y\right\rVert}_{2}^{2})$ and (ii) the inverse kernel $k(x,y)=\frac{\gamma}{\gamma+\mathinner{\!\left\lVert x-y\right\rVert}_{2}^{2}}$ for $\gamma>0$. Figures 4(a) and 4(c) show the accuracy and AUC of LOG, DH, and WMMD on various $n_{\rm u}$. The training sample size for the positive data is $n_{\rm p}=100$ and the class prior is $\pi_{+}=0.5$. The unlabeled sample size changes from 40 to 500 by 20. We repeat a random generation of training and test data 100 times. For comparison purposes, we add the $1-$Bayes risk for each unlabeled sample size. In every algorithm, using the Gaussian kernel achieves higher accuracy and AUC than using the inverse kernel in every $n_{\rm u}$.

Figures 4(b) and 4(d) show comparison of the accuracy and AUC as $\pi_{+}$ changes. The training sample size for the positive and the unlabeled data are $n_{\rm p}=100$ and $n_{\rm u}=400$, respectively. The class-prior $\pi_{+}$ changes from 0.05 to 0.95 by 0.05. The test sample size is $10^{3}$. We repeat a random generation of training and test data 100 times. Both kernels perform comparably for LOG and WMMD algorithm. The DH algorithm with the inverse kernel achieves higher accuracy when the class-prior is close to 0.5. But in terms of AUC, both kernels perform comparably in every $\pi_{+}$.