\mathrm{\RNum{1}}(\Delta)=&\frac{1}{n}\sum_{i=1}^n (x_i^\top \Delta)\nabla f(\Theta^\top x_i)^\top \left[\nabla A(f(\Theta^\top x_i))-\nabla A(f(\Theta^{*\top}x_i))\right]\\
=&\frac{1}{n}\sum_{i=1}^n (\Delta^\top x_i)^\top\nabla f(\Theta^\top x_i)^\top \nabla^2 A(f(\Theta_t^\top x_i))\nabla f(\Theta_t^\top x_i)(\Delta^\top x_i)\\
\geq&\frac{1}{n}\sum_{i=1}^n (\Delta^\top x_i)^\top\nabla f(\Theta_t^\top x_i)^\top \nabla^2 A(f(\Theta_t^\top x_i))\nabla f(\Theta_t^\top x_i)(\Delta^\top x_i)\\
&-\frac{1}{n}\sum_{i=1}^n \|\Delta^\top x_i\|_2^2\left\|\nabla f(\Theta_t^\top x_i)^\top\nabla^2 A(f(\Theta_t^\top x_i))\left[\nabla f(\Theta_t^\top x_i)-\nabla f(\Theta^\top x_i)\right]\right\|,
\end{split}
\end{equation}

where $\Theta_t=t\Theta+(1-t)\Theta^*$ lies between $\Theta$ and $\Theta^*$.
By Lemma <ref>, we have

\begin{equation}
\begin{split}
&(\Delta^\top x_i)^\top\nabla f(\Theta_t^\top x_i)^\top \nabla^2 A(f(\Theta_t^\top x_i))\nabla f(\Theta_t^\top x_i)(\Delta^\top x_i)\\
\geq &(h^{\MN}(R)+4C_xR)\|\Delta^\top x_i\|_2^2.
\end{split}
\end{equation}

On the other hand, since $\nabla f(v)=I_K-1_K\nabla A(v)^\top$, one can show that

\begin{equation}
\begin{split}
&\left\|\nabla f(\Theta_t^\top x_i)^\top\nabla^2 A(f(\Theta_t^\top x_i))\left[\nabla f(\Theta_t^\top x_i)-\nabla f(\Theta^\top x_i)\right]\right\|\\
\leq &\left\|\nabla f(\Theta_t^\top x_i)^\top\nabla^2 A(f(\Theta_t^\top x_i))1_K\right\|_2\left\|\nabla A(\Theta_t^\top x_i)-\nabla A(\Theta^\top x_i)\right\|\\
\leq&2\left\|\nabla A(\Theta_t^\top x_i)-\nabla A(\Theta^\top x_i)\right\|\\
\leq&4C_xR,
\end{split}
\end{equation}

where the third line is due to that

\begin{equation}
\begin{split}
&\left\|\nabla f(\Theta_t^\top x_i)^\top\nabla^2 A(f(\Theta_t^\top x_i))1_K\right\|_2\\
\leq&\|\nabla^2 A(f(\Theta_t^\top x_i))1_K\|_2+1_K^\top \nabla^2 A(f(\Theta_t^\top x_i))1_K\|\nabla A(\Theta_t^\top x_i)\|_2\\
\leq&1+\|\nabla A(\Theta_t^\top x_i)\|_2\\
\leq&2,
\end{split}
\end{equation}

and the last line is due to that

\begin{equation}
\begin{split}
&\left\|\nabla A(\Theta_t^\top x_i)-\nabla A(\Theta^\top x_i)\right\|_2\\
\leq&\sup_{\eta\in \bbR^K}\left[\sum_{i=1}^K\|(\nabla^2 A(\eta))_{i,:}\|^2_1\right]^{\frac{1}{2}}\|(\Theta_t-\Theta)^\top x_i\|_{\infty}\\
\leq &2\|(\Theta_t-\Theta)^\top x_i\|_{\infty}\\
\leq &2C_xR.
\end{split}
\end{equation}

\begin{equation}
\mathrm{\RNum{1}}(\Delta) \geq \frac{h^{\MN}(R_0)}{n}\sum_{i=1}^n \|\Delta^\top x_i\|_2^2,
\end{equation}

where we have applied the fact that $R\leq R_0$ and $h^{\MN}(\cdot)$ is a decreasing function.

\begin{equation}
A_1=\bbB_{\omega,2,1}(\sqrt{\rho})\cap \bbB_{2}(1)
\end{equation}

\begin{equation}
A_2=\left(1+\frac{2}{\min_j \omega_j}\right)\mathrm{cl}\{\mathrm{conv}\{\bbB_{\cG,0}(\lfloor \rho\rfloor)\cap \bbB_2(1)\}\}.
\end{equation}

Since both $A_1$ and $A_2$ are convex sets, it suffices to show that support function $\phi_{A_1}(U)\leq \phi_{A_2}(U)$ for any $U\in \bbR_{p\times K}$, where $\phi_{A_1}(U)=\sup_{V\in A_1}\langle V,U\rangle$ and $\phi_{A_2}(U)=\sup_{V\in A_2}\langle V,U\rangle$.
Let $\cG_{\lfloor \rho\rfloor}$ be the union of $\lfloor \rho\rfloor$ $\cG_j$'s such that $\|U_{\cR_j,\cC_j}\|_2$ are the largest. Then

\begin{equation}
\begin{split}
\phi_{A_1}(U)\leq &\|U_{\cG_{\lfloor \rho\rfloor}}\|_2+\sqrt{\rho}\|U_{\cG^c_{\lfloor \rho\rfloor}}\|_{\omega^{-1},2,\infty}\\
\leq&\|U_{\cG_{\lfloor \rho\rfloor}}\|_2+\frac{\sqrt{\rho}}{\sqrt{\lfloor \rho\rfloor}\min_j\omega_j}\|U_{\cG_{\lfloor \rho\rfloor}}\|_2\\
\leq&\left(1+\frac{2}{\min_j\omega_j}\right)\|U_{\cG_{\lfloor \rho\rfloor}}\|_2.
\end{split}
\end{equation}

On the other hand, one can also show that

\begin{equation}
\begin{split}
\phi_{A_2}(U)=&\left(1+\frac{2}{\min_j\omega_j}\right)\sup_{V\in \bbB_{\cG,0}(\lfloor \rho\rfloor)\cap\bbB_2(1)}\langle V,U\rangle\\
=&\left(1+\frac{2}{\min_j\omega_j}\right)\|U_{\cG_{\lfloor \rho\rfloor}}\|_2.
\end{split}
\end{equation}

Therefore, we have $A_1\subset A_2$.

By assumption <ref>, $\{(x_i\Delta_{:,j}\}_{i=1}^n$ are independent sub-Gaussian random variables with parameter $\sigma\|\Delta_{:,j}\|_2$, which implies that $\|(x_i^\top \Delta_{:,j})^2\|_{\psi_1}\leq C\sigma^2\|\Delta_{:,j}\|_2^2$. Thus

\begin{equation}\label{eq:rec_dev_col}
\bbP\left(\left|\Delta_{:,j}^\top \left(\frac{1}{n}X^\top X-\Sigma\right)\Delta_{:,j}\right|\geq t\|\Delta_{:,j}\|_2^2\right)\leq 2\exp\left\{-cn\min\left\{\frac{t^2}{\sigma^4},\frac{t}{\sigma^2}\right\}\right\}.
\end{equation}

Since $g(\Delta,\Delta)=\sum_{j=1}^K\Delta_{:,j}^\top \left(\frac{1}{n}X^\top X-\Sigma\right)\Delta_{:,j}$, (<ref>) implies that

\begin{equation}
\bbP\left(\left|g(\Delta,\Delta)\right|\geq t\right)\leq 2K\exp\left\{-cn\min\left\{\frac{t^2}{\sigma^4},\frac{t}{\sigma^2}\right\}\right\}.
\end{equation}

Let $\{\widetilde{U}_i\}_{i=1}^n$ be independent copies of $\{U_i\}_{i=1}^n$, and $\{\varepsilon_i\}_{i=1}^n$ be independent Rademacher random variables. Then one can show that

\begin{equation*}
\begin{split}
\leq&\bbE_{U,\widetilde{U}}\left(\sup_{\gamma\in\Gamma}\sum_{i=1}^n(\gamma(U_i)-\gamma(\widetilde{U}_i))\right)\\
\end{split}
\end{equation*}

where the 4th line is due to that

\begin{equation*}
\sum_{i=1}^n(\gamma(U_i)-\gamma(\widetilde{U}_i))\leq \sup_{\gamma\in\Gamma}\sum_{i=1}^n(\gamma(U_i)-\gamma(\widetilde{U}_i))
\end{equation*}
for any $\gamma\in \Gamma$, and the last line is due to that $U_i,\widetilde{U}_i, i=1,\dots,n$ are all independent, $U_i$ and $\widetilde{U}_i$ are identically distributed. Furthermore, we have

\begin{equation*}
\begin{split}
\leq&\bbE\left(\sup_{\gamma\in\Gamma}\sum_{i=1}^n\varepsilon_i\gamma(U_i)\right)+\bbE\left(\sup_{\gamma\in\Gamma}\sum_{i=1}^n(-\varepsilon_i)\gamma(\widetilde{U}_i)\right)\\
\end{split}
\end{equation*}

where we have utilized the symmetricity of $\varepsilon_i$.

First note that

\begin{equation}
\begin{split}
\nabla h_i(u)=&\left[\nabla A(\Theta^{*\top}x_i)-\nabla A(\Theta^{*\top}x_i+u)\right]1_K^\top \epsilon_i\\
&-\nabla^2 A(\Theta^{*\top}x_i+u)^\top u1_K^\top \epsilon_i,
\end{split}
\end{equation}

and $|1_K^\top \epsilon_i|=|1_K^\top \delta_i-1_K^\top\bbE\delta_i|\leq 1$. Recall the definition of $A(u)=\log(1+\sum_{k=1}^Ke^{u_i})$, one can show that $\|\nabla A(u)\|_2\leq 1$. As shown in the proof of Lemma <ref>, $\nabla^2 A(u)$ is positive definite with each entry

\begin{equation*}
(\nabla^2 A(u))_{jk}=\frac{e^{u_j}}{1+\sum_{l=1}^Ke^{u_l}}\left(\ind{j=k}-\frac{e^{u_k}}{1+\sum_{l=1}^Ke^{u_l}}\right).
\end{equation*}

\begin{equation*}
\begin{split}
\lambda_{\max}(\nabla^2 A(u))=&\sup_{\|v\|_2=1}v^\top \nabla^2A(u)v\\
\leq&\sup_{\|v\|_2=1}\sum_{j=1}^K\frac{v_j^2e^{u_j}}{1+\sum_{l=1}^Ke^{u_l}}\\
\leq&1,
\end{split}
\end{equation*}

which implies $\|\nabla^2A(u)v\|_2\leq \|v\|_2$ for any $u,v\in \bbR^K$. Meanwhile, for any $u\in \bbR^K$,

\begin{equation*}
\begin{split}
\|\nabla^2A(u)u\|_2^2=&\sum_{i=1}^K\left(\frac{u_ie^{u_i}}{1+\sum_{k=1}^Ke^{u_k}}-\frac{\sum_{j=1}^Ku_j e^{u_i+u_j}}{(1+\sum_{k=1}^Ke^{u_k})^2}\right)^2\\
=&\frac{\sum_{i=1}^K \left(u_ie^{u_i}+\sum_{j=1}^K(u_i-u_j)e^{u_i+u_j}\right)^2}{(1+\sum_{k=1}^Ke^{u_k})^4}\\
\leq&\frac{K^2\sum_{i=1}^K e^{4u_i}}{(1+\sum_{k=1}^Ke^{u_k})^4}\\
\leq &K^2,
\end{split}
\end{equation*}

where we have applied the fact that $x\leq e^x$ for any $x\in \bbR$ on the third line. Therefore,

\begin{equation}
\begin{split}
\|\nabla h_i(u)\|_2\leq& 2+\|\nabla^2A(\Theta^{*\top} x_i+u)^\top u\|_2\\
\leq& 2+\|\nabla^2A(\Theta^{*\top} x_i+u)^\top \Theta^*x_i\|_2\\
&+\|\nabla^2A(\Theta^{*\top} x_i+u)^\top (\Theta^*x_i+u)\|_2\\
\leq& 2+\|\Theta^*x_i\|_2+K\\
\leq &K+\sqrt{K}R^*C_x+2,
\end{split}
\end{equation}

where we used the fact that $\|\Theta^*\|_{1,\infty}\leq R^*$ and $\|x_i\|_{\infty}\leq C_x$ on the last line.

Let $E\in \bbR^{n\times K}$ such that $E_{ij}=\widetilde{\varepsilon}_{ij}$, then

\begin{equation*}
\bbE\left\|\sum_{i=1}^nx_i\widetilde{\varepsilon}_i^\top\right\|_{\omega^{-1},2,\infty}=\bbE\max_{1\leq j\leq J}\frac{1}{\omega_j}\|X_{\cR_j}^\top E_{\cC_j}\|_2.
\end{equation*}

First note that by similar arguments from the proof of Lemma <ref>, one can show that

\begin{equation}
\bbE\left[\|X_{\cR_j}^\top E_{\cC_j}\|_2|X\right]\leq \left(c_j\tr(X_{\cR_j}X_{\cR_j}^\top)\right)^{\frac{1}{2}}\leq \sqrt{r_jc_j}\|X_{\cR_j}\|.
\end{equation}

Meanwhile, applying Talagrand's contraction inequality (see e.g., Theorem 5.2.16 in [Vershynin, 2018]) shows us that

\begin{equation}\label{eq:RSC_term2_tail1}
\bbP\left(\|X_{\cR_j}^\top E_{\cC_j}\|_2>C\|X_{\cR_j}\|(\sqrt{r_jc_j}+t_j)|X\right)\leq \exp\{-ct_j^2\}.
\end{equation}

To obtain a tail probability bound without conditioning on $X$, we apply bounds on spectral norm of matrices with independent isotropic sub-Gaussian rows to $X_{\cR_j}\Sigma_{\cR_j,\cR_j}^{-\frac{1}{2}}$, then for any $t_j>0$,

\begin{equation}\label{eq:RSC_term2_tail2}
\bbP\left(\|X_{\cR_j}\|> \lambda^{\frac{1}{2}}_{\max}(\Sigma)(C\sqrt{n}+t_j)\right)\leq 2\exp\{-ct_j^2\}.
\end{equation}

Combining (<ref>) and (<ref>), one can show that

\begin{equation}
\begin{split}
&\bbP\left(\|X_{\cR_j}^\top E_{\cC_j}\|_2>C\lambda^{\frac{1}{2}}_{\max}(\Sigma)(\sqrt{r_jc_jn}+(\sqrt{n}+\sqrt{r_jc_j})t_j+t_j^2)\right)\\
\leq &3\exp\{-nt_j^2\},
\end{split}
\end{equation}

which is equivalent to

\begin{equation}
\begin{split}
&\bbP\left(\|X_{\cR_j}^\top E_{\cC_j}\|_2>C\lambda^{\frac{1}{2}}_{\max}(\Sigma)\sqrt{r_jc_jn}+t_j\right)\\
\leq &3\exp\left\{-\min\left\{\frac{t_j^2}{\lambda_{\max}(\Sigma)(n+r_jc_j)},\frac{t_j}{\lambda^{\frac{1}{2}}_{\max}(\Sigma)}\right\}\right\}.
\end{split}
\end{equation}

Therefore, we can bound $\bbE\max_{1\leq j\leq J}\frac{1}{\omega_j}\|X_{\cR_j}^\top E_{\cC_j}\|_2$ as follows:

\begin{equation}
\begin{split}
&\bbE\max_{1\leq j\leq J}\frac{1}{\omega_j}\|X_{\cR_j}^\top E_{\cC_j}\|_2\\
\leq&\int_{0}^{\infty}\bbP\left(\max_{1\leq j\leq J}\frac{1}{\omega_j}\|X_{\cR_j}^\top E_{\cC_j}\|_2>t\right)\mathrm{d}t\\
\leq&C\lambda^{\frac{1}{2}}_{\max}(\Sigma)\sqrt{mn}+C\lambda^{\frac{1}{2}}_{\max}(\Sigma)\sqrt{\log J(m+n)}\\
&+\int_{C\lambda^{\frac{1}{2}}_{\max}(\Sigma)\sqrt{\log J(m+n)}}^{\infty}3J\exp\left\{-c\min\left\{\frac{t}{\lambda^{\frac{1}{2}}_{\max}(\Sigma)},\frac{t^2}{\lambda_{\max}(\Sigma)(m+n)}\right\}\right\}\mathrm{d}t\\
\leq&C\lambda^{\frac{1}{2}}_{\max}(\Sigma)\sqrt{mn}+C\lambda^{\frac{1}{2}}_{\max}(\Sigma)\sqrt{\log J(m+n)}\\
\leq&C\lambda^{\frac{1}{2}}_{\max}(\Sigma)(\sqrt{mn}+\sqrt{(\log J(m+n)}+\sqrt{m+n})\\
\leq&C\lambda^{\frac{1}{2}}_{\max}(\Sigma)\sqrt{(m+\log J)n},
\end{split}
\end{equation}

where we have applied Assumption <ref> on the last line.

§.§ Supporting Lemmas for the Ordinal-PU Model

Let $u_i=u(x_i,\theta)$, $u_i^*=u(x_i,\theta^*)$. First note that $\exists t\in [0,1]$ such that

\begin{equation*}
\begin{split}
=&\frac{1}{n}\sum_{i=1}^n \begin{pmatrix}x_i^\top \Delta_{1:p}-\Delta_{p+1}\\
\vdots\\
-\Delta_{p+K}\end{pmatrix}^\top\nabla f^{\ON}(u_i^t)^\top\nabla^2 A(f^{\ON}(u_i^t))\nabla f^{\ON}(u_i) \begin{pmatrix}x_i^\top \Delta_{1:p}-\Delta_{p+1}\\
\vdots\\
\end{split}
\end{equation*}

where $u_i^t=u_i^*+t(u_i-u_i^*)$. Then some calculations show that

\begin{equation*}
\begin{split}
\mathrm{\RNum{1}}(\Delta)\geq&\frac{1}{n}\sum_{i=1}^n\Bigg[\lambda_{\min}\left(\nabla f^{\ON}(u_i^t)^\top \nabla^2 A(f^{\ON}(u_i^t))\nabla f^{\ON}(u_i^t)\right)\\
&-\left\|\nabla f^{\ON}(u_i^t)^\top\nabla^2 A(f^{\ON}(u_i^t)) (\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i))\right\|\Bigg]\\
\end{split}
\end{equation*}

To further lower bound the term above, two key steps here are to (i) lower bound the minimum eigenvalue of $\nabla f^{\ON}(u_i^t)^\top \nabla f^{\ON}(u_i^t)$ and to (ii) upper bound $\|\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i)\|$.

The second key step is why we need to constrain the distance between $\theta$ and $\theta^*$ in Lemma <ref>: by focusing on $\theta\in\theta^*+S(R_0,r_0)$ for some $R_0, r_0>0$, we can upper bound $\|\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i)\|$ appropriately so that we can show that $\mathrm{\RNum{1}}(\Delta)\geq\frac{\gamma}{n}\sum_{i=1}^n\left((x_i^\top\Delta_{1:p}-\Delta_{p+1})^2+\sum_{j=2}^{K}\Delta_{p+j}^2\right)$ for some $\gamma>0$.
The proof of Lemma <ref> and Lemma <ref> are included in Section <ref>.
As shown in the proof of Lemma <ref>, for any vector $\eta\in \bbR^K$,

\begin{equation*}
\lambda_{\min}(\nabla^2 A(\eta))\geq \frac{\min_j e^{\eta_j}}{\left(1+\sum_{l=1}^Ke^{\eta_l}\right)^2}.
\end{equation*}

Recall that $\Delta\in S(R_0,r_0)$, which implies that $\max\{\|\Delta_{1:p}\|_1,\|\Delta_{(p+1):(p+K)}\|_1\}\leq R_0$, and $-\min_{2\leq j\leq K}\Delta_{p+j}\leq r_0<r^*$. Since $\exp\{f_j(u_i^t)\}=\frac{n_j}{\pi_jn_u}p_j(u_i^t)$, where $p_j$ is as defined in (<ref>), then by Lemma <ref>, one can show that

\begin{equation*}
\min_j\exp\{f_j(u_i^t)\}=\min_j\frac{n_j}{\pi_jn_u}p_j(u_i^t)\geq \min_j\frac{n_j}{\pi_jn_u}(h^{\ON}(R_0,r_0))^{-1}.
\end{equation*}

Thus we have

\begin{equation*}
\lambda_{\min}(\nabla^2 A(f^{\ON}(u_i^t)))\geq \min_j\frac{n_j}{4\pi_jn_uh^{\ON}(R_0,r_0)}.
\end{equation*}

Now we recall Lemma <ref> and Lemma <ref> in the main paper.
Here we give a more slight modification of Lemma <ref> that better suits our purposes:

For any $\Delta\in S(R_0,r_0)$,

\|\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i)\|\leq \sqrt{5}[(h^{\ON}(R_0,r_0))^2+h^{\ON}(R_0,r_0)](C_x+1)R_0K,

where $h^{\ON}(\cdot)$ is as defined in (<ref>).

By Lemma <ref>,

\begin{equation}
\begin{split}
&\lambda_{\min}(\nabla f^{\ON}(u_i^t)^\top \nabla^2 A(f^{\ON}(u_i^t))\nabla f^{\ON}(u_i^t))\\
\geq&\lambda_{\min}(\nabla f^{\ON}(u_i^t)^\top\nabla f^{\ON}(u_i^t))\lambda_{\min}(\nabla^2 A(f^{\ON}(u_i^t)))\\
\geq &\min_j\frac{n_j}{\pi_jn_u}\frac{e^{2(C_x+1)(R^*+R_0)}}{16Kh^{\ON}(R_0,r_0)(1+e^{(C_x+1)(R^*+R_0)})^4}.
\end{split}
\end{equation}

Now we upper bound $\left\|\nabla f^{\ON}(u_i^t)\nabla^2 A(f^{\ON}(u_i^t)) (\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i))\right\|$. By (<ref>),

\begin{equation*}
\begin{split}
\|\nabla f^{\ON}(u_i^t)v\|_2^2=&\sum_{j=1}^{K-1}\frac{(\alpha_j(u_i^t)\widetilde{v}_j-\alpha_{j+1}(u_i^t)\widetilde{v}_{j+1})^2}{p_j^2(u_i^t)}+\frac{\alpha_K^2(u_i^t)\widetilde{v}_K^2}{p_K^2(u_i^t)}\\
\leq&\sum_{j=1}^{K-1}\frac{2\alpha_j^2(u_i^t)\widetilde{v}^2_j+2\alpha_{j+1}^2(u_i^t)\widetilde{v}^2_{j+1}}{p_j^2(u_i^t)}+\frac{\alpha_K^2(u_i^t)\widetilde{v}_K^2}{p_K^2(u_i^t)}\\
\leq&\frac{2(1+e^{(C_x+1)(R^*+R_0)})^4}{(r^*-r_0)^2e^{2(C_x+1)(R^*+R_0)}}\sum_{j=1}^{K-1}(\widetilde{v}^2_j+\widetilde{v}^2_{j+1})\\
\leq&(h^{\ON}(R_0,r_0))^2\left(2\sum_{j=1}^{K-1}(\widetilde{v}_j^2+\widetilde{v}_{j+1}^2)+\widetilde{v}_K^2\right)\\
\leq&(h^{\ON}(R_0,r_0))^2\left(2\sum_{j=1}^{K-1}(j\sum_{l=1}^jv_l^2+(j+1)\sum_{l=1}^{j+1}v_{l}^2)+K\sum_{l=1}^Kv_l^2\right)\\
\leq&(h^{\ON}(R_0,r_0))^2\sum_{l=1}^K\left(2l+2\sum_{j=l}^{K-1}(2j+1)+K\right)v_l^2\\
\leq&3(h^{\ON}(R_0,r_0))^2K^2\|v\|_2^2,
\end{split}
\end{equation*}

which implies that

\begin{equation*}
\|\nabla f^{\ON}(u_i^t)\|\leq \sqrt{3}h^{\ON}(R_0,r_0)K.
\end{equation*}

Furthermore, by Lemma <ref>,

\begin{equation*}
\begin{split}
\|\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i)\|\leq \sqrt{5}((h^{\ON}(R_0,r_0))^2+h^{\ON}(R_0,r_0))(C_x+1)R_0K.
\end{split}
\end{equation*}

While for $\|\nabla^2 A(f^{\ON}(u_i^t))\|$, note that for any $v\in \bbR^K$,

\begin{equation*}
\begin{split}
v^\top \nabla^2 A(f^{\ON}(u_i^t))v=&\sum_{j=1}^K(\nabla A(f^{\ON}(u_i^t)))_jv_j^2-\left(\sum_{j=1}^K(\nabla A(f^{\ON}(u_i^t)))_jv_j\right)^2\\
\leq&\sum_{j=1}^K(\nabla A(f^{\ON}(u_i^t)))_jv_j^2\\
\leq&\max_j\frac{e^{f^{\ON}(u_{ij}^t)}}{1+\sum_{k=1}^Ke^{f^{\ON}(u_{ik}^t)}}\|v\|_2^2\\
\leq&\|v\|_2^2,
\end{split}
\end{equation*}

and hence $\|\nabla^2 A(f^{\ON}(u_i^t))\|=\lambda_{\max}(\nabla^2 A(f^{\ON}(u_i^t)))\leq 1$.

\begin{equation*}
\begin{split}
&\|\nabla f^{\ON}(u_i^t)^\top\nabla^2A(f^{\ON}(u_i^t))(\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i))\|\\
\leq&4((h^{\ON}(R_0,r_0))^3+(h^{\ON}(R_0,r_0))^2)(C_x+1)R_0K^2.
\end{split}
\end{equation*}

By Lemma <ref>, we know that

\begin{equation*}
R_0\leq \min_j\frac{n_j}{\pi_jn_u}\frac{\left(1+e^{(C_x+1)(R^*+R_0)}\right)^{-2}}{512(C_x+1)K^3(h^{\ON}(R_0,r_0))^3(1+h^{\ON}(R_0,r_0))},
\end{equation*}

which then implies

\begin{equation*}
\begin{split}
&\|\nabla f^{\ON}(u_i^t)^\top\nabla^2A(f^{\ON}(u_i^t))(\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i))\|\\
\leq&\min_j\frac{n_j}{\pi_jn_u}\frac{\left(1+e^{(C_x+1)(R^*+R_0)}\right)^{-2}}{128Kh^{\ON}(R_0,r_0)}\\
\leq & \min_j\frac{n_j}{\pi_jn_u}\frac{e^{2(C_x+1)(R^*+R_0)}}{32Kh^{\ON}(R_0,r_0)\left(1+e^{(C_x+1)(R^*+R_0)}\right)^{4}},
\end{split}
\end{equation*}

where the last line is due to the fact that $e^{(C_x+1)(R^*+R_0)}\geq 1$.

\begin{equation*}
\begin{split}
&\lambda_{\min}(\nabla f^{\ON}(u_i^t)^\top \nabla^2A(f^{\ON}(u_i^t))\nabla f^{\ON}(u_i^t))\\
&-\|\nabla f^{\ON}(u_i^t)^\top \nabla^2A(f^{\ON}(u_i^t))(\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i))\|\\
\geq&\min_j\frac{n_j}{\pi_jn_u}\frac{e^{2(C_x+1)(R^*+R_0)}}{32Kh^{\ON}(R_0,r_0)(1+e^{(C_x+1)(R^*+R_0)})^{4}},
\end{split}
\end{equation*}

which finishes the proof.

By assumption <ref>, $\{(x_i^\top \widetilde{\Delta}_{1:p}\}_{i=1}^n$ are independent sub-Gaussian random variables with parameter $\sigma\|\widetilde{\Delta}_{1:p}\|_2$, which implies that

\begin{equation*}
\|(x_i^\top \widetilde{\Delta}_{1:p})^2\|_{\psi_1}\leq C\sigma^2\|\widetilde{\Delta}_{1:p}\|_2^2\leq C\sigma^2,
\end{equation*}

\begin{equation*}
\|x_i^\top \widetilde{\Delta}_{1:p}|\widetilde{\Delta}_{p+1}|\|_{\psi_2}\leq \sigma|\widetilde{\Delta}_{p+1}|\|\widetilde{\Delta}_{1:p}\|_2\leq \frac{\sigma}{2}.
\end{equation*}

Meanwhile, note that

\begin{equation*}
g(\widetilde{\Delta},\widetilde{\Delta})=\frac{1}{n}\sum_{i=1}^n(x_i^\top\widetilde{\Delta}_{1:p})^2-\widetilde{\Delta}_{1:p}^\top \Sigma\widetilde{\Delta}_{1:p}-2\widetilde{\Delta}_{p+1}\widetilde{\Delta}_{1:p}^\top x_i
\end{equation*}

\begin{equation}
\begin{split}
\bbP\left(|g(\widetilde{\Delta},\widetilde{\Delta})|\geq t\right)\leq&\bbP\left(\left|\frac{1}{n}\sum_{i=1}^n(x_i^\top\widetilde{\Delta}_{1:p})^2-\widetilde{\Delta}_{1:p}^\top \Sigma\widetilde{\Delta}_{1:p}\right|\geq \frac{t}{2}\right)\\
&+\bbP\left(\left|\widetilde{\Delta}_{p+1}\widetilde{\Delta}_{1:p}^\top x_i\right|\geq \frac{t}{4}\right)\\
\leq&\exp\left\{-cn\min\left\{\frac{t^2}{\sigma^4},\frac{t}{\sigma^2}\right\}\right\}+\exp\left\{-cn\frac{t^2}{\sigma^2}\right\}\\
\leq&2\exp\left\{-cn\min\left\{\frac{t^2}{\sigma^4},\frac{t^2}{\sigma^2},\frac{t}{\sigma^2}\right\}\right\}.
\end{split}
\end{equation}

For any $1\leq l\leq K$, some calculation shows that

\begin{equation*}
\begin{split}
\left|(\nabla q_i(u_i(\Delta)))_l\right|=&\bigg|\epsilon_i^\top\left(\nabla f^{\ON}(u_i(\theta^*)+u_i(\Delta))-\nabla f^{\ON}(u_i(\theta^*))\right)_{:,l}\\
&+\epsilon_i^\top (\nabla^2f^{\ON}(u_i(\theta^*)+u_i(\Delta)))_{:,:,l}u_i(\Delta)\bigg|\\
\leq&2\left\|\nabla f^{\ON}(u_i(\theta^*)+u_i(\Delta))-\nabla f^{\ON}(u_i(\theta^*))\right\|_{\infty}\\
\end{split}
\end{equation*}

By Lemma <ref>, we have

\begin{equation*}
\left|(\nabla q_i(u_i(\Delta)))_l\right|\leq h^{\ON}(R_0,r_0)+2\|u_i(\Delta)\|_1\left[\frac{2}{3}h^{\ON}(R_0,r_0)+\frac{1}{8}(h^{\ON}(R_0,r_0))^2\right].
\end{equation*}

Meanwhile, since
$\|u_i(\Delta)\|_1\leq C_x\|\Delta_{1:p}\|_1+\|\Delta_{(p+1):(p+K)}\|_1\leq (C_x+1)R_0$,
we obtain that

\begin{equation*}
\begin{split}
\|\nabla q_i(u_i(\Delta))\|_2\leq &\sqrt{K}h^{\ON}(R_0,r_0)+2\sqrt{K}(C_x+1)R_0\left[\frac{2}{3}h^{\ON}(R_0,r_0)+\frac{1}{8}(h^{\ON}(R_0,r_0))^{2}\right]\\
\end{split}
\end{equation*}

Firstly, we note that Assumption <ref> suggests that

\begin{equation*}
\begin{split}
R_0&\leq 80\frac{\left(1+e^{(C_x+1)(R^*+0.01)}\right)^{-6}}{512(C_x+1)K^3(1+\frac{2}{r^* - r_0})^3(2+\frac{2}{r^*-r_0})}\\
&\leq \frac{80}{1024K^3}\\
&< 0.01,
\end{split}
\end{equation*}

where we have applied $K\geq 2$ in the last line. Hence Assumption <ref> also implies

\begin{equation*}
\begin{split}
R_0\leq \min_j\frac{n_j}{\pi_jn_u}\frac{\left(1+e^{(C_x+1)(R^*+R_0)}\right)^{-6}}{512(C_x+1)K^3(1+\frac{2}{r^* - r_0})^3(2+\frac{2}{r^*-r_0})}.
\end{split}
\end{equation*}

On the other hand, we know that $h^{\ON}(R_0,r_0) \leq (1+e^{(C_x+1)(R^*+R_0)})(1+\frac{2}{r^*-r_0})$. This bound then implies the second inequality in Lemma <ref>.

§.§ Proof of Properties of $f^{\ON}$ for the Ordinal-PU model

In this section, we prove Lemma <ref> and Lemma <ref>, which tackles some of the major challenges brought by the complicated form of ordinal log-likelihood losses.
Some calculations show that, for any vector $u\in\bbR^K$, $\nabla f^{\ON}(u)\in \bbR^{K\times K}$ satisfies

\begin{equation*}
\nabla f^{\ON}(u)_{jk}=\frac{\alpha_j(u)\ind{k\leq j}-\alpha_{j+1}(u)\ind{k\leq j+1}}{p_j(u)}.
\end{equation*}

where functions $p:\bbR^K\rightarrow \bbR^K$ and $\alpha:\bbR^K\rightarrow \bbR^{K+1}$ are defined as follows:

\begin{equation}\label{eq:ordinal_alpha_def}
\alpha_j(u)=\begin{cases}e^{\sum_{l=1}^ju_{l}}\left(1+e^{\sum_{l=1}^{j}u_{l}}\right)^{-2},&j\leq K,\\
\end{cases}
\end{equation}

\begin{equation}\label{eq:ordinal_p_def}
(1+e^{\sum_{l=1}^{j+1}u_{l}})^{-1}-(1+e^{\sum_{l=1}^ju_{l}})^{-1},&1\leq j<K,\\
\end{cases}
\end{equation}

where $\alpha_j(u),p_j(u)$ refers to the $j$th coordinate of $\alpha(u)$ and $p(u)$.

By the definition of function $f^{\ON}$ in (<ref>), one can show that for any $1\leq j,k\leq K$,

\begin{equation}\label{eq:ordinal_f_grad}
(\nabla f^{\ON}(u_i^*))_{jk}=p_j(u_i^*)^{-1}\left[\alpha_j(u_i^*)\ind{k\leq j}-\alpha_{j+1}(u_i^*)\ind{k\leq j+1}\right].
\end{equation}

For any vector $v\in \bbR^K$, we have

\begin{equation*}
\begin{split}
\nabla f^{\ON}(u_i^t)v=&\begin{pmatrix}
\frac{\alpha_1(u_i^t)}{p_1(u_i^t)}&-\frac{\alpha_2(u_i^t)}{p_1(u_i^t)}&0&\hdots&0\\
\vdots&\vdots&\ddots&\ddots&0\\
\end{pmatrix}\begin{pmatrix}v_1\\v_1+v_2\\\sum_{j=1}^3v_j\\\vdots\\\sum_{j=1}^Kv_j\end{pmatrix}\\
\end{split}
\end{equation*}

where $\widetilde{v}_k=\sum_{j=1}^kv_j$. Then one can show that

\begin{equation}\label{eq:ordinal_RSC_f_grad_lwrbnd}
\begin{split}
\left\|\nabla f^{\ON}(u_i^t)v\right\|_2^2=&\sum_{j=1}^{K-1}\frac{(\alpha_j(u_i^t)\widetilde{v}_j-\alpha_{j+1}(u_i^t)\widetilde{v}_{j+1})^2}{p_j^2(u_i^t)}+\frac{\alpha_K^2(u_i^t)\widetilde{v}_K^2}{p_K^2(u_i^t)}\\
\geq&\frac{1}{K}\sum_{k=1}^{K-1}\sum_{j=k}^{K-1}\frac{(\alpha_j(u_i^t)\widetilde{v}_j-\alpha_{j+1}(u_i^t)\widetilde{v}_{j+1})^2}{p_j^2(u_i^t)}+\frac{\alpha_K^2(u_i^t)\widetilde{v}_K^2}{p_K^2(u_i^t)},
\end{split}
\end{equation}

By the Cauchey-Schwarz inequality,

\begin{equation*}
\begin{split}
\geq &\left(\sum_{j=k}^{K-1}(\alpha_j(u_i^t)\widetilde{v}_j-\alpha_{j+1}(u_i^t)\widetilde{v}_{j+1})+\alpha_K(u_i^t)\widetilde{v}_K\right)^2 =  \alpha^2_k(u^t_i)\widetilde{v}^2_k,
\end{split}
\end{equation*}

which implies

\begin{equation*}
\begin{split}
\left\|\nabla f^{\ON}(u_i^t)v\right\|_2^2\geq&\frac{1}{K}\Bigg[\frac{\alpha_K^2(u_i^t)\widetilde{v}_K^2}{p_K^2(u_i^t)}+\sum_{k=1}^{K-1}\left(\sum_{j=k}^{K}p_j^2(u_i^t)\right)^{-1}\alpha^2_k(u^t_i)\widetilde{v}^2_k\Bigg]\\
\geq&\frac{1}{K}\sum_{k=1}^K\alpha^2_k(u^t_i)\widetilde{v}^2_k.
\end{split}
\end{equation*}

In addition, since

\begin{equation*}
\begin{split}
\|v\|_2^2=\sum_{j=1}^{K-1}(\widetilde{v}_{j+1}-\widetilde{v}_j)^2+\widetilde{v}_1^2\leq 4\sum_{j=1}^K\widetilde{v}_j^2\leq \frac{4}{\min_{1\leq k\leq K}\alpha_k^2(u_i^t)}\sum_{j=1}^K\alpha_j^2(u_i^t)\widetilde{v}_j^2,
\end{split}
\end{equation*}

we have

\left\|\nabla f^{\ON}(u_i^t)v\right\|_2^2\geq \frac{\min_k\alpha_k^2(u_i^t)}{4K}\|v\|_2^2.

Recall the definition of $\alpha_k(u)$ in (<ref>) and the fact that

\begin{equation*}
\left|\sum_{j=1}^ku_{ij}^t\right|\leq (C_x+1)(R^*+R_0),
\end{equation*}

one can show that for $1\leq k\leq K$,

\begin{equation*}
\alpha_k(u_i^t)\geq \frac{e^{(C_x+1)(R^*+R_0)}}{(1+e^{(C_x+1)(R^*+R_0)})^2}.
\end{equation*}

\begin{equation*}
\lambda_{\min}(\nabla f^{\ON}(u_i^t)^\top \nabla f^{\ON}(u_i^t))\geq \frac{e^{2(C_x+1)(R^*+R_0)}}{4K(1+e^{(C_x+1)(R^*+R_0)})^4}.
\end{equation*}

While for bounding $\|\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i)\|$, we first recall each entry of $\nabla f^{\ON}$ in (<ref>), and note that

\begin{equation}\label{eq:ordinal_rsc_f_grad_diff}
\begin{split}
&\|\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i)\|\\
\leq &\|\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i)\|_F\\
\leq&\sqrt{2K^2+3K-1}\max_{j,k}\left|\frac{\alpha_k(u_i^t)}{p_j(u_i^t)}-\frac{\alpha_k(u_i)}{p_j(u_i)}\right|.
\end{split}
\end{equation}

One can show that

\begin{equation}\label{eq:ordinal_rsc_f_grad_diff0}
\begin{split}
\leq&\max_j\left|p_j^{-1}(u_i^t)-p_j^{-1}(u_i)\right|+\max_{j,k}p_j^{-1}(u_i)\left|\alpha_k(u_i^t)-\alpha_k(u_i)\right|,
\end{split}
\end{equation}

where we have applied the fact that $0<\alpha_k(w)< 1$ for any $w\in \bbR^K$ on the second line. Meanwhile, for any $w\in \bbR^K$, $1\leq j<K$,

\begin{equation*}
\nabla_{w} (p_j^{-1}(w))=\frac{e^{\sum_{l=1}^{j+1}w_l}}{p_j^{2}(w)(1+e^{\sum_{l=1}^{j+1}w_l})^{2}}(1_{j+1}^\top,0_{K-j-1}^\top)^\top-\frac{e^{\sum_{l=1}^{j}w_l}}{p_j^{2}(w)(1+e^{\sum_{l=1}^{j}w_l})^{2}}(1_{j}^\top,0_{K-j}^\top)^\top,
\end{equation*}

\begin{equation*}
\nabla_w (p_K^{-1}(w))=-p_K^{-2}(w)\frac{e^{\sum_{l=1}^Kw_l}}{(1+e^{\sum_{l=1}^K w_l})^{2}}1_K
\end{equation*}

Since $0<\frac{e^x}{(1+e^x)^2}<1$ for $x\in \bbR$, we have

\begin{equation*}
\left\|\nabla_w (p_j^{-1}(w))\right\|_{\infty}\leq p_j^{-2}(w)\leq (h^{\ON}(R,r))^2,
\end{equation*}

for any $w$ lying between $u_i^t$ and $u_i$, which implies that

\begin{equation}\label{eq:ordinal_rsc_f_grad_diff1}
\max_j|p_j^{-1}(u_i^t)-p_j^{-1}(u_i)|\leq (h^{\ON}(R,r))^2\|u_i^t-u_i\|_1.
\end{equation}

While for bounding $p_j^{-1}(u_i)|\alpha_k(u_i^t)-\alpha_k(u_i)|$, we apply the following lemma which gives several bounds for functions $p(\cdot)$, $\alpha(\cdot)$, $f^{\ON}(\cdot)$:

The functions $p(\cdot),\alpha(\cdot)$ and $f^{\ON}(\cdot)$ defined in (<ref>), (<ref>) and (<ref>) satisfy the following bounds:

\begin{equation}
\begin{split}
0<\alpha_j(u)\leq \frac{1}{4},\quad p_j(u)\geq\min\{\frac{-\max_{k>1}u_ke^{\|u\|_1}}{(1+e^{\|u\|_1})^2},\frac{1}{1+e^{\|u\|_1}}\},
\end{split}
\end{equation}

if $u_j<0$ for $j>1$.
In particular, if $u=u(x,\theta)$ for any $\|x\|_{\infty}\leq C_x$, and $\theta\in \theta^*+S(R,r)$ for
$$S(R,r)=\{\Delta:\|\Delta_{1:p}\|_1,\|\Delta_{(p+1):(p+K)}\|_1\leq R,\min_{2\leq j\leq K}\Delta_{p+j}\geq -r\},$$

\begin{equation}
\begin{split}
p_j(u)\geq (h^{\ON}(R,r))^{-1},\quad\|\nabla f^{\ON}(u)\|_{\infty}\leq \frac{1}{4}h^{\ON}(R,r)\\
\|\nabla^2f^{\ON}(u)\|_{\infty}\leq\frac{2}{3}h^{\ON}(R,r)+\frac{1}{8}(h^{\ON}(R,r))^2.
\end{split}
\end{equation}

Lemma <ref> is proved subsequently.
Hence by Lemma <ref>, $p_j^{-1}(u_i)\leq h^{\ON}(R_0,r_0)$. In addition, let $g(x)=\frac{e^x}{(1+e^x)^2}$, then $\alpha_k(u_i^t)-\alpha_k(u_i)=g(\sum_{j=1}^ku^t_{ij})-g(\sum_{j=1}^ku_{ij})$. Since

\begin{equation*}
|g'(x)|=\frac{|e^x-e^{2x}|}{(1+e^x)^3}\leq 1,
\end{equation*}

\begin{equation}\label{eq:ordinal_rsc_f_grad_diff2}
\begin{split}
\leq &h^{\ON}(R,r)\|u_i^t-u_i\|_1.
\end{split}
\end{equation}

Since $\|u_i^t-u_i\|_1\leq \|u_i^*-u_i\|_1\leq (C_x+1)R$, combining (<ref>), (<ref>), (<ref>) and (<ref>) leads us to

\begin{equation*}
\begin{split}
\|\nabla f^{\ON}(u_i^t)-\nabla f^{\ON}(u_i)\|\leq \sqrt{5}((h^{\ON}(R,r))^2+h^{\ON}(R,r))(C_x+1)RK,
\end{split}
\end{equation*}

Let $g(x)=\frac{x}{(1+x)^2}$ then one can show that $g(x)\leq g(1)=\frac{1}{4}$ for any $x\geq 0$, which directly implies that $0<\alpha_j(u)\leq\frac{1}{4}$.
While for $p_j(u)$, by (<ref>) one can show that

\begin{equation*}
\end{equation*}

For $1\leq j<K$, $p_j(u)=-u_{j+1}g(\xi)$ for some $\xi\in [e^{\sum_{l=1}^{j+1}u_l},e^{\sum_{l=1}^ju_l}]$. Noting the fact that $g(x)$ is increasing when $x<1$ and decreasing when $x>1$, and $g(\frac{1}{x})=g(x)$, we know that
$$g(\xi)\geq \min\{g(e^{\sum_{l=1}^{j+1}u_l}), g(e^{\sum_{l=1}^ju_l})\}\geq \frac{e^{\|u\|_1}}{(1+e^{\|u\|_1})^2},$$

which implies $p_j(u)\geq \frac{-\max_{k>1}u_ke^{\|u\|_1}}{(1+e^{\|u\|_1})^2}$. In addition, recall that we have defined $u$ as $u=u(x,\theta)$, and we assumed $\|x\|_{\infty}\leq C_x$, $\theta\in S(R,r)$. Thus one has

\begin{equation*}
\|u(x,\theta)\|_1\leq C_x\|\theta_{1:p}\|_1+\|\theta_{(p+1):(p+K)}\|_1\leq (C_x+1)(R^*+R),
\end{equation*}

and for $k>1$,

\begin{equation*}
u_k(x,\theta)\leq -\min_{k>1}\theta_{p+k}\leq r-r^*.
\end{equation*}

Thus by the definition of $h^{\ON}(\cdot,\cdot)$ in (<ref>),

\begin{equation*}
p_j(u)\geq (h^{\ON}(R,r))^{-1}
\end{equation*}

for $1\leq j\leq K$.

Now we start bounding $\|\nabla f^{\ON}(u)\|_{\infty}$ and $\|\nabla^2f^{\ON}(u)\|_{\infty}$. Some calculation shows that, for any vector $u\in\bbR^K$, $\nabla f^{\ON}(u)\in \bbR^{K\times K}$ satisfies

\begin{equation*}
\nabla f^{\ON}(u)_{jk}=\frac{\alpha_j(u)\ind{k\leq j}-\alpha_{j+1}(u)\ind{k\leq j+1}}{p_j(u)}.
\end{equation*}

Since $0<\alpha_j(u)\leq \frac{1}{4}$, $p_j(u)\geq (h^{\ON}(R,r))^{-1}$, we have $\|\nabla f^{\ON}(u)\|_{\infty}\leq \frac{1}{4}h^{\ON}(R,r)$. Meanwhile, one can show that

\begin{equation*}
\begin{split}
(\nabla^2 f^{\ON}(u))_{jkl}=&\left(\frac{(\nabla \alpha_j(u))_l}{p_j(u)}-\frac{\alpha_j(u)(\nabla p_j(u))_l}{p_j^2(u)}\right)\ind{k\leq j}\\
&-\left(\frac{(\nabla \alpha_{j+1}(u))_l}{p_j(u)}-\frac{\nabla p_j(u))_l\alpha_{j+1}(u)}{p_j^2(u)}\right)\ind{k\leq j+1},
\end{split}
\end{equation*}

and $|(\nabla p_j(u))_l|=|\alpha_j(u)\ind{l\leq j}-\alpha_{j+1}(u)\ind{l\leq j+1}|\leq \frac{1}{4}$,

\begin{equation*}
\begin{split}
|(\nabla \alpha_j(u))_l|=\frac{\left|e^{\sum_{m=1}^ju_j}-e^{2\sum_{m=1}^ju_j}\right|}{(1+e^{\sum_{m=1}^ju_j})^3}\ind{l\leq j}\leq\frac{1}{3}.
\end{split}
\end{equation*}

Thus we have

\begin{equation*}
\|\nabla^2f^{\ON}(u)\|_{\infty}\leq \max_j\left(\frac{2}{3}h^{\ON}(R,r)+\frac{1}{8}(h^{\ON}(R,r))^2\right).
\end{equation*}

§ REGULARIZED EM ALGORITHMS, ADDITIONAL CONVERGENCE GUARANTEES, AND PROOFS OF CONVERGENCE

EM algorithm for Multinomial-PU model: An alternative estimation method is the EM algorithm, which has been widely adopted in the literature on missing/hidden variable models. Here we present a regularized EM algorithm for our model.
Based on the full log-likelihood function
$$\log L_f^{\MN}(\Theta,b;\{(x_i,y_i,z_i,s_i=1)\}_{i=1}^n)$$
presented in Lemma <ref>, we summarize our regularized EM algorithm for the multinomial-PU model in Algorithm <ref>. In the E-step of Algorithm <ref>, we estimate the unobserved $\ind{y_i=j}$ by
for $j=0,\dots,K$ at the $m$th iteration. Since $z_i=k>0$ implies $y_i=k$, one can show that
if $k>0$. On the other hand, when $z_i=0$,

\begin{equation*}
\begin{split}
\end{split}
\end{equation*}

where the last line can be derived from the log-likelihood function for the full data presented $\log L_f^{\MN}$ in Lemma <ref>. In the M-step, we minimize the regularized full log-likelihood loss
$$-\frac{1}{n}\log L_f^{\MN}(\Theta,b;\{(x_i,\widehat{y}_i(\Theta^m,b^m),z_i,s_i=1)\}_{i=1}^n)+P^{\MN}_{\lambda}(\Theta),$$

where we use $\widehat{y}_{ik}(\Theta^m,b^m)$ as a surrogate for $\ind{y_i=k}, k=1,\dots, K$.

To solve the M-step, one can apply the proximal gradient descent algorithm [Wright et al., 2009].
For the initialization, we still consider any $\Theta^0$, $b^0$ such that the regularized observed log-likelihood loss $\mathcal{F}_n^{\MN}$ is no larger than any intercept-only model. Without any prior knowledge, one can simply choose the $(0_{p\times K},\arg\min_b\mathcal{F}_n^{\MN}(0_{p\times K},b))$ as the initializer.

Input: $\Theta^0$, $b^0$ such that $\mathcal{F}^{\MN}_n(\Theta^0,b^0)\leq\min_b\mathcal{F}_n^{\MN}(0_{p\times K},b)$

E-step: calculate $\widehat{y}(\Theta^m,b^m)\in \bbR^{n\times K}$ as follows, where each entry $\widehat{y}_{ij}$ is an estimate for $\ind{y_i=j}$:

\[\hat{y}_{ij}(\Theta^m,b^m)= \begin{cases}
\frac{e^{x_i^{\top}\Theta_j^m+b_j^m}}{1+\sum\limits_{k=1}^Ke^{x_i^{\top}\Theta_k^m+b_k^m}} & \mathrm{if}\ z_i=0\\
1 & \mathrm{if}\ z_i=j\\
0 & \mathrm{else}
\end{cases}\]

M-step: obtain $\Theta^{m+1}, b^{m+1}$ by

\begin{equation*}
\begin{split}
&-\sum_{j=1}^K\hat{y}_{ij}(\Theta^m,b^m)(x_i^\top \Theta_j+b_j)\Bigg]+P_{\lambda}^{\MN}(\Theta)
\end{split}
\end{equation*}

EM for the Multinomial-PU Model

EM algorithm for Ordinal-PU model: We also propose the following regularized EM algorithm for the ordinal-PU model.

Input: $\theta^0$ such that $\mathcal{F}^{\ON}_n(\theta^0)\leq\min_{\theta_{1:p}=0_{p\times 1}}\mathcal{F}^{\ON}_n(\theta)$

E-step: calculate $\widehat{y}(\theta^m)\in\bbR^{n\times K}$ as follows, where each entry $[\widehat{y}(\theta^m)]_{ij}$ is an estimate for $\ind{y_i=j}$:

\[\hat{y}_{ij}(\theta^m)= \begin{cases}
\frac{1}{1+e^{x^\top \theta^m_{1:p}-\sum_{l=1}^{j+1}\theta^m_{p+l}}}-\frac{1}{1+e^{x^\top \theta^m_{1:p}-\sum_{l=1}^{j}\theta^m_{p+l}}},  &\mathrm{if}\ z_i=0,\\
1,&\mathrm{if}\ z_i=j,\\
0,& \mathrm{else.}
\end{cases}\]

M step: obtain $\theta^{m+1}$ by

\begin{equation*}
\begin{split}
\theta_j^{m+1}\in\argmin\limits_{\theta}&\frac{1}{n}\sum_{i=1}^n\Bigg[-\sum_{j=1}^K\widehat{y}_{ij}(\theta^m)\log r_j(x_i, \theta)\\
\end{split}
\end{equation*}

EM for the Ordinal-PU Model

Similarly to Algorithm <ref>, in the E-step, we still estimate $\ind{y_i=j}$ by $\bbP_{\theta^m}(y_i=j|x_i,z_i,s_i=1)$ for $j=0,\dots,K$ at the $m$th iteration, which satisfies $\bbP_{\theta^m}(y_i=j|x_i,z_i,s_i=1)=\ind{j=z_i}$ if $z_i>0$. While if $z_i=0$, by the full log-likelihood function $\log L_f^{\MN}$ presented in Lemma <ref>,

\begin{equation*}
\begin{split}
=&(1+e^{x^\top \theta^m_{1:p}-\sum_{l=1}^{j+1}\theta^m_{p+l}})^{-1}-(1+e^{x^\top \theta^m_{1:p}-\sum_{l=1}^{j}\theta^m_{p+l}})^{-1}.
\end{split}
\end{equation*}

The M-step minimizes
$$-\frac{1}{n}\log L_f^{\ON}(\theta;\{(x_i,\widehat{y}_i(\theta^m),z_i,s_i=1)\}_{i=1}^n)+P^{\ON}_{\lambda}(\theta),$$
and can also be solved by the proximal gradient descent algorithm [Wright et al., 2009]

Convergence Properties: In fact, we can show the same convergence properties for both the PGD algorithms and the regularized EM algorithms. They all converge to stationary points of the corresponding penalized log-likelihood losses.

PGD for solving penalized MLE v.s. EM algorithm: From the theoretical perspective, the PGD algorithm applied on the penalized log-likelihood losses and the EM algorithms are similar since they have the same convergence guarantees. However, as suggested by our numerical experiments, the PGD algorithms are much faster while enjoying similar statistical errors to the EM algorithms under both models. Hence we recommend the practitioners to use the PGD algorithms.

If the parameter iterates $\{\theta^m\}_m)$ are generated by the proximal gradient descent algorithm (<ref>) with proper choices of step sizes or Algorithm <ref>, they would satisfy the following:

(i) The sequence $\{\theta^m\}_m$ has at least one limit point.

(ii)There exists $R^{\ON}>0$ such that all limit points of $\{\theta^m\}_m$ belong to $\Gamma^{\ON}$, the set of first order stationary points of the optimization problem $\min_{\|\theta\|_2\leq R^{\ON}}\mathcal{F}_n^{\ON}(\theta)$.

(iii)The sequence of function values $\{\mathcal{F}_n^{\ON}(\theta^m)\}_m$ is non-increasing, and $\mathcal{F}_n^{\ON}(\theta^{m+1})<\mathcal{F}_n^{\ON}(\theta^{m})$ holds if $\theta^m\notin \Gamma^{\ON}$. There exists $\widetilde{\theta}\in\Gamma^{\ON}$ such that $\{\mathcal{F}_n^{\ON}(\theta^m)\}_m$ converges monotonically to $\mathcal{F}_n^{\ON}(\widetilde{\theta})$.

Here we first prove that the sequence of function values is non-increasing for both models and algorithms. Specifically, consider the sequence $\{\cF_n^{\MN}(\Theta^m,b^m)\}_m$ where the parameter iterates $\{\Theta^m,b^m\}_m$ are generated by Algorithm <ref>. For any $m\geq 0$, one can show that

\begin{align*}
=&-\frac{1}{n}\sum_{i=1}^n\log \bbP(z_i|x_i,\Theta,b)\\
=&-\frac{1}{n}\sum_{i=1}^n\bbE_{y_i|x_i,z_i,\Theta^m,b^m}\left[\log\frac{\bbP(y_i,z_i|x_i,\Theta,b)}{\bbP(y_i|x_i,z_i,\Theta^m,b^m)}-\log \frac{\bbP(y_i|x_i,z_i,\Theta,b)}{\bbP(y_i|x_i,z_i,\Theta^m,b^m)}\right]\\
\leq&-\frac{1}{n}\sum_{i=1}^n\bbE_{y_i|x_i,z_i,\Theta^m,b^m}\left[\log\frac{\bbP(y_i,z_i|x_i,\Theta,b)}{\bbP(y_i|x_i,z_i,\Theta^m,b^m)}\right]\\
=&-\frac{1}{n}\log L_f^{\MN}(\Theta,b;\{(x_i,\widehat{y}_i(\Theta^m,b^m),z_i)\}_{i=1}^n)\\
\end{align*}

where the fifth line is due to that

\begin{align*}
&\bbE_{y_i|x_i,z_i,\Theta^m,b^m}\left[\log \frac{\bbP(y_i|x_i,z_i,\Theta,b)}{\bbP(y_i|x_i,z_i,\Theta^m,b^m)}\right]\\
\leq&\log \bbE_{y_i|x_i,z_i,\Theta^m,b^m}\left[\frac{\bbP(y_i|x_i,z_i,\Theta,b)}{\bbP(y_i|x_i,z_i,\Theta^m,b^m)}\right]\\
\end{align*}

and the equality holds if and only if $\Theta=\Theta^m$, $b=b^m$.
Here we have omitted data point $(x_i,z_i,s_i=1)$ to $(x_i,z_i)$ for simplicity. Therefore, we have

\begin{align*}
\leq &-\frac{1}{n}\log L_f^{\MN}(\Theta^{m+1},b^{m+1};\{(x_i,\widehat{y}_i(\Theta^m,b^m),z_i)\}_{i=1}^n)+P_{\lambda}^{\MN}(\Theta^{m+1})\\
&+\frac{1}{n}\sum_{i=1}^n\bbE_{y_i|x_i,z_i,\Theta^m,b^m}\left[\log P(y_i|x_i,z_i,\Theta^m,b^m)\right]\\
\leq &-\frac{1}{n}\log L_f^{\MN}(\Theta^m,b^m;\{(x_i,\widehat{y}_i(\Theta^m,b^m),z_i)\}_{i=1}^n)+P_{\lambda}^{\MN}(\Theta^m)\\
&+\frac{1}{n}\sum_{i=1}^n\bbE_{y_i|x_i,z_i,\Theta^m,b^m}\left[\log P(y_i|x_i,z_i,\Theta^m,b^m)\right]\\
\end{align*}

where the second inequality is due to the M-step in Algorithm <ref>. If $\cF_n^{\MN}(\Theta^{m+1},b^{m+1})=\cF_n^{\MN}(\Theta^m,b^m)$, then $(\Theta^{m},b^{m})$ must be a minimizer of
$$-\frac{1}{n}\log L_f^{\MN}(\Theta,b;\{(x_i,\widehat{y}_i(\Theta^m,b^m),z_i)\}_{i=1}^n)+P_{\lambda}^{\MN}(\Theta)$$
over $\bbR^{p\times K}\times \bbR^K$. Since

\begin{align*}
&\left(\nabla \cL^{\MN}(\Theta,b;\{(x_i,z_i)\}_{i=1}^n)\right)|_{\Theta^m,b^m}\\
=&\left(-\frac{1}{n}\nabla \log L_f^{\MN}(\Theta,b;\{(x_i,\widehat{y}_i(\Theta^m,b^m),z_i)\}_{i=1}^n)\right)|_{\Theta^m,b^m},
\end{align*}

$\cF_n^{\MN}(\Theta^{m+1},b^{m+1})=\cF_n^{\MN}(\Theta^m,b^m)$ implies that $(\Theta^m,b^m)\in \Gamma^{\MN}$. Using the same arguments we can show that $\{\cF_n^{\ON}(\theta^m)\}_{m}$ is also non-increasing if $\{\theta^m\}_m$ are generated by Algorithm <ref>, and $\cF_n^{\ON}(\theta^{m+1})<\cF_n^{\ON}(\theta^{m})$ if $\theta^m\in \Gamma^{\ON}$. On the other hand, if $\{\Theta^m,b^m\}_m$ (resp. $\{\theta^m\}_m$) are generated by the proximal gradient descent algorithm with each update (<ref>) (resp. (<ref>)), then existed results <cit.> suggest that with appropriately chosen step sizes $\eta_m$, $\{\cF_n^{\MN}(\Theta^m,b^m)\}$ (resp. $\{\cF_n^{\ON}(\theta^m)\}$) is non-increasing. If $\cF_n^{\MN}(\Theta^{m+1},b^{m+1})=\cF_n^{\MN}(\Theta^m,b^m)$, $(\Theta^m,b^m)$ must be a minimizer of

\eta_mP^{\MN}_{\lambda}(\Theta)+\frac{1}{2}\|[\Theta;b]-([\Theta^m;b^m]-\eta_m\nabla \cL_n^{\MN}(\Theta^m,b^m))\|_2^2.

Since any subgradient of the loss above evaluated at $(\Theta^m,b^m)$ satisfies

\begin{align*}
&\left(\nabla \eta_mP^{\MN}_{\lambda}(\Theta)+\frac{1}{2}\|[\Theta;b]-([\Theta^m;b^m]-\eta_m\nabla \cL_n^{\MN}(\Theta^m,b^m))\|_2^2\right)|_{\Theta^m,b^m}\\
=&\eta_m(\nabla P^{\MN}_{\lambda}(\Theta^m)+\nabla \cL_n^{\MN}(\Theta^m,b^m)))\\
=&\eta_m(\nabla \cF_n^{\MN}(\Theta,b))|_{\Theta^m,b^m},
\end{align*}

we also have $\cF_n^{\MN}(\Theta^{m+1},b^{m+1})<\cF_n^{\MN}(\Theta^m,b^m)$ if $(\Theta^m,b^m)\notin \Gamma^{\MN}$. Similarly we have $\cF_n^{\ON}(\theta^{m+1})<\cF_n^{\ON}(\theta^m)$ if $\theta^m\notin \Gamma^{\ON}$ when $\{\theta^m\}_m$ are generated by the proximal gradient descent algorithm (<ref>).

Now we show that $\lim_{\|\Theta\|_2+\|b\|_2\rightarrow +\infty}\cF_n^{\MN}(\Theta,b)=+\infty$ and $\lim_{\|\theta\|_2\rightarrow +\infty}\cF_n^{\ON}(\theta)=+\infty$, and hence
$$\cF_n^{\MN}(\Theta,b)\leq \argmin_{b}\cF_n^{\MN}(0_{p\times K},b),\quad \cF_n^{\ON}(\theta)\leq \argmin_{\theta_{1:p}=0_{p\times 1}}\cF_n^{\ON}(\theta)$$
implies $\|\Theta\|_2+\|b\|_2\leq R^{\MN}$ and $\|\theta\|_2\leq R^{\ON}$ for some $R^{\MN}>0$, $R^{\ON}>0$. By the initialization conditions for our proximal gradient descent algorithms and EM algorithms (Algorithm <ref>, <ref>) and the non-increasing property of $\{\cF_n^{\MN}(\Theta^m,b^m)\}$ and $\{\cF_n^{\ON}(\theta^m)\}$, it is guaranteed that all parameter iterates $\{(\Theta^m,b^m)\}_m$ belong to set $\cC^{\MN}=\{[\Theta;b]:\|\Theta\|_2+\|b\|_2\leq R^{\MN}\}$, and $\{\theta^m\}_m$ belong to $\cC^{\ON}=\{\theta:\|\theta\|_2\leq R^{\ON}\}$. Both $\cC^{\MN}$ and $\cC^{\ON}$ are compact sets, which implies (i) in Lemma <ref> and Lemma <ref>.

Due to the continuity of our loss functions $\cF_n^{\MN}$, $\cF^{\ON}$ and the fact that all parameter iterates belong to compact sets $\cC^{\MN}$, $\cC^{\ON}$, we can apply the global convergence theorem in [Zangwill, 1969] and obtain (ii) and (iii) in Proposition <ref> and Proposition <ref>.

§ DERIVATION OF LOG-LIKELIHOOD FUNCTIONS IN THE CASE-CONTROL SETTING

First we derive the conditional distribution of $z$ given $x, y$ and $s=1$ in the case-control setting, which was presented earlier in Section <ref>. Then we will present the proofs of
Lemma <ref> and Lemma <ref>.

First note that in the case-control setting, $n_u$ samples are randomly drawn from the whole population and $n_k$ samples are drawn from the population with label $k$ for $1\leq k\leq K$. Hence we have

\begin{equation*}
\frac{\bbP(z=k|s=1)}{\bbP(z=0|s=1)}=\frac{n_k}{n_u},
\end{equation*}

which implies that

\begin{equation*}
\bbP(z=k|y=k,s=1)=\frac{\bbP(z=k,s=1|y=k)}{\bbP(s=1|y=k)}=\frac{n_k}{n_k+\pi_kn_u}, \quad \bbP(z=0|y=k,s=1) = \frac{\pi_k n_u}{n_k + \pi_k n_u}.
\end{equation*}

First we derive the conditional probability mass function of any data point $(x,y,z,s=1)$ given $x$ and $s=1$. Note that there are $2K+1$ possible combinations of $(y,z)$: $y=z=0$, $y=z=j$, and $y=j, z=0$ for $1\leq j\leq K$. Then one can show that for any $1\leq j\leq K$,

\begin{equation}
\begin{split}
\bbP(y=j,z=j|x,s=1)=&\frac{\bbP(y=j|x)\bbP(z=j,s=1|y=j)}{\sum_{k=1}^K\bbP(y=k|x)\bbP(z=k,s=1|y=k)+\bbP(z=0,s=1|x)}\\
\end{split}
\end{equation}

where we have applied the fact that $\bbP(y=j|x)=\frac{e^{x^\top\Theta_j+b_j}}{1+\sum\limits_{k=1}^{K}e^{x^\top\Theta_k+b_k}}$ and
$$\frac{\bbP(z=j,s=1|y=j)}{\bbP(z=k,s=1|y=k)}=\frac{n_j/\pi_j}{n_k/\pi_k},\quad \frac{\bbP(z=j,s=1|y=j)}{\bbP(z=0,s=1|x)}=\frac{n_j/\pi_j}{n_u}.$$
Applying the above facts also leads us to

\begin{equation*}
\begin{split}
\bbP(y=j,z=0|x,s=1)=&\frac{\bbP(y=j|x)\bbP(z=0,s=1|y=j)}{\sum_{k=1}^K\bbP(y=k|x)\bbP(z=k,s=1|y=k)+\bbP(z=0,s=1|x)}\\
\bbP(y=0,z=0|x,s=1)=&\frac{\bbP(y=0|x)\bbP(z=0,s=1|y=0))}{\sum_{k=1}^K\bbP(y=k|x)\bbP(z=k,s=1|y=k)+\bbP(z=0,s=1|x)}\\
\end{split}
\end{equation*}

Thus the full log-likelihood function of $\{x_i,y_i,z_i,s_i=1\}_{i=1}^n$ is as follows:

\begin{equation}
\begin{split}
&\log L^{\MN}_f(\Theta,b;\{x_i,y_i,z_i,s_i=1\}_{i=1}^n)\\
=&\sum_{i=1}^n\left[\sum_{k=1}^K\ind{y_i=k}(x_i^\top \Theta_k+b_k)-\log\left(1+\sum\limits_{k=1}^{K}(1+\frac{n_k}{\pi_kn_u})e^{x_i^\top\Theta_k+b_k}\right)\right]\\
&+\sum_{i=1}^n\sum_{k=1}^K\ind{y_i=z_i=k}\log \frac{n_k}{\pi_kn_u}.
\end{split}
\end{equation}

To show the log-likelihood function for the observed data $\{x_i,z_i,s_i=1\}_{i=1}^n$, note that the conditional probability mass function for $z$ given $x,s=1$ is as follows:

\begin{equation}
\begin{split}
\bbP(z=j|x,s=1)=&\bbP(y=j,z=j|x,s=1)\\
=&\frac{\frac{n_j}{\pi_jn_u}e^{x^\top\Theta_j+b_j}}{1+\sum\limits_{k=1}^{K}(1+\frac{n_k}{\pi_kn_u})e^{x^\top\Theta_k+b_k}},\quad j>0,\\
\bbP(z=0|x,s=1)=&\sum_{j=1}^K\bbP(y=j,z=0|x,s=1)+\bbP(y=0,z=0|x,s=1)\\
\end{split}
\end{equation}

Therefore, the log-likelihood function of the observed data $\{x_i,z_i,s_i=1\}_{i=1}^n$ is

\begin{equation}
\begin{split}
&\log L^{\MN}(\Theta,b;\{x_i,z_i,s_i=1\}_{i=1}^n)\\
=&\sum_{i=1}^n\Bigg[\sum_{j=1}^K\ind{z_i=j}(x_i^\top \Theta_j+b_j+\log (\frac{n_j}{\pi_jn_u}))+\ind{z_i=0}\log(1+\sum_{j=1}^Ke^{x_i^\top\theta_j+b_j})\\
=&\sum_{i=1}^n\Bigg[\sum_{j=1}^K\ind{z_i=j}(x_i^\top \Theta_j+b_j+\log \frac{n_j}{\pi_j n_u}-\log(1+\sum_{j=1}^Ke^{x_i^\top\Theta_j+b_j}))\\
&-\log\left(1+\frac{\sum\limits_{k=1}^{K}\frac{n_k}{\pi_kn_u}e^{x_i^\top\Theta_k+b_k}}{1+\sum_{k=1}^K e^{x_i^\top \Theta_k+b_k}}\right)\Bigg].
\end{split}
\end{equation}

Similarly from the previous set-up, we can derive the distribution of $y,z$ given $x$ as

\begin{equation}
\begin{split}
\bbP(y=j,z=j|x,s=1)=&\frac{\bbP(y=j,z=j,s=1|x)}{\bbP(s=1|x)}\\
=&\frac{\frac{n_j}{\pi_jn_u}r_j(x,\theta)}{1+\sum\limits_{k=1}^{K}(\frac{n_k+\pi_kn_u}{\pi_kn_u})r_k(x,\Theta)},\quad j>0,\\
\bbP(y=j,z=0|x,s=1)=&\frac{r_j(x,\theta)}{1+\sum\limits_{k=1}^{K}(\frac{n_k+\pi_kn_u}{\pi_kn_u})r_k(x,\theta)},\\
\bbP(y=0,z=0|x,s=1)=&\frac{1}{1+\sum\limits_{k=1}^{K}(\frac{n_k+\pi_kn_u}{\pi_kn_u})r_k(x,\theta))},
\end{split}
\end{equation}

\begin{equation}
\begin{split}
=&(1+e^{x^\top \theta_{1:p}-\theta_{p+1}})\left[(1+e^{x^\top \theta_{1:p}-\sum_{l=1}^{j+1}\theta_l})^{-1}-(1+e^{x^\top \theta_{1:p}-\sum_{l=1}^j\theta_l})^{-1}\right].
\end{split}
\end{equation}

On the other hand, the distribution of $z$ given $x$ is

\begin{equation}
\begin{split}
\bbP(z=j|x,s=1)=&\frac{\frac{n_j}{\pi_jn_u}r_j(x,\theta)}{1+\sum\limits_{k=1}^{K}(\frac{n_k+\pi_kn_u}{\pi_kn_u})r_k(x,\theta)},\quad j>0,\\
\bbP(z=0|x,s=1)=&\frac{1+\sum_{k=1}^Kr_k(x,\theta)}{1+\sum\limits_{k=1}^{K}(\frac{n_k+\pi_kn_u}{\pi_kn_u})r_k(x,\theta)}.
\end{split}
\end{equation}

Thus the full log-likelihood function of $\{(x_i,y_i,z_i,s_i=1)\}_{i=1}^n$ is as follows:

\begin{equation}
\begin{split}
&\log L_f^{\ON}(\theta;\{(x_i,y_i,z_i,s_i=1)\}_{i=1}^n)\\
=&\sum_{i=1}^n\left[\sum_{j=1}^K\ind{y_i=j}\log r_j(x_i, \theta)-\log\left(1+\sum\limits_{k=1}^{K}\frac{n_k+\pi_kn_u}{\pi_kn_u}r_k(x_i,\theta)\right)\right]\\
&+\sum_{i=1}^n\sum_{j=1}^K\ind{y_i=z_i=j}\log \frac{n_j}{\pi_jn_u},
\end{split}
\end{equation}

and the log-likelihood of the observed data $\{(x_i,z_i,s_i=1)\}_{i=1}^n$ is

\begin{equation*}
\begin{split}
&\log L^{\ON}(\theta;\{(x_i,z_i,s_i=1)\}_{i=1}^n)\\
=&\sum_{i=1}^n\Bigg[\sum_{j=1}^K\ind{z_i=j}(\log r_j(x_i,\theta)-\log (\frac{n_j}{\pi_jn_u}))+\ind{z_i=0}\log(1+\sum_{j=1}^Kr_j(x_i,\theta))\\
=&\sum_{i=1}^n \left[\sum_{k=1}^K \ind{z_i=k}[f(\log r(x_i,\theta))]_k-\log\left(1+\sum_{j=1}^Ke^{[f(\log r(x_i,\theta))]_k}\right)\right],
\end{split}
\end{equation*}

where $f:\bbR^K\rightarrow \bbR^K$ satisfies $(f^{\ON}(u))_k=u_k+\log \frac{n_k}{\pi_kn_u}-\log(1+\sum_{j=1}^Ke^{u_j})$.

§ ALGORITHMS FOR THE PU MODELS UNDER THE SINGLE-TRAINING-SET SCENARIO

Here we consider the single-training-set scenario [Elkan and Noto, 2008], a different setting from the case-control scenario for PU models. In particular, we still assume that the multinomial and ordinal distributions for the true labeled data $(x_i,y_i)$, while letting the observation $z_i$ be unlabeled with some constant probability. Formally, for any $j,k\geq 0$,

\begin{align}\label{eq:random-missing}
\bbP(z_i=k|x_i,y_i=j)=\begin{cases}\pi^{\st}_j,&\text{if $j=k>0$,}\\
1-\pi^{\st}_j,& \text{if $j>0$, $k=0$,}\\
1,&\text{if $j=k=0$,}\\
\end{cases}
\end{align}

where $\pi^{\st}\in(0,1)^K$ indicates the probabilities of positive samples being labeled, and the superscript $\st$ refers to the single-training-set scenario.

§.§ Multinomial-PU Model

The following lemma presents the log-likelihood functions for the observed data $\{(x_i,z_i)\}$ and for the full data $\{(x_i,y_i,z_i)\}$ under the multinomial-PU model and the single-training-set scenario.

When the observed responses $z_i$ is generated according to (<ref>), the log-likelihood function for the observed presence-only data $\{x_i,z_i\}_{i=1}^n$ under the multinomial model is

\begin{equation}\label{eq:MN-log-likelihood-st}
\begin{split}
&\log L^{\MN,\st}(\Theta,b;\{(x_i,z_i)\}_{i=1}^n)\\
=&\sum_{i=1}^n \left[\sum_{k=1}^K \ind{z_i=k}[\widetilde{f}(\Theta^\top x_i+b+\log (1-\pi^{st}))]_k-\log\left(1+\sum_{k=1}^Ke^{[\widetilde{f}(\Theta^\top x_i+b+\log (1-\pi^{st}))]_k}\right)\right],
\end{split}
\end{equation}

where $\widetilde{f}:\bbR^K\rightarrow \bbR^K$ satisfies $(\widetilde{f}(u))_k=u_k+\log \frac{\pi^{\st}_k}{1-\pi^{\st}_k}-\log(1+\sum_{j=1}^Ke^{u_j})$. Here $\log (1-\pi^{\st})=(\log(1-\pi^{\st}_1),\dots,\log(1-\pi^{\st}_K))^\top\in\bbR^K$.
The log-likelihood function for the full data $\{x_i,y_i,z_i\}_{i=1}^n$ is

\begin{equation}\label{eq:MN-full-log-likelihood-st}
\begin{split}
&\log L_f^{\MN,\st}(\Theta,b;\{(x_i,y_i,z_i)\}_{i=1}^n)\\
\bigg[\sum_{k=1}^K \ind{y_i=k}(x_i^\top \Theta_k+b_k+\log(1-\pi^{\st}_k))-\log\left(1+\sum_{k=1}^Ke^{x_i^\top\Theta_k+b_k}\right)\\
\end{split}
\end{equation}

where $\Theta_k$ is the $k$th column of $\Theta$.

Under the multinomial-PU model with single-training-set scenario, one can show that for any data point $(x,y,z)$ and $1\leq j\leq K$,

\begin{align*}
\bbP(y=j,z=j|x)=&\bbP(y=j|x)\bbP(z=j|y=j)=\frac{\pi^{\st}_je^{x^\top\Theta_j+b_j}}{1+\sum_{k=1}^Ke^{x^\top\Theta_k+b_k}},\\
\bbP(y=j,z=0|x)=&\bbP(y=j|x)\bbP(z=0|y=j)=\frac{(1-\pi^{\st}_j)e^{x^\top\Theta_j+b_j}}{1+\sum_{k=1}^Ke^{x^\top\Theta_k+b_k}},\\
\bbP(y=0,z=0|x)=&\bbP(y=0|x)=\frac{1}{1+\sum_{k=1}^Ke^{x^\top\Theta_k+b_k}}.
\end{align*}

While for the observed PU data $(x,z)$, we have

\begin{align*}
\bbP(z=j|x)=&\bbP(y=j,z=j|y=j)=\frac{\pi^{\st}_je^{x^\top\Theta_j+b_j}}{1+\sum_{k=1}^Ke^{x^\top\Theta_k+b_k}},\\
\bbP(z=0|x)=&\sum_{j=1}^K\bbP(y=j,z=0|y=j)+\bbP(y=0,z=0|x)=\frac{1+\sum_{j=1}^K(1-\pi^{\st}_j)e^{x^\top\Theta_j+b_j}}{1+\sum_{k=1}^Ke^{x^\top\Theta_k+b_k}}.
\end{align*}

Hence the full log-likelihood for samples $\{(x_i,y_i,z_i)\}_{i=1}^n$ can be written as follows:

\begin{align*}
&\log L_f^{\MN,\st}(\Theta,b;\{(x_i,y_i,z_i)\}_{i=1}^n)\\
=&\sum_{i=1}^n\left[\sum_{j=1}^K\ind{y_i=j}(x_i^\top \Theta_j+b_j+\log(1-\pi^{\st}_j))-\log(1+\sum_{k=1}^Ke^{x_i^\top\Theta_k+b_k})\right]\\
\end{align*}

While for the log-likelihood for the observed PU data $\{(x_i,z_i)\}_{i=1}^n$ is

\begin{align*}
&\log L^{\MN,\st}(\Theta,b;\{(x_i,z_i)\}_{i=1}^n)\\
=&\sum_{i=1}^n\Bigg[\sum_{j=1}^K\ind{z_i=j}(x_i^\top\Theta_j+b_j+\log \pi^{\st}_j-\log(1+\sum_{k=1}^K(1-\pi^{\st}_k)e^{x_i^\top\Theta_k+b_k})))\\
&-\log(1+\frac{\sum_{k=1}^K\pi^{\st}_ke^{x_i^\top\Theta_k+b_k}}{1+\sum_{k=1}^K(1-\pi^{\st}_k)e^{x_i^\top\Theta_k+b_k}}) \Bigg]\\
=&\sum_{i=1}^n \Bigg[\sum_{k=1}^K \ind{z_i=k}[\widetilde{f}(\Theta^\top x_i+b+\log (1-\pi^{st}))]_k-\log\left(1+\sum_{j=1}^Ke^{[\widetilde{f}(\Theta^\top x_i+b+\log (1-\pi^{st}))]_k}\right)\Bigg],
\end{align*}

Similarly to the estimation under the case-control case (discussed in Section <ref>), one can apply the proximal gradient descent algorithm on
$$\cF_n^{\MN,\st}(\Theta,b)=-\frac{1}{n}\log L^{\MN,\st}(\Theta,b;\{(x_i,z_i)\}_{i=1}^n)+P_{\lambda}^{\MN}(\Theta)$$
for the estimation of $\Theta$ and $b$. We can also derive the EM algorithm under this setting (see Algorithm <ref>).

Input: $\Theta^0$, $b^0$ such that $\mathcal{F}^{\MN,\st}_n(\Theta^0,b^0)\leq\min_b\mathcal{F}_n^{\MN,\st}(0_{p\times K},b)$

E-step: calculate $\widehat{y}(\Theta^m,b^m)\in \bbR^{n\times K}$ as follows, where each entry $\widehat{y}_{ij}$ is an estimate for $\ind{y_i=j}$:

\[\hat{y}_{ij}(\Theta^m,b^m)= \begin{cases}
\frac{e^{x_i^{\top}\Theta_j^m+b_j^m+\log(1-\pi^{\st}_j)}}{1+\sum\limits_{k=1}^Ke^{x_i^{\top}\Theta_k^m+b_k^m+\log(1-\pi^{\st}_k)}} & \mathrm{if}\ z_i=0\\
1 & \mathrm{if}\ z_i=j\\
0 & \mathrm{else}
\end{cases}\]

M-step: obtain $\Theta^{m+1}, b^{m+1}$ by

\begin{equation*}
\begin{split}
&-\sum_{j=1}^K\hat{y}_{ij}(\Theta^m,b^m)(x_i^\top \Theta_j+b_j)\Bigg]+P_{\lambda}^{\MN}(\Theta)
\end{split}
\end{equation*}

EM for the Multinomial-PU Model under the single-training-set scenario

Theoretical properties: One key observation is that the multinomial-PU model under the single-training-set scenario can be viewed as a simple reparameterization of the case-control setting. More specifically, to obtain the log-likelihood functions under the single-training-set scenario, one can simply substitute $\frac{n_k}{\pi_kn_u}$ in $\log L_f^{\MN}$ and $\log L^{\MN}$ by $\frac{\pi^{\st}_k}{1-\pi^{\st}_k}$, and change $b_k$ to $b_k+\log(1-\pi^{\st}_k)$ for $1\leq k\leq K$. Hence all the optimization and statistical guarantees presented in Proposition <ref> and Theorem <ref> still hold for this setting, as long as we change the condition on $\frac{n_k}{\pi^{\st}_kn_u}$ to the condition on $\frac{\pi^{\st}_k}{1-\pi^{\st}_k}$.

§.§ Ordinal-PU Model

The following lemma presents the log-likelihood functions for the observed data $\{(x_i,z_i)\}$ and for the full data $\{(x_i,y_i,z_i)\}$ under the ordinal-PU model and the single-training-set scenario.

When the observed responses $z_i$ is generated according to (<ref>), the log-likelihood function for the observed presence-only data $\{x_i,z_i\}_{i=1}^n$ under the ordinal-PU model is

\begin{equation}\label{eq:ON-log-likelihood-st}
\begin{split}
&\log L^{\ON,\st}(\theta;\{(x_i,z_i)\}_{i=1}^n)\\
=&\sum_{i=1}^n \left[\sum_{k=1}^K \ind{z_i=k}[\widetilde{f}(\log r(x_i,\theta)+\log(1-\pi^{\st}))]_k-\log\left(1+\sum_{k=1}^Ke^{[\widetilde{f}(\log r(x_i,\theta)+\log(1-\pi^{\st}))]_k}\right)\right],
\end{split}
\end{equation}

where $\widetilde{f}$ is defined in Lemma <ref> and $r:\bbR^p\times \bbR^{p+K}\rightarrow \bbR^{K}$ is defined in Lemma <ref>.
The log-likelihood function for the full data $\{x_i,y_i,z_i\}_{i=1}^n$ is

\begin{equation}\label{eq:ON-full-log-likelihood-st}
\begin{split}
&\log L_f^{\ON,\st}(\theta;\{(x_i,y_i,z_i)\}_{i=1}^n)\\
\bigg[\sum_{k=1}^K \ind{y_i=k}(\log r_k(x_i,\theta)+\log(1-\pi^{\st}_k)-\log\left(1+\sum_{k=1}^Kr_k(x_i,\theta)\right)\\
\end{split}
\end{equation}

Similarly to the proof of Lemma <ref>, we first write out the joint probability mass for $(x,y,z)$ as follows:

\begin{align*}
\bbP(y=j,z=j|x)=&\frac{\pi^{\st}_j r_j(x,\theta)}{1+\sum_{k=1}^K r_k(x,\theta)},\\
\bbP(y=j,z=0|x)=&\frac{(1-\pi^{\st}_j) r_j(x,\theta)}{1+\sum_{k=1}^K r_k(x,\theta)},\\
\bbP(y=0,z=0|x)=&\frac{1}{1+\sum_{k=1}^K r_k(x,\theta)}.
\end{align*}

The joint probability mass for $(x,z)$ would then become:

\begin{align*}
\bbP(z=j|x)=&\frac{\pi^{\st}_j r_j(x,\theta)}{1+\sum_{k=1}^K r_k(x,\theta)},\\
\bbP(z=0|x)=&\frac{1+\sum_{k=1}^K(1-\pi^{\st}_k) r_k(x,\theta)}{1+\sum_{k=1}^K r_k(x,\theta)}.
\end{align*}

Therefore, one can directly write out the log-likelihood functions as shown in Lemma <ref>.

$$\cF^{\ON,\st}_n(\theta)=-\frac{1}{n}\log L^{\ON,\st}(\theta;\{(x_i,z_i)\}_{i=1}^n)+P^{\ON}_{\lambda}(\theta_{1:p}),$$
then we can estimate $\theta^*$ by applying the proximal gradient descent algorithm on $\cF^{\ON,\st}_n(\theta)$. The regularized EM algorithm under this setting is summarized in Algorithm <ref>.

Input: $\theta^0$ such that $\mathcal{F}^{\ON,\st}_n(\theta^0)\leq\min_{\theta_{1:p}=0_{p\times 1}}\mathcal{F}^{\ON,\st}_n(\theta)$

E-step: calculate $\widehat{y}(\theta^m)\in\bbR^{n\times K}$ as follows, where each entry $[\widehat{y}(\theta^m)]_{ij}$ is an estimate for $\ind{y_i=j}$:

\[\hat{y}_{ij}(\theta^m)= \begin{cases}
\frac{(1-\pi^{\st}_j)r_j(x_i,\theta^m)}{1+\sum_{k=1}^K(1-\pi^{\st}_k)r_k(x_i,\theta^m)},  &\mathrm{if}\ z_i=0,\\
1,&\mathrm{if}\ z_i=j,\\
0,& \mathrm{else.}
\end{cases}\]

M step: obtain $\theta^{m+1}$ by

\begin{equation*}
\begin{split}
\theta_j^{m+1}\in\argmin\limits_{\theta}&\frac{1}{n}\sum_{i=1}^n\Bigg[-\sum_{j=1}^K\widehat{y}_{ij}(\theta^m)\log r_j(x_i, \theta)\\
\end{split}
\end{equation*}

EM for the Ordinal-PU Model under the single-training-set scenario

Theoretical properties: Using similar arguments as the proof of Proposition <ref>, we can still show the convergence of these two algorithms to stationary points of the optimization problem $\min_{\|\theta\|_2\leq C} \cF_n^{\ON,\st}(\theta)$ for some constant $C>0$. Although we do not have a rigorous statistical error bound for the stationary points under this setting, we conjecture that they would satisfy similar statistical properties from Theorem <ref>. A rigorous proof for this conjecture is left as future work.

§ NUMERIC DETAILS

Controlling the true prevalence of different categories: In our comparative simulation studies, we controlled the prevalence of different categories in an approximate way by setting the intercepts appropriately. In particular, for the experiments investigating the effect of the prevalence, we choose the intercepts to ensure that the intercept-only model have prevalence $\pi_0 \in \{0.1,\,0.2,\dots,0.7\}$ and $\pi_1 = \pi_2 = \frac{1-\pi_0}{2}$. Since all features $X_j$ are generated from a mean zero Gaussian distribution, we would expect that the prevalence of the full models are not too different from the prevalence of the intercept-only models. The prevalence of the full models are then estimated from the generated data and presented in the figures. While for the experiments investigating the effect of sampling ratio $\frac{n_u}{n}$, we ensure the intercept-only models have $30\%$ unlabeled data, $35\%$ samples of label $1$ and $2$. The resulting prevalence of the full multinomial model is $\pi_0=0.2584,\,\pi_1 = 0.3677,\,\pi_2 = 0.3739$, and that of the full ordinal model is $\pi_0=0.3753,\,\pi_1 = 0.2153,\,\pi_2 = 0.4095$. Hence the whole population is reasonably balanced.

§ ACKNOWLEDGMENTS

We would like to thank David Neiman for some preliminary empirical studies and real data exploration. GR and LZ were partially supported by NSF-DMS 1811767, NIH R01 GM131381-01. GR was also partially supported by NSF-DMS 1839338.

[Agresti, 2010]
Alan Agresti.
Analysis of ordinal categorical data, volume 656.
John Wiley & Sons, 2010.

[Archer et al., 2014]
Kellie J Archer, Jiayi Hou, Qing Zhou, Kyle Ferber, John G Layne, and Amanda E
ordinalgmifs: An r package for ordinal regression in high-dimensional
data settings.
Cancer informatics, 13:0 CIN–S20806, 2014.

[Beck, 2017]
Amir Beck.
First-order methods in optimization.
SIAM, 2017.

[Bertens et al., 2016]
Loes CM Bertens, Karel GM Moons, Frans H Rutten, Yvonne van Mourik, Arno W
Hoes, and Johannes B Reitsma.
A nomogram was developed to enhance the use of multinomial logistic
regression modeling in diagnostic research.
Journal of clinical epidemiology, 71:0 51–57, 2016.

[Du Plessis et al., 2015]
Marthinus Du Plessis, Gang Niu, and Masashi Sugiyama.
Convex formulation for learning from positive and unlabeled data.
In International conference on machine learning, pages
1386–1394. PMLR, 2015.

[Dua and Graff, 2017]
Dheeru Dua and Casey Graff.
UCI machine learning repository, 2017.
URL <http://archive.ics.uci.edu/ml>.

[Elkan and Noto, 2008]
Charles Elkan and Keith Noto.
Learning classifiers from only positive and unlabeled data.
In Proceedings of the 14th ACM SIGKDD international conference
on Knowledge discovery and data mining, pages 213–220, 2008.

[Fan et al., 2010]
Jianqing Fan, Rui Song, et al.
Sure independence screening in generalized linear models with
The Annals of Statistics, 380 (6):0
3567–3604, 2010.

[Fernandes et al., 2017]
Kelwin Fernandes, Jaime S Cardoso, and Jessica Fernandes.
Transfer learning with partial observability applied to cervical
cancer screening.
In Iberian conference on pattern recognition and image
analysis, pages 243–250. Springer, 2017.

[Fields, 2001]
Peter A Fields.
Protein function at thermal extremes: balancing stability and
Comparative Biochemistry and Physiology Part A: Molecular &
Integrative Physiology, 1290 (2-3):0 417–431, 2001.

[Friedman et al., 2010]
Jerome Friedman, Trevor Hastie, and Robert Tibshirani.
A note on the group lasso and a sparse group lasso.
arXiv preprint arXiv:1001.0736, 2010.

[Kakade et al., 2010]
Sham Kakade, Ohad Shamir, Karthik Sindharan, and Ambuj Tewari.
Learning exponential families in high-dimensions: Strong convexity
and sparsity.
In Proceedings of the thirteenth international conference on
artificial intelligence and statistics, pages 381–388. JMLR Workshop and
Conference Proceedings, 2010.

[Krishnapuram et al., 2005]
Balaji Krishnapuram, Lawrence Carin, Mário AT Figueiredo, and Alexander J
Sparse multinomial logistic regression: Fast algorithms and
generalization bounds.
IEEE transactions on pattern analysis and machine
intelligence, 270 (6):0 957–968, 2005.

[Lancaster and Imbens, 1996]
Tony Lancaster and Guido Imbens.
Case-control studies with contaminated controls.
Journal of Econometrics, 710 (1-2):0 145–160,

[Lee and Wang, 2020]
Chanwoo Lee and Miaoyan Wang.
Tensor denoising and completion based on ordinal observations.
In International Conference on Machine Learning, pages
5778–5788. PMLR, 2020.

[Li et al., 2020]
Yuan Li, Benjamin Mark, Garvesh Raskutti, Rebecca Willett, Hyebin Song, and
David Neiman.
Graph-based regularization for regression problems with alignment and
highly correlated designs.
SIAM journal on mathematics of data science, 20
(2):0 480–504, 2020.

[Little and Rubin, 2019]
Roderick JA Little and Donald B Rubin.
Statistical analysis with missing data, volume 793.
John Wiley & Sons, 2019.

[Liu et al., 2003]
Bing Liu, Yang Dai, Xiaoli Li, Wee Sun Lee, and Philip S Yu.
Building text classifiers using positive and unlabeled examples.
In Third IEEE International Conference on Data Mining, pages
179–186. IEEE, 2003.

[Loh and Wainwright, 2013]
Po-Ling Loh and Martin J Wainwright.
Regularized m-estimators with nonconvexity: Statistical and
algorithmic theory for local optima.
arXiv preprint arXiv:1305.2436, 2013.

[Loh et al., 2017]
Po-Ling Loh, Martin J Wainwright, et al.
Support recovery without incoherence: A case for nonconvex
Annals of Statistics, 450 (6):0 2455–2482,

[Maurer, 2016]
Andreas Maurer.
A vector-contraction inequality for rademacher complexities.
In International Conference on Algorithmic Learning Theory,
pages 3–17. Springer, 2016.

[McCullagh, 1980]
Peter McCullagh.
Regression models for ordinal data.
Journal of the Royal Statistical Society: Series B
(Methodological), 420 (2):0 109–127, 1980.

[McDiarmid, 1989]
Colin McDiarmid.
On the method of bounded differences.
Surveys in combinatorics, 1410 (1):0 148–188,

[Negahban et al., 2009]
Sahand Negahban, Bin Yu, Martin J Wainwright, and Pradeep K Ravikumar.
A unified framework for high-dimensional analysis of $ m $-estimators
with decomposable regularizers.
In Advances in neural information processing systems, pages
1348–1356, 2009.

[Romero et al., 2015]
Philip A Romero, Tuan M Tran, and Adam R Abate.
Dissecting enzyme function with microfluidic-based deep mutational
Proceedings of the National Academy of Sciences, 1120
(23):0 7159–7164, 2015.

[Sivapriya et al., 2015]
TR Sivapriya, AR Kamal, and P Thangaiah.
Ensemble merit merge feature selection for enhanced multinomial
classification in alzheimer’s dementia.
Computational and mathematical methods in medicine, 2015,

[Song and Raskutti, 2020]
Hyebin Song and Garvesh Raskutti.
Pulasso: High-dimensional variable selection with presence-only data.
Journal of the American Statistical Association, 1150
(529):0 334–347, 2020.

[Song et al., 2020]
Hyebin Song, Ran Dai, Garvesh Raskutti, and Rina Foygel Barber.
Convex and non-convex approaches for statistical inference with
class-conditional noisy labels.
Journal of Machine Learning Research, 210
(168):0 1–58, 2020.

[Song et al., 2021]
Hyebin Song, Bennett J Bremer, Emily C Hinds, Garvesh Raskutti, and Philip A
Inferring protein sequence-function relationships with large-scale
positive-unlabeled learning.
Cell systems, 120 (1):0 92–101, 2021.

[Tibshirani, 1996]
Robert Tibshirani.
Regression shrinkage and selection via the lasso.
Journal of the Royal Statistical Society: Series B
(Methodological), 580 (1):0 267–288, 1996.

[Vaart and Wellner, 1997]
AW van der Vaart and Jon A Wellner.
Weak convergence and empirical processes with applications to
Journal of the Royal Statistical Society-Series A Statistics in
Society, 1600 (3):0 596–608, 1997.

[Van de Geer et al., 2008]
Sara A Van de Geer et al.
High-dimensional generalized linear models and the lasso.
Annals of Statistics, 360 (2):0 614–645,

[Vershynin, 2010]
Roman Vershynin.
Introduction to the non-asymptotic analysis of random matrices.
arXiv preprint arXiv:1011.3027, 2010.

[Vershynin, 2018]
Roman Vershynin.
High-dimensional probability: An introduction with applications
in data science, volume 47.
Cambridge university press, 2018.

[Wainwright, 2019]
Martin J Wainwright.
High-dimensional statistics: A non-asymptotic viewpoint,
volume 48.
Cambridge University Press, 2019.

[Ward et al., 2009]
Gill Ward, Trevor Hastie, Simon Barry, Jane Elith, and John R Leathwick.
Presence-only data and the em algorithm.
Biometrics, 650 (2):0 554–563, 2009.

[Wright et al., 2009]
Stephen J Wright, Robert D Nowak, and Mário AT Figueiredo.
Sparse reconstruction by separable approximation.
IEEE Transactions on signal processing, 570
(7):0 2479–2493, 2009.

[Wurm et al., 2017]
Michael J Wurm, Paul J Rathouz, and Bret M Hanlon.
Regularized ordinal regression and the ordinalnet r package.
arXiv preprint arXiv:1706.05003, 2017.

[Xu et al., 2017]
Jason Xu, Eric C Chi, and Kenneth Lange.
Generalized linear model regression under distance-to-set penalties.
arXiv preprint arXiv:1711.01341, 2017.

[Zangwill, 1969]
Willard I Zangwill.
Nonlinear programming: a unified approach, volume 52.
Prentice-Hall Englewood Cliffs, NJ, 1969.