Monotone single index models with general non-symmetric design

We assume the corrupted sample $Y$ is independent of $X$ given $\tilde{Y}$, which corresponds to a missing-at-random assumption. The corruption structure is completely specified by a $2\times 2$ confounding matrix $Q$ such that $Q_{k+1,l+1}=\mathbb{P}(Y=k|\tilde{Y}=l)$, whose off-diagonal terms are smaller than diagonals. Under these assumptions, $\mathbb{E}[Y|X]=\mathbb{P}(Y=1|X)$ is given by

where we let $g_{\star}(s):=q(s)Q_{22}+(1-q(s))Q_{21}$. Note $g_{\star}$ is monotonic since $g_{\star}(s)=q(s)(Q_{22}-Q_{21})+Q_{21}$ and we assume that $Q_{22}>Q_{21}$. In many practical situations the confounding matrix $Q$ is unknown, which limits the parametric approach even with the assumption of a specific $q$, such as a sigmoid function $q(u)=1/(1+\exp(-u)$.

More concretely, let $(\tilde{X},\tilde{Y})\in\mathbb{R}^{p}\times\mathbb{R}$ be a random sample from the population. As in Example 2, we assume $\tilde{Y}$ is related to $\tilde{X}$ via $\mathbb{E}[\tilde{Y}|\tilde{X}]=g(\tilde{X}^{T}u_{*})$ for some unknown parameter vector $u_{*}\in\mathbb{R}^{p}$ and a monotonic function $g$. Instead of $(\tilde{X},\tilde{Y})$, we observe the pair $(X,Y)$ which is defined as follows. If a sample is labeled ($Y=1$), an associated response is positive ($\tilde{Y}=1$) and the sample is taken from positive sub-population, i.e. $X\overset{d}{=}(\tilde{X}|\tilde{Y}=1$). On the contrary, if a sample is unlabeled ($Y=0$), we know that the sample is randomly drawn from the population, i.e. $X\overset{d}{=}\tilde{X}$, but an associated outcome is unknown.

Regarding $Y$ as a corrupted version of $\tilde{Y}$, the PU problem can be viewed as a special case of Example 2, where only positive samples $(\tilde{Y}=1)$ are corrupted. In fact, the confounding matrix Q is given in the following form :

where $\gamma$ is the ratio of the number of positive to unlabeled samples. Therefore, following the same lines in (1.1) and noting that $\tilde{Y}|X\overset{d}{=}\tilde{Y}|\tilde{X}$, we have $\mathbb{E}[Y|X]=g_{\star}(X^{T}u_{\star})$ for $g_{\star}$ defined in Example 2 with $Q$ defined in (1.3).

We note that the conditional mean function $g_{\star}$ in (1.1) is known up to a contamination proportion $\pi=\mathbb{P}(\tilde{Y}=1)$. If the contamination proportion is known from alternative source, we may use a parametric approach to estimate $u_{\star}$, as in . However, if the contamination proportion is unknown, reliable estimation of such proportion is very challenging and heavily relies on model specification ([Hastie2013-wq]). The single index model framework translates the problem of unknown $\pi$ to unknown $g_{\star}$ and consequently provides a framework for addressing the PU-learning problem with unknown $\pi$.

Another interesting aspect of this PU problem is that the distribution of $X$ cannot be symmetric regardless of distribution of original design $\tilde{X}$. This is because $X$ has a mixture distribution of two distributions $\tilde{X}|\tilde{Y}=1$ and $\tilde{X}$, and by the monotonicity assumption of $g$, the distribution of $\tilde{X}|\tilde{Y}=1$ cannot be symmetric and precludes existing methods relying on the symmetric design assumption.

The single index model has been an active area in statistics for a number of decades, since its introduction in Econometrics and Statistics [Ichimura1993-gx, Horowitz1996-jv, Klein1993-is]. There is an extensive literature about estimation of the index vector $u_{\star}$, which can be broadly classified into two categories : methods which directly estimate the index vector while avoiding a “nuisance’ ’infinite-dimensional function, and methods which involve estimation of both the infinite-dimensional function and the index vector.

Methods of the first type usually require strong design assumption such as Gaussian or an elliptically symmetric distribution. [Duan1991-nk] proposed an estimator of $u_{\star}$ using sliced inverse regression in the context of sufficient dimension reduction. $\sqrt{n}$ consistency and asymptotically normality of the estimator are established in the low-dimensional regime, under the elliptically symmetrical design assumption. This result is extended in [Lin2018-rk] in high-dimensional regime under similar design assumptions. [Plan2013-ew] studied sparse recovery of the index vector $u_{\star}$ when $X$ is Gaussian, and established a non-asymptotic $\mathcal{O}((s\log(2p/s)/n)^{1/2})$ mean-squared error bound, where $s$ is a sparsity level, i.e. $s:=\|u_{\star}\|_{0}$. [Ai2014-lo] studied the error bound of the same estimator when the design is not Gaussian and showed that the estimator is biased and the bias depends on the size of $\|u_{\star}\|_{\infty}$. The proposed estimator in [Plan2013-ew] is generalized in several directions; [Yang2017-cl] propose an estimator based on the score function of the covariate, which relies on prior knowledge of covariate distribution. [Chen2017-py] proposed estimators based on U-statistics which also include [Plan2013-ew] as a special case with standard Gaussian design assumption.

Methods of the second type usually do not require strong design assumption, and our method also falls into this category. One popular approach is via M-estimation or solving estimation equations. Under smoothness assumptions of $g_{\star}$, methods for a link function estimation were proposed via Kernel estimation([Ichimura1993-gx, Hardle1993-oj, Klein1993-is, Delecroix2003-fv, Cui2011-go]), local polynomials ([Carroll1997-iz, Chiou1998-zw]) , or splines ([Yu2002-hq, Kuchibhotla2016-ef]). An estimator for the index vector is then obtained by minimizing certain criterion function such as squared loss ([Ichimura1993-gx, Hardle1993-oj, Yu2002-hq]) or solving an estimating equation ([Chiou1998-zw, Cui2011-go]). $\sqrt{n}$ consistency and asymptotic normality of those proposed estimators were established for those estimators in the low-dimensional regime. Another approach is the average derivative method ([Powell1989-je, Hristache2001-om]), which takes advantage of the fact $\frac{d}{dx}E[Y|X]=u_{\star}g^{\prime}(X^{T}u_{\star})$. An estimator is then obtained through a non-parametric estimation of $\frac{d}{dx}E[Y|X]$, also under smoothness assumption of $g_{\star}$. There are also studies for single index model in high-dimensional regime, using PAC-Bayesian approach ([Alquier2013-ij]) or incorporating penalties ([Wang2012-mb, Radchenko2015-wa] )

There have been an increasing number of studies where the link function is restricted to have a particular shape, such as a monotone shape, but is not necessarily a smooth function. [Kalai2009-ew] and [Kakade2011-yk] investigated single index problem in low-dimensional regime under monotonic and Lipschitz continuous link assumption, and proposed iterative perceptron-type algorithms with prediction error guarantees. In particular, an isotonic regression is used for the estimation of $g_{\star}$ then update $u_{\star}$ based on estimated $g_{\star}$. [Ganti2017-dq] extended [Kalai2009-ew, Kakade2011-yk] to the high-dimensional setting via incorporating a projection step onto a sparse subset of the parameter space. They empirically demonstrated good predictive performance, but no theoretical guarantee were provided. [Balabdaoui2016-sd] study the least square estimators under monotonicity constraint and establish $n^{1/3}$ consistency of the least-square estimator. Also in [Balabdaoui2017-at], a $\sqrt{n}$ consistent estimator for the index vector $u_{\star}$, based on solving a score function, is proposed.

Our work focuses on estimation of the index vector in high-dimensional regime where an unknown function is assumed to be Lipschitz-continuous and monotonic. We provide an efficient algorithm via iterative hard-thresholding and a non-asymptotic $\ell_{2}$ and mean function error bound.

Our major contributions can be summarized as follows:

• •

  Develop a scalable projection-based iterative approach. For each iteration $u=u^{t}$, we first perform isotonic regression of $Y$ on to $Xu^{t}$ using the PAVA algorithm, taking an orthogonal gradient step to update $u$, and then enforcing sparsity using a thresholding operator. The algorithm is stopped once the mean-squared error drops below a threshold $\tau$.

• •

  The main contribution is to provide theoretical guarantees for our algorithm, both for the parameter vector $u_{\star}$ and the mean function $g_{\star}(X^{T}u_{\star})$ that scales as $\biggr{(}\frac{s\log p}{n}\biggr{)}^{1/3}$. Note that the minimax rate for low-dimensional isotonic regression is $n^{-1/3}$ and our algorithm achieves this rate up to $\log p$.

• •

  Since our goal is largely to estimate the parameter vector $u_{\star}$, we also show that under additional assumptions including a $\beta$-min condition, we show that a one-step correction achieves the parametric low-dimensional $\frac{1}{\sqrt{n}}$ rate. The one-step correction is motivated by the earlier work of [Balabdaoui2017-at] which proves $\frac{1}{\sqrt{n}}$ rates in the low-dimensional case.

• •

  Finally we support our theoretical results with a simulation study which supports our theoretical rates and also shows that our algorithm compares favorable to existing approaches.

The remainder of the paper is organized as follows: Section LABEL:SecProb provides a problem setup and we develop our algorithm; in Section 3 we provide theoretical guarantees for our algorithm which are supported by experimental results in Section 4. Proofs are provided in Section 5.