# Two-level monotonic multistage recommender systems††thanks: Research
supported in part by NSF grants DMS-1712564, DMS-1721216, DMS-1952539, and NIH
grants 1R01GM126002, R01HL105397, R01AG069895.

Ben Dai1,3 and Xiaotong Shen1 and Wei Pan2
1Department of Statistics, University of Minnesota
2Division of Biostatistics, University of Minnesota
3Department of Statistics, The Chinese University of Hong Kong
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A recommender system learns to predict the user-specific preference or
intention over many items simultaneously for all users, making personalized
recommendations based on a relatively small number of observations. One
central issue is how to leverage three-way interactions, referred to as user-
item-stage dependencies on a monotonic chain of events, to enhance the
prediction accuracy. A monotonic chain of events occurs, for instance, in an
article sharing dataset, where a “follow” action implies a “like” action,
which in turn implies a “view” action. In this article, we develop a
multistage recommender system utilizing a two-level monotonic property
characterizing a monotonic chain of events for personalized prediction.
Particularly, we derive a large-margin classifier based on a nonnegative
additive latent factor model in the presence of a high percentage of missing
observations, particularly between stages, reducing the number of model
parameters for personalized prediction while guaranteeing prediction
consistency. On this ground, we derive a regularized cost function to learn
user-specific behaviors at different stages, linking decision functions to
numerical and categorical covariates to model user-item-stage interactions.
Computationally, we derive an algorithm based on blockwise coordinate descent.
Theoretically, we show that the two-level monotonic property enhances the
accuracy of learning as compared to a standard method treating each stage
individually and an ordinal method utilizing only one-level monotonicity.
Finally, the proposed method compares favorably with existing methods in
simulations and an article sharing benchmark.

Keywords: Monotonicity, Nonconvex Minimization, Explicit/implicit Feedback,
Nonnegative Latent Factor Model

##  1 Introduction

A multistage recommender system on a monotonic chain of events predicts a
user’s preference of a large collection of items based on only a few user-item
feedback at multiple stages, where a user’s positive feedback of subsequent
stages can only be observed given his/her positive feedback at present stages.
As a result, a user may exhibit an increasing level of intention and
preference given the feedback, ranging from the most implicit to the most
explicit. It has been widely used for personalized prediction in e-commerce
and social networks as well as individual drug responses over multiple phases
in personalized medicine (Suphavilai et al., 2018). In such a situation, the
objective of a multistage recommender system is to predict a user’s subsequent
behavior from his/her previous feedback with a high percentage of missing
observations. For instance, as displayed in Figure 1, the Deskdrop article
sharing dataset111https://www.kaggle.com/gspmoreira/articles-sharing-reading-
from-cit-deskdrop consists of a sequence of actions ranging from view, like,
to follow, where only less than 0.2% of values are observations with a total
of $73,000$ users and more than 3,000 articles. A user likes an article only
if this user has viewed it, and the user may follow an article only when
he/she has liked it. Thus, three pairwise predictions are made based on three
pairs of present and subsequent stages are performed, namely, the prediction
of if the user will like an article given that he has viewed, that if he/she
will follow given that he has viewed, and that if he/she will follow given
that he has liked.

Figure 1: Demonstration of multistage prediction for Deskdrop dataset, a user
may like an article only if this user has viewed it, and the user may follow
an article only if this user has liked it. Three stagewise predictions based
on different pairs of a given stage and a subsequent stage are considered:
given view to predict like, given view to predict follow, and given like to
predict follow.

One key characteristic of multistage prediction on a monotonic chain of events
is a two-level monotonic property, as precisely defined in (1). On the one
hand, monotonicity occurs over subsequent stages given a present stage. For
instance, based on the fact that a user has viewed the article, this user
won’t follow it if he/she doesn’t like it. On the other hand, monotonicity
exhibits over present stages given a subsequent stage in a reversed order.
That said, if the user has decided to follow only he/she viewed it, then the
user will follow an article when he/she has liked.

###  1.1 Literature review and our contributions

Despite recent progress in recommender systems, multistage personalized
prediction on a monotonic chain of events remains largely unexplored. A
standard method is to predict a user’s behaviors at a subsequent stage given a
present stage individually, ignoring the monotonic property, as described in
(14). In this sense, essentially all existing methods of recommender systems
for a single-stage are applicable, including latent factor models (Koren,
2008; Dai et al., 2019b), logistic matrix factorization (Johnson, 2014),
tensor factorization(Bi et al., 2018), and deep neural networks (He et al.,
2017). The worth of note is that ordinal regression/classification may be
adopted by treating subsequent stages as an ordered class given a present
stage, as in (15). For example, in Wan and McAuley (2018); Tran et al. (2012),
latent factor models are developed based on an ordered logit model. Yet, a
two-level monotonic property introduced by a monotonic chain of events is not
met, which considers only forward monotonicity of subsequent stages given the
present one but not the backward monotonicity over present stages given the
subsequent one. As to be seen, the two-level monotonicity yields not only a
high accuracy of prediction but also prediction consistency (1) among
different stage pairs. In practice, prediction consistency is indeed desirable
for decision-making.

This article develops a multistage recommender system that models user-item-
stage interactions and integrates the two-level monotonic property into
personalized prediction on a monotonic chain of events. Three key
contributions in this article can be summarized:

  * •

The proposed method can produce a recommendation for any subsequent stage
given observations at any present stage, which is highly demanded in real
applications; see Figure 1. Yet, most conventional recommender systems focus
on a fixed present stage.

  * •

A novel multistage loss function is proposed to treat multistage prediction
for any pair of present and subsequent stages. This loss function admits user-
item interactions observed at different stages and evaluates the prediction
accuracy on all subsequent stages.

  * •

The two-level monotonic property is fully accounted for by our nonnegative
additive latent factor model based on the Bayes rule in Lemma 1. As a result,
it substantially reduces model parameters and most importantly ensures
prediction consistency across different stages. By comparison, none of
aforementioned methods can guarantee prediction consistency, c.f., Tables 2
and 3.

  * •

An algorithm is developed to implement the proposed method based on blockwise
coordinate descent. Moreover, a learning theory is established to quantify the
generalization error to demonstrate the benefits of modeling user-item-stage
interactions based on the two-level monotonic property.

##  2 Multistage classification on a monotonic chain

Consider a recommender system in which triplets
$(\Delta_{ij},\bm{X}_{ij},\bm{Y}_{ij})$, with
$\bm{Y}_{ij}=(Y^{1}_{ij},\cdots,Y^{T}_{ij})$ are observed for user $i$ on item
$j$ at stage $t$; $1\leq i\leq n$, $1\leq j\leq m$, $1\leq t\leq T$. Here $n$
and $m$ are the number of users and items, respectively. $\Delta_{ij}=1/0$
indicates if $(\bm{X}_{ij},\bm{Y}_{ij})$ is observed or missing and
$\Omega_{t}=\\{(i,j):\Delta_{ij}=1,y^{t}_{ij}=1\\}$ is an index set of
observed positive feedback at stage $t$, where $Y^{t}_{ij}=\pm 1$ with 1/-1
indicating positive/negative feedback at stage $t$, and ⊤ denotes transpose.
Moreover,
$\bm{x}_{ij}=(\bm{u}_{i}^{\top},\bm{v}_{j}^{\top},\bm{s}_{i}^{\top},\bm{o}^{\top}_{j})^{\top}$
consists of both numerical and categorical features, where
$\bm{u}_{i}\in[0,1]^{p_{1}}$ and $\bm{v}_{j}\in[0,1]^{p_{2}}$ are user-
specific and item-specific numerical predictors, respectively, and
normalization is performed for each feature to scale to $[0,1]$ otherwise, and
moreover $\bm{s}_{i}=(s_{i1},\cdots,s_{id_{1}})^{\top}$ with
$s_{il}\in\\{1,\cdots,n_{l}\\}$ and $\bm{o}_{j}=(o_{j1},\cdots,o_{jd_{2}})$
with $o_{jl}\in\\{1,\cdots,m_{l}\\}$ are $d_{1}$-dimensional and
$d_{2}$-dimensional user-specific and item-specific categorical predictor
vectors. For instance, in the Deskdrop
dataset222https://www.kaggle.com/gspmoreira/articles-sharing-reading-from-cit-
deskdrop, $\bm{v}_{j}$ is a numerical embedding for the content, $\bm{s}_{i}$
consists of “person Id”, “user agent”, “user region”, “user country”,
$\bm{o}_{j}$ consists of “content Id”, “author Id”, “language”, and
$Y^{1}_{ij}$ and $Y^{2}_{ij}$ indicate if article $j$ is liked and followed by
user $i$, respectively. One important characteristic of this kind of data is
that subsequent events will not occur if one present event does not occur. For
instance, if a user does not view an item, then subsequent events of like and
follow will not occur.

On this ground, we define a monotonic behavior chain as follows:

$\displaystyle Y^{t}_{ij}=-1\text{ if }Y^{t-1}_{ij}=-1;\quad 1\leq i\leq
m,1\leq j\leq n,t=1,\cdots,T,$ (1)

which says that event $Y^{t-1}_{ij}=-1$ implies that any subsequent event
$Y^{t}_{ij}=-1$ must occur, which encodes a certain causal relation as defined
by the local Markov property (Edwards, 2012) in a directed acyclic graphical
model.

###  2.1 Two-level monotonicity

Multistage classification learns to predict the outcome of $Y_{ij}^{t}$ at a
subsequent stage $t$ given observations at a present stage
$Y_{ij}^{t^{\prime}}$ and $\bm{x}_{ij}$; $0\leq t^{\prime}<t\leq T$. With the
convention, we set $Y^{0}_{ij}=1$. To predict, we introduce a decision
function
$\phi_{ij}^{t^{\prime}t}=\phi_{ij}^{t^{\prime}t}(\bm{x}_{ij},y^{t^{\prime}}_{ij})$
for classification, which depends on observations $\bm{x}_{ij}$ and
$y^{t^{\prime}}_{ij}$. For user $i$ on item $j$ across $T$ stages, decision
functions can be expressed as $\bm{\phi}_{ij}=(\phi_{ij}^{t^{\prime}t})_{0\leq
t^{\prime}<t\leq T}$. For all users over items across $T$ stages, decision
functions are $\bm{\phi}=(\bm{\phi}_{ij})_{1\leq i\leq n;1\leq j\leq m}$.

To evaluate the overall performance of $\bm{\phi}$, we define the multistage
misclassification error

$\displaystyle
e(\bm{\phi})=(nm)^{-1}\sum_{i=1}^{n}\sum_{j=1}^{m}e_{ij}(\bm{\phi}_{ij}),\quad
e_{ij}(\bm{\phi}_{ij})=\sum_{0\leq t^{\prime}<t\leq
T}w_{t^{\prime}t}\mathbb{E}_{ij}\big{(}\Delta_{ij}\mathbb{I}\big{(}Y_{ij}^{t}\bm{\phi}_{ij}^{t^{\prime}t}(\bm{X}_{ij},Y^{t^{\prime}}_{ij})\leq
0\big{)}\big{)},$ (2)

where $e_{ij}(\bm{\phi}_{ij})$ is the pairwise misclassification error for
user $i$ on item $j$ across $T$ stages, and
$\mathbb{P}_{ij}(\cdot)=\mathbb{P}_{ij}(\cdot\mid\bm{X}_{ij}=\bm{x}_{ij})$ and
$\mathbb{E}_{ij}(\cdot)=\mathbb{E}(\cdot\mid\bm{X}_{ij}=\bm{x}_{ij})$ denote
the conditional probability and expectation given $\bm{X}_{ij}=\bm{x}_{ij}$.
In (2), $w_{t^{\prime}t}\geq 0$ is a pre-specified weight for predicting the
outcome at subsequent stage $t$ based on present stage $t^{\prime}$; $0\leq
t^{\prime}<t\leq T$, reflecting the relative importance with respect to
prediction at different stages in the overall evaluation. In particular, (2)
reduces to the next-stage prediction and last-stage prediction when
$w_{t^{\prime}t}=I(t-t^{\prime}=1)$ and $w_{t^{\prime}t}=I(t=T)$,
respectively. Moreover, if missing occurs completely at random, or missing
pattern $\Delta_{ij}$ is independent of $\bm{Y}_{ij}$, then (2) is
proportional to the standard misclassification error.

Lemma 1 gives the multistage Bayes decision function
$\bar{\bm{\phi}}_{ij}=\mathop{\rm
argmin}_{\bm{\phi}_{ij}}e_{ij}(\bm{\phi}_{ij})$ for user $i$ on item $j$ in
(2) subject to a constraint that predicted $Y_{ij}^{t}$ by the Bayes rule
satisfies the monotonic behavior chain property (1).

######  Lemma 1 (Multistage Bayes-rule)

The optimal multistage pairwise decision function ${\bar{\bm{\phi}}}_{ij}$
minimizing (2) can be written as: for $1\leq i\leq n$, $1\leq j\leq m$, $0\leq
t^{\prime}<t\leq T$,

$\displaystyle\bar{\phi}_{ij}^{t^{\prime}t}(\bm{x}_{ij},y^{t^{\prime}}_{ij})=\left\\{\begin{array}[]{ll}-1,\quad\text{if
}y^{t^{\prime}}_{ij}=-1,\\\ \mathop{\rm
sign}\big{(}\bar{f}_{ij}^{t^{\prime}t}(\bm{x}_{ij})\big{)}=\mathop{\rm
sign}\big{(}\mathbb{P}_{ij}\big{(}Y^{t}_{ij}=1|Y^{t^{\prime}}_{ij}=1,\Delta_{ij}=1\big{)}-1/2\big{)},\quad\text{if
}y^{t^{\prime}}_{ij}=1;\end{array}\right.$

The multistage Bayes rule for predicting $Y^{t}_{ij}$ by $\bm{x}_{ij}$ and
$Y^{t^{\prime}}_{ij}$ is $\mathop{\rm
sign}\big{(}\bar{\phi}_{ij}^{t^{\prime}t}(\bm{x}_{ij},Y^{t^{\prime}}_{ij})\big{)}$.
Here $\bar{f}_{ij}^{t^{\prime}t}(\bm{x}_{ij})$ satisfies the two-level
monotonic property (1) iff $Y^{t}_{ij}$ follows (1), that is, for $0\leq
t^{\prime}<t\leq T$,

Forward: $\displaystyle\mathop{\rm
sign}(\bar{f}_{ij}^{t^{\prime}(t+1)}(\bm{x}_{ij}))=-1,\text{ if }\mathop{\rm
sign}(\bar{f}_{ij}^{t^{\prime}t}(\bm{x}_{ij}))=-1,$ Backward:
$\displaystyle\mathop{\rm
sign}(\bar{f}_{ij}^{t^{\prime}t}(\bm{x}_{ij}))=-1,\text{ if }\mathop{\rm
sign}(\bar{f}_{ij}^{(t^{\prime}+1)t}(\bm{x}_{ij}))=-1.$ (4)

Moreover, there exists $\bar{h}^{r}_{ij}(\bm{x}_{ij})\geq 0$ such that
$\bar{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})$ can be written in an additive form:

$\bar{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})=\bar{h}^{0}_{ij}(\bm{x}_{ij})-\sum_{r=t^{\prime}+1}^{t}\bar{h}^{r}_{ij}(\bm{x}_{ij}),\text{
with }\bar{h}^{r}_{ij}\geq 0;\ 0\leq r\leq T.$ (5)

The two-level monotonic property (1) says that $\mathop{\rm
sign}(\bar{f}_{ij}^{t^{\prime}t}(\bm{x}_{ij}))$ is decreasing in $t$ for any
fixed $t^{\prime}$ (forward monotonicity) whereas is increasing in
$t^{\prime}$ for any fixed $t$ (backward monotonicity). Note that (5)
guarantees (1).

In view of Lemma 1, we introduce our multistage pairwise decision functions to
mimic the additive multistage Bayes rule: for some
$h_{ij}^{r}(\bm{x}_{ij})\geq 0$; $0\leq r\leq T$.

$\displaystyle\phi_{ij}^{t^{\prime}t}(\bm{x}_{ij},y_{ij}^{t^{\prime}})=\left\\{\begin{array}[]{ll}-1,\quad\text{if
}y^{t^{\prime}}_{ij}=-1,\\\
f_{ij}^{t^{\prime}t}(\bm{x}_{ij})=h_{ij}^{0}(\bm{x}_{ij})-\sum_{r=t^{\prime}+1}^{t}h_{ij}^{r}(\bm{x}_{ij}),\quad\text{if}\
y^{t^{\prime}}_{ij}=1.\end{array}\right.$ (8)

###  2.2 Two-level monotonic multistage classification

Based on the representation of decision function in (8), we rewrite
$e_{ij}(\bm{\phi}_{ij})$ in (2) as follows:

$\displaystyle e_{ij}(\mathbf{f}_{ij})$ $\displaystyle=\sum_{0\leq
t^{\prime}<t\leq
T}w_{t^{\prime}t}\mathbb{E}_{ij}\big{(}\Delta_{ij}\mathbb{I}\big{(}Y^{t}_{ij}f^{t^{\prime}t}_{ij}(\bm{X}_{ij})\leq
0\big{)}\mathbb{I}(Y^{t^{\prime}}_{ij}=1)\big{)}.$ (9)

Note that the indicator function $\mathbb{I}(\cdot)$ in (9) is difficult to
treat in optimization. Therefore, we replace it by a surrogate loss $V(u)$ for
large-margin classification, in which $V$ is a function of the corresponding
functional margin $Y^{t}_{ij}f^{t^{\prime}t}_{ij}(\bm{x}_{ij})$. This
includes, but is not limited to, the hinge loss $V(u)=(1-u)_{+}$ (Cortes and
Vapnik, 1995), the import vector machine $V(z)=\log(\exp(-z)/(1-\exp(-z)))$
(Zhu and Hastie, 2002), and the $\psi$-loss $V(u)=\min(1,(1-u)_{+})$ (Shen et
al., 2003). On this ground, we propose a multistage large-margin loss
function:

$\displaystyle L(\mathbf{f}_{ij},\bm{Z}_{ij})$ $\displaystyle=\sum_{0\leq
t^{\prime}<t\leq
T}w_{t^{\prime}t}\Delta_{ij}V\big{(}Y^{t}_{ij}f^{t^{\prime}t}_{ij}(\bm{x}_{ij})\big{)}\mathbb{I}(Y^{t^{\prime}}_{ij}=1),$
(10)

where $\bm{Z}_{ij}=(\Delta_{ij},\bm{Y}_{ij})$ and
$\mathbf{f}_{ij}=(f_{ij}^{t^{\prime}t})_{0\leq t^{{}^{\prime}}<t\leq T}$.

Lemma 2 says that a minimizer of the cost
$l(\mathbf{f})=(nm)^{-1}\sum_{ij}\mathbb{E}_{ij}L(\mathbf{f}_{ij},\bm{Z}_{ij})$
with respect to $\mathbf{f}$ satisfies the Bayes rule in Lemma 1.

######  Lemma 2 (Multistage Fisher consistency)

The minimizer of $l(\mathbf{f})$ is Fisher-consistent in that it satisfies the
Bayes rule in Lemma 1 if the surrogate loss $V(\cdot)$ is Fisher consistency
in binary classification.

Next we parametrize our decision functions based on an additive latent factor
model with to incorporate the collaborative information across users, items
and stages. In particular, we define
$h^{r}_{ij}(\bm{x}_{ij})=\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}\bm{q}_{r}$ in (5), and set
$\bm{q}_{0}=\bm{1}$ to avoid the over-parametrization. Moreover,
$a(\bm{u}_{i},\bm{s}_{i})$ and $b(\bm{v}_{j},\bm{o}_{j})$ are proposed to link
our decision function linearly to numerical predictors $\bm{u}_{i}$ and
$\bm{v}_{j}$, as well as additive latent factors structured by categorical
predictors $\bm{s}_{i}$ and $\bm{o}_{j}$. Then, the proposed prediction
function can be written as

$\displaystyle
f_{ij}^{t^{\prime}t}(\bm{x}_{ij})=a^{\top}(\bm{u}_{i},\bm{s}_{i})b(\bm{v}_{j},\bm{o}_{j})-\sum_{r=t^{\prime}+1}^{t}\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}\bm{q}_{r},$ (11) $\displaystyle\
a(\bm{u}_{i},\bm{s}_{i})=\bm{A}\bm{u}_{i}+\sum_{l=1}^{d_{1}}\bm{a}_{ls_{il}},\
b(\bm{v}_{j},\bm{o}_{j})=\bm{B}\bm{v}_{j}+\sum_{l=1}^{d_{2}}\bm{b}_{lo_{jl}},$

subject to $\bm{A}\geq\bm{0}$, $\bm{B}\geq\bm{0}$, $\bm{q}_{r}\geq\bm{0}$;
$r=1,\cdots,T$, $\bm{a}_{lh}\geq\bm{0}$; $l=1,\cdots,d_{1}$;
$h=1,\cdots,n_{l}$; $\bm{b}_{lh}\geq\bm{0}$; $l=1,\cdots,d_{2}$;
$h=1,\cdots,m_{l}$, where $\circ$ is the Hadamard product,
$\bm{A}\in\mathbb{R}^{K\times p_{1}}$ and $\bm{B}\in\mathbb{R}^{K\times
p_{2}}$ are two matrices, which transform the user-specific and item-specific
numerical features to $K$-dimensional latent vectors, and $\bm{a}_{ls_{il}}$
and $\bm{b}_{lo_{jl}}$are two $K$-dimensional latent factor for $s_{il}$ and
$o_{jl}$.

Representation (11) is highly interpretable, it leverages all feature-
interaction between users and items based on an additive model. For example,
in Deskdrop dataset, the interaction effect of “userId”-“articleId”,
“userId”-“authorId”, “userAgent”-“articleId”, “userAgent”-“authorId”, are all
captured in (11). Furthermore, nonnegative constraints for $\bm{A}\geq\bm{0}$,
$\bm{B}\geq\bm{0}$, $\bm{a}_{lh}\geq\bm{0}$ and $\bm{b}_{lh}\geq\bm{0}$ are
enforced to ensure the two-level monotonic property (1).

Figure 2: Architecture of proposed nonnegative additive latent factor models
for continuous and categorical user-item specific features as well as the
user-item-stage effect.

Given observations $(\bm{x}_{ij},y^{t}_{ij})_{(i,j)\in\Omega_{0},1\leq t\leq
T}$ and positive index set
$\Omega_{t^{\prime}}=\\{(i,j):\Delta_{ij}=1,y^{t^{\prime}}_{ij}=1\\}$ at stage
$t^{\prime}$, we propose a regularized multistage large-margin loss

$\displaystyle\min_{\mathbf{f}\in\mathcal{F}}\tilde{l}_{nm}(\mathbf{f});$
$\displaystyle\quad\tilde{l}_{nm}(\mathbf{f})=(nm)^{-1}\sum_{0\leq
t^{\prime}<t\leq
T}\sum_{(i,j)\in\Omega_{t^{\prime}}}w_{t^{\prime}t}V\big{(}y^{t}_{ij}f^{t^{\prime}t}_{ij}(\bm{x}_{ij})\big{)}+J_{\lambda}(\mathbf{f}),$
(12)

subject to $\bm{A}\geq\bm{0}$, $\bm{B}\geq\bm{0}$, $\bm{a}_{lh}\geq 0$; $1\leq
l\leq d_{1}$; $1\leq h\leq n_{l}$, $\bm{b}_{lh}\geq 0$; $1\leq l\leq
d_{2}$,$1\leq h\leq m_{l}$, where
$J_{\lambda}(\mathbf{f})=\lambda_{1}\big{(}\|\bm{A}\|_{F}^{2}+\|\bm{B}\|_{F}^{2}\big{)}+\lambda_{2}\big{(}\|\bm{a}\|_{2}^{2}+\|\bm{b}\|_{2}^{2}\big{)}+\lambda_{3}\|\bm{q}\|_{2}^{2}$,
$\bm{a}=(\bm{a}^{\top}_{1},\cdots,\bm{a}^{\top}_{d_{1}})^{\top}$,
$\bm{b}=(\bm{b}^{\top}_{1},\cdots,\bm{b}^{\top}_{d_{2}})^{\top}$,
$\bm{q}=(\bm{q}^{\top}_{1},\cdots,\bm{q}^{\top}_{T})^{\top}$,
$\bm{a}_{l}=(\bm{a}_{l1}^{\top},\cdots,\bm{a}_{ln_{l}}^{\top})^{\top}$,
$\bm{b}_{l}=(\bm{b}_{l1}^{\top},\cdots,\bm{b}_{lm_{l}}^{\top})^{\top}$,
$\lambda_{1},\lambda_{2},\lambda_{3}\geq 0$ are tuning parameters controlling
the trade-off between learning and regularization, and $\mathcal{F}$ is a
space of candidate decision functions in (11), which can be written as

$\displaystyle\mathcal{F}=\Big{\\{}\mathbf{f}=(\mathbf{f}^{t^{\prime}t}_{ij}):\mathbf{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})=a^{\top}(\bm{u}_{i},\bm{s}_{i})b(\bm{v}_{j},\bm{o}_{j})-\sum_{r=t^{\prime}+1}^{t}\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}\bm{q}_{r},$
$\displaystyle\bm{a}_{lh}\in\mathbb{R}^{K}_{+},\quad
l=1,\cdots,d_{1},h=1,\cdots,n_{l};\quad\bm{b}_{lh}\in\mathbb{R}^{K}_{+},\quad
l=1,\cdots,d_{2},h=1,\cdots,m_{l},$
$\displaystyle\bm{q}_{t}\in\mathbb{R}_{+}^{K},\quad
t=1,\cdots,T;\quad\bm{A}\in\mathbb{R}_{+}^{K\times
p_{1}},\bm{B}\in\mathbb{R}_{+}^{K\times p_{2}}\Big{\\}},$ (13)

where $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ are defined in (11).

Minimization (12) in $(\bm{A},\bm{B},\bm{q},\bm{b},\bm{q})$ yields an
estimated $(\hat{\bm{A}},\hat{\bm{B}},\hat{\bm{a}},\hat{\bm{b}},\hat{\bm{q}})$
thus $\widehat{\mathbf{f}}_{ij}^{t^{\prime}t}$ by (11) and
$\widehat{\bm{\phi}}_{ij}^{t^{\prime}t}$ by (8). Finally, the outcome of
$Y_{ij}^{t}$ is predicted by $\mathop{\rm
sign}(\widehat{\bm{\phi}}_{ij}^{t^{\prime}t}(\bm{x}_{ij},Y_{ij}^{t^{\prime}}))$;
$1\leq i\leq n$, $1\leq j\leq m$, $0\leq t^{\prime}<t\leq T$.

###  2.3 Connection with existing frameworks

This section compares the proposed method with a standard method treating each
stage individually and an ordinal method treating the stage as an ordinal
class.

Standard. A standard method treats $(t^{\prime},t)$-pairwise classification
separately and then combine the prediction results for $0\leq t^{\prime}<t\leq
T$. In fact, it estimates
$\mathbf{f}^{t^{\prime}t}=(f_{ij}^{t^{\prime}t})_{1\leq i\leq n;1\leq j\leq
m}$ by solving:

$\displaystyle\min_{\mathbf{f}^{t^{\prime}t}\in\mathcal{F}^{t^{\prime}t}}\
(nm)^{-1}\sum_{(i,j)\in\Omega_{t^{\prime}}}V\big{(}y^{t}_{ij}f^{t^{\prime}t}_{ij}(\bm{x}_{ij})\big{)}+J_{\lambda}(\mathbf{f}^{t^{\prime}t});\quad
0\leq t^{\prime}<t\leq T,$ (14)

where $\mathcal{F}^{t^{\prime}t}$ is a parameter space of candidate pairwise
decision functions. A standard method includes latent factor models (Koren,
2008), gradient boosting (Cheng et al., 2014), and deep neural networks (He et
al., 2017), and $J_{\lambda}(\mathbf{f}^{t^{\prime}t})$ regularizes
$\mathbf{f}^{t^{\prime}t}$ with a nonnegative tuning parameter $\lambda\geq
0$.

Ordinal. An ordinal method treats each stage $t^{\prime}$ separately, which
estimates
$\mathbf{f}^{t^{\prime}}=(\mathbf{f}^{t^{\prime}t})_{t^{\prime}<t\leq T}$ by
solving:

$\displaystyle\min_{\mathbf{f}^{t^{\prime}}\in\mathcal{F}^{t^{\prime}}}\
(nm)^{-1}\sum_{(i,j)\in\Omega_{t^{\prime}}}\sum_{t=t^{\prime}+1}^{T}w_{t^{\prime}t}V\big{(}y^{t}_{ij}f^{t^{\prime}t}_{ij}(\bm{x}_{ij})\big{)}+J_{\lambda}(\mathbf{f}^{t^{\prime}}),\quad\text{for
}0\leq t^{\prime}<T,$ (15)

subject to $f^{t^{\prime}(t^{\prime}+1)}_{ij}\geq
f^{t^{\prime}(t^{\prime}+2)}_{ij}\geq\cdots\geq f^{t^{\prime}T}_{ij}$, where
$\mathcal{F}^{t^{\prime}}$ is a parameter space of candidate stagewise
decision functions and $J_{\lambda}(\mathbf{f}^{t^{\prime}})$ regularizes
$\mathbf{f}^{t^{\prime}}$. For example, Crammer and Singer (2002) formulates a
parallel decision function
$f_{ij}^{t^{\prime}t}(\bm{x}_{ij})=\bm{\beta}^{\top}_{t^{\prime}}\bm{x}_{ij}+\beta_{t^{\prime}t}$
and imposes the constraint
$\beta_{t^{\prime}(t^{\prime}+1)}\geq\beta_{t^{\prime}(t^{\prime}+2)}\cdots\geq\beta_{t^{\prime}T}$.
Moreover, since $\bm{x}_{ij}$ is pre-scaled in $[0,1]$, we can set
$\mathbf{f}^{t^{\prime}t}_{ij}=(\bm{\beta}_{t^{\prime}0}-\sum_{t=t^{\prime}+1}^{T}\bm{\beta}_{t^{\prime}t})^{\top}\bm{x}_{ij}$,
where the forward monotonicity is ensured by positive constraints
$\bm{\beta}_{t^{\prime}t}\geq\bm{0}$, $t=t^{\prime}+1,\cdots,T$. Note that
when $V(y^{t}_{ij}f^{t^{\prime}t}_{ij})$ is replaced as a negative log-
likelihood function, then (15) is a formulation for ordinal regression
McCullagh and Nelder (2019); Bhaskar (2016).

Ordinal classification incorporates the forward monotonicity but ignores
backward monotonicity. Whereas such monotonicity helps to reduce the size of
the parameter space, by comparison, the proposed method can further reduce
parameters utilizing the backward monotonicity. In contrast, a standard method
does not leverage any level of monotonicity. Most critically, both standard
and ordinal methods fail to yield prediction consistency.

##  3 Large-scale computation

This section develops a computational scheme to solve nonconvex minimization
(12). For illustration, consider the hinge loss $V(u)=(1-u)_{+}$. The scheme
minimizes a nonconvex cost function (12) by solving a sequence of relaxed
convex subproblems via a block successive minimization (BSM). The scheme uses
blockwise descent to alternate the following convex subproblems. For each
subproblem, we decompose (12) into an equivalent form of many small
optimizations for parallelization and for alleviating the memory requirement.
Let $\bm{q}_{t^{\prime}t}=\bm{1}_{K}-\sum_{r=t^{\prime}+1}^{t}\bm{q}_{r}$, we
solve (12) as follows.

User-effect block $\bm{A}$. This convex optimization solves for $\bm{A}$:

$\displaystyle\min_{\text{vec}(\bm{A})\geq 0}$ $\displaystyle\
(nm)^{-1}\sum_{0\leq t^{\prime}<t\leq
T}\sum_{(i,j)\in\Omega_{t^{\prime}}}w_{t^{\prime}t}\xi_{ij}^{t^{\prime}t}+\lambda_{1}\sum_{h=1}^{n_{l}}\|\text{vec}(\bm{A})\|_{2}^{2},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq
1-y^{t}_{ij}\Big{(}\text{vec}(\bm{A})^{\top}\big{(}\bm{u}_{i}\otimes\big{(}b(\bm{v}_{j},\bm{o}_{j})\circ\bm{q}_{t^{\prime}t}\big{)}\big{)}+\big{(}\sum^{d_{1}}_{l=1}\bm{a}_{ls_{il}}\big{)}^{\top}\big{(}b(\bm{v}_{j},\bm{o}_{j})\circ\bm{q}_{t^{\prime}t}\big{)}\Big{)},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq 0;(i,j)\in\Omega_{t^{\prime}},\ 0\leq
t^{\prime}<t\leq T,$ (16)

where $\otimes$ is the Kronecker product and $\text{vec}(\bm{A})$ is the
column vectorization of matrix $\bm{A}$.

User-effect block $\bm{a}$. This convex optimization solves for
$(\bm{a}_{1},\cdots,\bm{a}_{d_{1}})$:

$\displaystyle\min_{\bm{a}_{l}\geq\bm{0}}$ $\displaystyle\
(nm)^{-1}\sum_{h=1}^{n_{l}}\sum_{0\leq t^{\prime}<t\leq
T}\sum_{\\{i:\bm{s}_{il}=h\\}}\sum_{j\in\\{j:(i,j)\in\Omega_{t^{\prime}}\\}}w_{t^{\prime}t}\xi_{ij}^{t^{\prime}t}+\lambda_{2}\sum_{h=1}^{n_{l}}\|\bm{a}_{lh}\|_{2}^{2},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq
1-y^{t}_{ij}\Big{(}\bm{a}^{\top}_{ls_{il}}\big{(}b(\bm{v}_{j},\bm{o}_{j})\circ\bm{q}_{t^{\prime}t}\big{)}+\big{(}\bm{A}\bm{u}_{i}+\sum_{l^{\prime}\neq
l}\bm{a}_{l^{\prime}s_{il^{\prime}}}\big{)}^{\top}\big{(}b(\bm{v}_{j},\bm{o}_{j})\circ\bm{q}_{t^{\prime}t}\big{)}\Big{)},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq 0;(i,j)\in\Omega_{t^{\prime}};\ 0\leq
t^{\prime}<t\leq T.$ (17)

Note that (3) can be separately solved for each $\bm{a}_{lh}$ in a parallel
fashion,

$\displaystyle\min_{\bm{a}_{lh}\geq\bm{0}}$ $\displaystyle\
(nm)^{-1}\sum_{0\leq t^{\prime}<t\leq
T}\sum_{\\{i:\bm{s}_{il}=h\\}}\sum_{j\in\\{j:(i,j)\in\Omega_{t^{\prime}}\\}}w_{t^{\prime}t}\xi_{ij}^{t^{\prime}t}+\lambda_{2}\|\bm{a}_{lh}\|_{2}^{2},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq
1-y^{t}_{ij}\Big{(}\bm{a}^{\top}_{lh}\big{(}b(\bm{v}_{j},\bm{o}_{j})\circ\bm{q}_{t^{\prime}t}\big{)}+\big{(}\bm{A}\bm{u}_{i}+\sum_{l^{\prime}\neq
l}\bm{a}_{l^{\prime}s_{il^{\prime}}}\big{)}^{\top}\big{(}b(\bm{v}_{j},\bm{o}_{j})\circ\bm{q}_{t^{\prime}t}\big{)}\Big{)},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq 0;\ (i,j)\in\Omega_{t^{\prime}};\
0\leq t^{\prime}<t\leq T.$ (18)

Item-effect block $\bm{B}$. This convex optimization solves for $\bm{B}$:

$\displaystyle\min_{\text{vec}(\bm{B})\geq 0}$ $\displaystyle\
(nm)^{-1}\sum_{0\leq t^{\prime}<t\leq
T}\sum_{(i,j)\in\Omega_{t^{\prime}}}w_{t^{\prime}t}\xi_{ij}^{t^{\prime}t}+\lambda_{1}\sum_{h=1}^{n_{l}}\|\text{vec}(\bm{B})\|_{2}^{2},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq
1-y^{t}_{ij}\Big{(}\text{vec}(\bm{B})^{\top}\big{(}\bm{v}_{j}\otimes\big{(}a(\bm{u}_{i},\bm{s}_{j})\circ\bm{q}_{t^{\prime}t}\big{)}\big{)}+\big{(}\sum^{d_{1}}_{l=1}\bm{b}_{ls_{il}}\big{)}^{\top}\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ\bm{q}_{t^{\prime}t}\big{)}\Big{)},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq 0;(i,j)\in\Omega_{t^{\prime}},\ 0\leq
t^{\prime}<t\leq T.$ (19)

Item-effect block $\bm{b}$. For $l=1,\cdots,d_{2}$,
$(\bm{b}_{lh})_{h=1}^{m_{l}}$ is solved in a parallel fashion,

$\displaystyle\min_{\bm{b}_{lh}\geq\bm{0}}$ $\displaystyle\
(nm)^{-1}\sum_{0\leq t^{\prime}<t\leq
T}\sum_{\\{j:\bm{o}_{jl}=h\\}}\sum_{i\in\\{i:(i,j)\in\Omega_{t^{\prime}}\\}}w_{t^{\prime}t}\xi_{ij}^{t^{\prime}t}+\lambda_{2}\|\bm{b}_{lh}\|_{2}^{2},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq
1-y^{t}_{ij}\Big{(}\bm{b}^{\top}_{lh}\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ\bm{q}_{t^{\prime}t}\big{)}+\big{(}\bm{B}\bm{v}_{j}+\sum_{l^{\prime}\neq
l}\bm{b}_{l^{\prime}o_{jl^{\prime}}}\big{)}^{\top}\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ\bm{q}_{t^{\prime}t}\big{)}\Big{)},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq 0;j\in\\{j:\bm{o}_{jl}=h\\},\
(i,j)\in\Omega_{t^{\prime}};\ 0\leq t^{\prime}<t\leq T.$ (20)

Stage-effect block $\bm{q}$. This convex optimization solves for $\bm{q}_{r}$,
for $r=1,\cdots,T$, given the present values for all other variables:

$\displaystyle\min_{\bm{q}_{r}\geq\bm{0}}$ $\displaystyle\
(nm)^{-1}\sum_{0\leq t^{\prime}<r\leq t\leq
T}\sum_{(i,j)\in\Omega_{t^{\prime}}}w_{t^{\prime}t}\xi_{ij}^{t^{\prime}t}+\lambda_{3}\|\bm{q}_{r}\|_{2}^{2},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq
1-y^{t}_{ij}\Big{(}-\bm{q}_{r}^{\top}\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{)}+(\bm{1}-\sum_{r^{\prime}=t^{\prime}+1;r^{\prime}\neq
r}^{t}\bm{q}_{r})^{\top}\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{)}\Big{)},$
$\displaystyle\xi_{ij}^{t^{\prime}t}\geq 0;(i,j)\in\Omega_{t^{\prime}};\ 0\leq
t^{\prime}<t\leq T.$ (21)

As a technical note, (3)-(3) are standard SVMs with fixed intercept and
positive constrains for model parameters, which can be solved via one
efficient implementation in our python package
varsvm333https://pypi.org/project/varsvm/ based on the coordinate decent
algorithm (Wright, 2015).

The aforementioned scheme is summarized in Algorithm 1.

Algorithm 1 (Parallelized version: hinge loss)

Step 1 (Initialization). Initialize the values of ($\bm{A}$, $\bm{B}$,
$\bm{a}$, $\bm{b}$, $\bm{q}$) and specify the tolerance error.

Step 2 (Update $\bm{A}$ and $\bm{a}$). Update $\bm{A}$ by solving (3) given
present values of the other variables. For $l=1,\cdots,d_{1}$, update
$\bm{a}_{lh}$; $h=1,\cdots,n_{l}$ in a parallel fashion by solving (3) given
present values of the other variables.

Step 3 (Update $\bm{B}$ and $\bm{b}$). Update $\bm{B}$ by solving (3) given
present values of the other variables. For $l=1,\cdots,d_{2}$, update
$\bm{b}_{lh}$; $h=1,\cdots,m_{l}$ in a parallel fashion by solving (3) given
present values of the other variables.

Step 4 (Update $\bm{q}$). For $r=1,\cdots,T$, update $\bm{q}_{r}$ by solving
(3) given present values of the other variables.

Step 5. Iterate Steps 2-4 until the decrement of the cost function in (12) is
less than the tolerance error.

Algorithm 1 returns an estimate
$(\widehat{\bm{A}},\widehat{\bm{B}},\widehat{\bm{a}},\widehat{\bm{b}},\widehat{\bm{q}})$
at termination, which in turn yields an estimated decision functions
$\widehat{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})=\widehat{a}^{\top}(\bm{u}_{i},\bm{s}_{i})\widehat{b}(\bm{v}_{j},\bm{o}_{j})-\sum_{r=t^{\prime}+1}^{t}\big{(}\widehat{a}(\bm{u}_{i},\bm{s}_{i})\circ\widehat{b}(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}\widehat{\bm{q}}_{r}$,
and $\widehat{\phi}^{t^{\prime}t}_{ij}=-1$ if $Y^{t^{\prime}}_{ij}=-1$,
$\widehat{f}^{t^{\prime}t}_{ij}$ otherwise. Finally, the prediction at the
$t$-th stage given the $t^{\prime}$-th stage is $\mathop{\rm
sign}\big{(}\widehat{\phi}^{t^{\prime}t}_{ij}(y^{t^{\prime}}_{ij},\bm{x}_{ij})\big{)}$
for new user-item covariates $\bm{x}_{ij}$.

To implement Algorithm 1, we develop a Python library to solve an SVM-type of
optimization, including weighted SVMs, drifted SVMs, and non-negative SVMs,
based on coordinate descent of the dual problem (Wright, 2015). This
implementation has been available in varsvm library in Python on the GitHub
repository varsvm444https://github.com/statmlben/variant-svm, and the source
code is released with the MIT license using GitHub and the release is
available via PyPI. The development undergoes an integration of Ubuntu Linux
and Mac OS X, in Python 3.

######  Lemma 3 (Convergence of Algorithm 1)

A solution
$(\widehat{\bm{A}},\widehat{\bm{a}},\widehat{\bm{B}},\widehat{\bm{b}},\widehat{\bm{q}})$
of Algorithm 1 is a stationary point of $\tilde{l}_{nm}$ in (12) in that

$\displaystyle\tilde{\bm{A}}$ $\displaystyle=\mathop{\rm
argmin}_{\bm{A}}\tilde{l}_{nm}(\bm{A},\tilde{\bm{a}},\tilde{\bm{B}},\tilde{\bm{b}},\tilde{\bm{q}}),\quad\tilde{\bm{B}}=\mathop{\rm
argmin}_{\bm{B}}\tilde{l}_{nm}(\tilde{\bm{A}},\tilde{\bm{a}},\bm{B},\tilde{\bm{b}},\tilde{\bm{q}}),$
$\displaystyle\tilde{\bm{a}}_{l}$ $\displaystyle=\mathop{\rm
argmin}_{\bm{a}_{l}}\tilde{l}_{nm}(\tilde{\bm{A}},\tilde{\bm{a}}_{1},\cdots,\bm{a}_{l},\cdots,\tilde{\bm{a}}_{d_{1}},\tilde{\bm{B}},\tilde{\bm{b}},\tilde{\bm{q}});\
1\leq l\leq d_{1};$ $\displaystyle\tilde{\bm{b}}_{l}$
$\displaystyle=\mathop{\rm
argmin}_{\bm{b}_{l}}\tilde{l}_{nm}(\tilde{\bm{A}},\tilde{\bm{a}},\tilde{\bm{B}},\tilde{\bm{b}}_{1},\cdots,\bm{b}_{l},\cdots,\tilde{\bm{b}}_{d_{2}},\tilde{\bm{q}});\
1\leq l\leq d_{2};$ $\displaystyle\tilde{\bm{q}}_{t}$
$\displaystyle=\mathop{\rm
argmin}_{\bm{q}_{t}}\tilde{l}_{nm}(\tilde{\bm{A}},\tilde{\bm{a}},\tilde{\bm{B}},\tilde{\bm{b}},\tilde{\bm{q}}_{1},\cdots,\bm{q}_{t},\cdots,\tilde{\bm{q}}_{T});\
1\leq t\leq T.$

The successive update in Algorithm 1 suffices to ensure its convergence
because the cost function in (12) is strictly blockwise convex, that is, (3),
(3), (3), (3), and (3) are strictly convex in terms of its parameters within
each block, yielding a unique minimizer at each step. This aspect differs from
the maximum block improvement Chen et al. (2012) that guarantees the
convergence of blockwise coordinate descent for a general objective function.
One benefit of successive updating is that it significantly reduces
computational complexity. Note that the solution of Algorithm 1 can be a
global minimizer when certain additional assumptions are made Haeffele and
Vidal (2015).

##  4 Theory

This section investigates the generalization aspect of the proposed multistage
recommender $\hat{\mathbf{f}}$ in terms of the accuracy of classification, as
measured by the classification regret, defined as
$e(\hat{\mathbf{f}})-e(\bar{\mathbf{f}})$, where $e(\cdot)$ is the
generalization error defined in (2).

Assume that $(\Delta_{ij},Y_{ij})$ given $\bm{X}_{ij}$; $1\leq i\leq n$;
$1\leq j\leq m$, are conditionally independent, although
$(\Delta_{ij},\bm{X}_{ij},Y_{ij})$’s may not be independent as multiple items
can be purchased for a same user or multiple users may purchase a same item.
Let $\bar{\mathbf{f}}\in\mathcal{F}$ be a Bayes decision function, with
$(\bar{\bm{A}},\bar{\bm{B}},\bar{\bm{a}},\bar{\bm{b}},\bar{\bm{q}})$, and the
number of latent factors be $\bar{K}$. Let the truncated $V$-loss be
$V_{B}(u)=\min(V(u),B)$, where $B\geq
V\big{(}Y^{t}_{ij}\bar{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})\big{)}$ for any
$i=1,\cdots,n;j=1,\cdots,m;0\leq t^{\prime}<t\leq T$, and the corresponding
loss defined by $V_{B}(u)$ be

$\displaystyle L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})=\sum_{0\leq t^{\prime}<t\leq
T}w_{t^{\prime}t}\Delta_{ij}V_{B}\big{(}Y^{t}_{ij}f^{t^{\prime}t}_{ij}(\bm{X}_{ij})\big{)}\mathbb{I}(Y^{t^{\prime}}_{ij}=1),$

and the truncated cost function is
$l_{B}(\mathbf{f})=(nm)^{-1}\sum_{i,j}\mathbb{E}_{ij}L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})$.
The following technical assumptions are assumed:

Assumption A (Conversion). There exist constants $0\leq\alpha\leq\infty$ and
$c_{1}>0$ such that for all $0<\epsilon\leq B$ and $\mathbf{f}\in\mathcal{F}$,

$\sup_{\\{\mathbf{f}\in\mathcal{F};\
l_{B}(\mathbf{f})-l_{B}(\bar{\mathbf{f}})\leq\epsilon\\}}e(\mathbf{f})-e(\bar{\mathbf{f}})\leq
c_{1}\epsilon^{\alpha}.$

Assumption B (Variance property). There exist constants $0\leq\mu\leq 2$ and
$c_{2}>0$ such that for any $\epsilon\geq 0$, $n,m\geq 1$, and
$\mathbf{f}\in\mathcal{F}$,

$\sup_{\\{\mathbf{f}\in\mathcal{F}:\
l_{B}(\mathbf{f})-l_{B}(\bar{\mathbf{f}})\leq\epsilon\\}}(nm)^{-1}\sum_{i=1}^{n}\sum_{j=1}^{m}\mathbb{E}_{ij}\Big{(}\big{(}L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})-L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})\big{)}^{2}\Big{)}\leq
c_{2}\epsilon^{\mu}.$

Assumptions A and B are two local smooth conditions regarding a connection
between two different metrics and the relation between the mean and variance
of the regret function in a neighborhood of $\bar{\mathbf{f}}$. Similar
assumptions have been used in binary classification can be founded in Zhang
and Liu (2014), for example, $\alpha=\mu=1$ when $V(\cdot)$ is the $\psi$-loss
Shen et al. (2003) or the hinge loss under low-noise assumption.

Let $K$ and $\bar{K}$ be the numbers of latent factor of $\hat{\mathbf{f}}$
and $\bar{\mathbf{f}}$, respectively.

###### Theorem 4

Let $\hat{\mathbf{f}}$ be a global minimizer of (12) with $K\geq\bar{K}$. If
Assumptions A and B hold, then there exist constants $c_{3}>0$ and $c_{4}>0$
such that

$\mathbb{P}\big{(}e(\hat{\mathbf{f}})-e(\bar{\mathbf{f}})\geq\epsilon^{2\alpha}_{nm}\big{)}\leq
c_{3}\exp\big{(}-c_{4}nm(\lambda\bar{J})^{2-\min(1,\mu)}\big{)},$

where $\bar{J}=\max(J(\bar{\mathbf{f}}),1)$,
$J_{\lambda}(\bar{\mathbf{f}})\sim\epsilon^{2}_{nm}/2$, and

$\epsilon^{2}_{nm}=O\Big{\\{}\Big{(}\frac{(\Lambda+T)K}{nm}\log\big{(}\frac{nm\bar{w}B\bar{J}(p_{1}+p_{2}+d_{1}+d_{2})(6T^{3}+T^{2}K)}{(\Lambda+T)K}\big{)}\Big{)}^{\frac{1}{2-\min(1,\mu)}}\Big{\\}},$
(22)

where $\bar{w}=\sum_{0\leq t^{\prime}<t\leq T}w^{2}_{t^{\prime}t}$,
$(\Lambda+T)K$ is the total number of model parameters, and
$\Lambda=p_{1}+p_{2}+\sum_{l=1}^{d_{1}}n_{l}+\sum_{l=1}^{d_{2}}m_{l}$.

When $nm\geq(p_{1}+p_{2}+d_{1}+d_{2})(\Lambda+T)K(6T^{3}+T^{2}K)$, (22)
suggests that the convergence rate for the proposed method is
$O\big{(}\big{(}\frac{(\Lambda+T)K}{nm}\log(\frac{nm}{(\Lambda+T)K})\big{)}^{\frac{\alpha}{2-\min(1,\mu)}}\big{)}$,
which becomes
$O\big{(}\frac{(\Lambda+T)K}{nm}\log(\frac{nm}{(\Lambda+T)K})\big{)}$ when
$\alpha=\mu=1$, and matches with the minimax rate of the standard
classification in (Tsybakov et al., 2004).

As indicated in (22), the convergence rate is improved by the two-level
monotonic chain property in (1) and the additive property in (5). As
illustrated in Table 1, based on the factorization model in (11) while
ignoring the two-level monotonicity (1), a standard method and an ordinal
method involve $\Lambda KT(T+1)/2$ and $\Lambda KT$ parameters, which are more
than $(\Lambda+T)K$ and significantly impedes the performance when $T$ becomes
large. Furthermore, due to the monotonic property in (1), the effective sample
size for $Y^{t^{\prime}}_{ij}=1$ decreases exponentially as $t^{\prime}$
increases. Then the standard and ordinal methods suffer from over-fitting at
late stages with a small sample size. In contrast, the proposed method
leverages the two-level monotonicity to reduce the dimension of the underlying
problem.

Method |  | Forward monotonicity |  | Backward monotonicity |  | Consistency |  | #Parameters |
---|---|---|---|---|---|---|---|---|---
Proposed |  | ✓ |  | ✓ |  | ✓ |  | $(\Lambda+T)K$ |
Standard |  | ✗ |  | ✗ |  | ✗ |  | $\Lambda KT(T+1)/2$ |
Ordinal |  | ✓ |  | ✗ |  | ✗ |  | $\Lambda KT$ |

Table 1: Monotonicity and model parameters of the proposed, standard, and
ordinal frameworks, denoting (12), (14) and (15). Here forward and backward
monotonicity is defined in (1), consistency denotes if the methods can provide
a consistent prediction result, and #Parameters denotes the number of model
parameters based on (5).

##  5 Numerical examples

This section examines the proposed method in (12) with the hinge loss
$V(u)=(1-u)_{+}$ and compares it with standard methods, including latent
factor models (SVD++; (Koren, 2008)), gradient boosting (GradBoost) (Friedman,
2001), support vector machine (SVM) (Cortes and Vapnik, 1995), deep neural
network (DeepNN) (Schmidhuber, 2015), ordinal method (15), namely ordinal
support vector machine (OSVM), in a simulation and an article sharing
benchmark.

For a standard method, SVD++, GradBoost, SVM, and DeepNN treat each
$(t^{\prime},t)$ pair separately for all possible pairs, where we use one-hot
encoder (Weinberger et al., 2009) to convert a categorical covariate to 0-1
numerical predictors for training. Moreover, for SVD++, the number of latent
factors is set as 20, and latent factors are estimated using the first columns
of $\bm{s}_{i}$ and $\bm{o}_{j}$ with the rating
$(Y^{t}_{ij}+1)/2\in\\{0,1\\}$ while the label is predicted as 1 if the
predicted rating exceeds 0.5 and -1 otherwise; for GradBoost, the minimum
number of samples to split an internal node is set to be 2 and the minimum
number of samples required to be at a leaf node is 1, and the number of
boosting stages is set to be 10. For DeepNN, a ReLU network is used with the
number of node in each layer being fixed at 32. For the ordinal method, we set
the OSVM classifier in (15) as
$\mathbf{f}^{t^{\prime}t}_{ij}=(\bm{\beta}_{t^{\prime}0}-\sum_{t_{0}=t^{\prime}+1}^{t}\bm{\beta}_{t^{\prime}t_{0}})^{\top}\bm{x}_{ij}$
since $\bm{x}_{ij}$ is pre-scaled in $[0,1]$, where the forward monotonicity
is ensured by positivity constraints $\bm{\beta}_{t^{\prime}t}\geq\bm{0}$,
$t=t^{\prime}+1,\cdots,T$.

For implementation, we use our Python library varsvm for the proposed method
and OSVM, the Python library sklearn555https://scikit-learn.org for GradBoost,
SVM, and DeepNN, and the Python library Surprise666http://surpriselib.com/ for
SVD++. For all methods, the prediction of $Y_{ij}^{t}$ is automatically set as
$-1$ when $Y^{t^{\prime}}_{ij}=-1$.

###  5.1 Simulations

In simulations, data
$\\{\bm{s}_{i},\bm{o}_{j},y^{t}_{ij}\\}_{i=1,\cdots,n;j=1,\cdots,m;t=0,\cdots,T}$
are generated to mimic a multistage recommender system on a monotonic chain,
with $d_{1}=3$, $d_{2}=2$, $n=n_{1}=50$, $n_{2}=30$, $n_{3}=50$,
$m=m_{1}=100$, and $m_{2}=40$. Then we set $K=20$, and for each $(l,h,t)$,
$\bm{a}_{lh}\sim\chi^{2}(1)$, $\bm{b}_{lh}\sim\chi^{2}(1)$, and
$\bm{q}_{t}\sim\chi^{2}(1)$ are independently sampled from the chi-squared
distribution on 1 degree of freedom, and $(s_{il})_{l=1}^{3}$ and
$(o_{jl})_{l=1}^{2}$ are uniformly sampled from $\\{1,\cdots,n_{l}\\}$ and
$\\{1,\dots,m_{l}\\}$, respectively. Moreover, we set $y^{0}_{ij}=1$ at stage
0, set
$p^{1}_{ij}=\bm{a}^{\top}(\bm{s}_{i})\bm{b}(\bm{o}_{j})-\big{(}\bm{a}(\bm{s}_{i})\circ\bm{b}(\bm{o}_{j})\big{)}^{\top}\bm{q}_{1}$,
$y^{1}_{ij}=\mathop{\rm sign}\big{(}p^{1}_{ij}+\sigma_{1}N(0,.1)\big{)}$ at
stage 1, and set
$p^{2}_{ij}=\bm{a}^{\top}(\bm{s}_{i})\bm{b}(\bm{t}_{j})-\big{(}\bm{a}(\bm{s}_{i})\circ\bm{b}(\bm{t}_{j})\big{)}^{\top}\bm{q}_{2}$,
$y^{2}_{ij}=\mathop{\rm sign}\big{(}p^{2}_{ij}+\sigma_{2}N(0,.1)\big{)}$ for
$y^{1}_{ij}=1$, and $y^{2}_{ij}=-1$ for $y^{1}_{ij}=-1$ at stage 2 to satisfy
the monotonic chain property (1), where $\sigma_{t}$ is the standard deviation
for $(p^{t}_{ij})_{(i,j)\in\Omega_{0}}$; $t=1,2$. To generate dataset, we
choose indices $(i,j)$ of $\Omega_{0}$ randomly with $|\Omega_{0}|=50000$, and
the resulting missing ratio is 99.98%, and $\Omega_{1}$ and $\Omega_{2}$ are
generated as $\Omega_{1}=\\{(i,j)\in\Omega_{0}:y^{1}_{ij}=1\\}$,
$\Omega_{2}=\\{(i,j)\in\Omega_{0}:y^{2}_{ij}=1\\}$. This leads to
$\\{\bm{s}_{i},\bm{o}_{j},\bm{y}_{ij}\\}_{(i,j)\in\Omega_{0}}$. Finally, the
data is split into training, validation, and testing sets with a partition
ratio of 10%, 10%, and 80%, respectively.

To avoid overfitting, for gradient boosting, the percentage of features are
tuned via a grid search based on the overall misclassification error on the
validation set, where the grid set is chosen as
$\\{.001,.005,.2,.4,.6,.8,1.\\}$. For the neural network, its depth is tuned
over a set $\\{1,2,3,5,7,9,11\\}$. For SVM, OSVM, and SVD++, the penalization
parameters are tuned based on validation set through
$\\{.001,.005,.01,.05,.1,.5,1,5,10,50,100,500\\}$. For the proposed method,
$K=20$ and three tuning parameters $(\lambda_{1},\lambda_{2},\lambda_{3})$ are
tuned by minimizing the misclassification error on a validation set via a grid
search over $\lambda_{1}\in\\{.001,.005,.01\\}$,
$\lambda_{2}\in\\{.01,.05,.1\\}$, and $\lambda_{3}\in\\{.0001,.0005,.001\\}$.
For evaluation, we set all stage weights $w_{t^{\prime}t}=1$.

As suggested by Table 2, the proposed multistage recommender performs the best
for one-step, two-step, and three-step forecasting, followed by DeepNN, SVM,
and GradBoost. In terms of the overall performance, this same phenomenon is
also observed and the amount of improvements of the proposed method are
$43\%$, 22%, $20\%$, and $13\%$ over Gradient Boosting, OSVM, SVM, and DeepNN.
Concerning inconsistency prediction, SVM and OSVM perform the worst, followed
by DeepNN and Gradient Boosting. Note that there is no single instance in
which prediction by the proposed method is inconsistent, which is ensured by
the two-level monotonicity property of the classier.

(observe $t^{\prime}$, predict $t$) |  | Proposed |  | SVD++ |  | GradBoost |  | SVM |  | DeepNN |  | OSVM
---|---|---|---|---|---|---|---|---|---|---|---|---
(Stage 0, Stage 1) |  | .043(.000) |  | .095(.001) |  | .062(.000) |  | .044(.000) |  | .044(.000) |  | .059(.000)
(Stage 0, Stage 2) |  | .120(.001) |  | .387(.002) |  | .217(.000) |  | .153(.000) |  | .146(.000) |  | .165(.000)
(Stage 0, Stage 3) |  | .148(.000) |  | .417(.003) |  | .272(.001) |  | .187(.000) |  | .185(.000) |  | .199(.000)
(Stage 1, Stage 2) |  | .097(.001) |  | .294(.002) |  | .158(.000) |  | .103(.000) |  | .100(.000) |  | .116(.000)
(Stage 1, Stage 3) |  | .133(.000) |  | .414(.002) |  | .246(.001) |  | .163(.000) |  | .160(.000) |  | .175(.000)
(Stage 2, Stage 3) |  | .095(.001) |  | .341(.003) |  | .157(.001) |  | .101(.000) |  | .099(.000) |  | .102(.000)
Overall |  | .106(.000) |  | .325(.001) |  | .185(.000) |  | .125(.000) |  | .122(.000) |  | .136(.000)
%Inconsist |  | .000%(.000) |  | .680%(.001) |  | 2.22%(.000) |  | 10.9%(.000) |  | 6.01%(.000) |  | 10.9%(.001)

Table 2:  Test errors of six competitors and their estimated standard
deviations in parentheses on simulated examples over 50 replications as well
as a proportion of inconsistent instances violating monotonicity (1),denoted
by %Inconsist. Here proposed, SVD++, GradBoost, SVM, DeepNN and OSVM denote
the proposed method in (12) with the hinge loss, the latent factor model
(Koren, 2008), gradient boosting (Friedman, 2001), support vector machine
(Cortes and Vapnik, 1995), deep learner (Schmidhuber, 2015), and conditional
ordinal support vector machine in (15). The best performer in each case is
bold-faced.

###  5.2 Benchmark: Deskdrop article sharing

A Deskdrop article sharing dataset contains a sample of 12 months logs and
about 73,000 logged user interactions on more than 3,000 public articles. In
particular, article-specific features are “article Id”, “author Id”, “plain
text”, as well as “language”, and user-specific features include: “user Id”,
“user agent”, and “user region”, and users’ actions are logged, including
view, like, and follow. The goal is to predict user-item interaction over
three stage pairs.

For the proposed method, $\bm{v}_{j}$ is the numerical embedding based on
Doc2Vec (Le and Mikolov, 2014) of the “plain text” of item $j$, $\bm{o}_{j}$
consists of “article Id”, “author Id”, and “language” of item $j$, and
$\bm{s}_{i}$ is composed of “user Id”, “user agent”, and “user region” of user
$i$. For other competing methods, we use all features based on a one-hot
encoder to convert categorical covariates to zero-one dummy predictors.
Moreover, we predict $(t^{\prime}=0,t=1)$ based on the proposed method with
$T=1$, and predict $(t^{\prime}=0,t=2)$ and $(t^{\prime}=1,t=2)$ based on the
proposed method with $T=2$.

For evaluation, we set all stage weights $w_{t^{\prime}t}=1$ and adjust the
class weights at each stage inversely proportional to the number of
observations in each class since feedback at each stage is highly imbalanced
for this data. Hence, all methods are fitted based on balanced class weights
or over- or under-sampling, except SVD++ which based on a regression approach.
Yet, for consistency, we focus on balanced zero-one loss for comparison since
all methods are fitted and tuned based on a balanced-weighted classification
loss.

The misclassification errors of all competing methods for each
$(t^{\prime},t)$ pair are computed based on a nested 15-fold cross-validation
with an outer 5-fold loop and an inner 3-fold loop (Cawley and Talbot, 2010).
For the proposed method, $K=7$ and three tuning parameters
$(\lambda_{1},\lambda_{2},\lambda_{3})$ are tuned by minimizing the
misclassification error on a validation set via a grid search over
$\lambda_{1}\in\\{.3,.5,.7,1\\}$, $\lambda_{2}\in\\{2.,2.5,3.\\}$, and
$\lambda_{3}\in\\{.001,.005\\}$. For the neural network, SVM, OSVM, and SVD++,
the tuning parameters are selected with the same manner as simulation.

As suggested by Table 3, the proposed multistage recommender outperforms
SVD++, GradBoost, SVM, DeepNN, and OSVM, in terms of the overall performance
with the amount of improvement ranging from 75.3% to 8.20%, in terms of the
class-balanced zero-one loss. Similarly, for one-step and two-step forward
prediction, it is either the best or nearly the best. Concerning inconsistency
prediction, all other methods have inconsistent cases ranging from 13.7% to
0.98%. Interestingly, the conditional ordinal method (OSVM) performs similarly
to the conditional individual method (SVM), indicating that it does not fully
account for the monotonicity property.

(observe $t^{\prime}$, predict $t$) |  | Proposed |  | SVD++ |  | GradBoost |  | SVM |  | DeepNN |  | OSVM
---|---|---|---|---|---|---|---|---|---|---|---|---
(Stage 0, Stage 1) |  | 0.136(.001) |  | 0.499(.003) |  | 0.166(.000) |  | 0.139(.000) |  | 0.154(.001) |  | 0.140(.000)
(Stage 0, Stage 2) |  | 0.191(.001) |  | 0.499(.007) |  | 0.218(.001) |  | 0.218(.001) |  | 0.223(.006) |  | 0.216(.001)
(Stage 1, Stage 2) |  | 0.043(.000) |  | 0.497(.007) |  | 0.049(.000) |  | 0.045(.000) |  | 0.139(.011) |  | 0.045(.000)
Overall |  | 0.123(.001) |  | 0.498(.006) |  | 0.144(.000) |  | 0.134(.000) |  | 0.172(.005) |  | 0.134(.000)
%Inconsist |  | 0.00%(.000) |  | 7.79%(.002) |  | 1.40%(.000) |  | 1.90%(.000) |  | 13.7%(.005) |  | 0.98%(.002)

Table 3: Class-balanced zero-one losses (CB-01) of six competitors and their
estimated standard deviations in parentheses on the Deskdrop benchmark, in
addition to %Inconsist denoting a proportion of inconsistent instances
violating the monotonicity (1). Here proposed, SVD++, GradBoost, SVM, DeepNN
and OSVM denote the proposed method in (12) with the hinge loss, the latent
factor model (Koren, 2008), gradient boosting (Friedman, 2001), support vector
machine (Cortes and Vapnik, 1995), deep learner (Schmidhuber, 2015), and
ordinal support vector machine (15). The best performance in each case is
bold-faced.

## Technical Proofs

Proof of Lemma 1. It suffices to show that
$\bar{\bm{\phi}}^{t^{\prime}t}_{ij}$ is a global minimizer of
$l_{ij}(\mathbf{f}_{ij})$. By (1),

$\displaystyle\mathbb{E}_{ij}\big{(}\mathbb{I}(Y^{t}_{ij}{\phi}^{t^{\prime}t}_{ij}(\bm{x}_{ij},Y^{t^{\prime}}_{ij})\leq
0)|\Delta_{ij}=1\big{)}$
$\displaystyle=\mathbb{P}_{ij}(Y^{t^{\prime}}_{ij}=1|\Delta_{ij}=1)\mathbb{E}_{ij}\big{(}\mathbb{I}(Y^{t}_{ij}\phi^{t^{\prime}t}_{ij}(\bm{x}_{ij},1)\leq
0)|\Delta_{ij}=1,Y^{t^{\prime}}_{ij}=1\big{)}$
$\displaystyle\quad+\mathbb{P}_{ij}\big{(}Y^{t^{\prime}}_{ij}=-1|\Delta_{ij}=1\big{)}\mathbb{E}_{ij}\big{(}\mathbb{I}(-\phi^{t^{\prime}t}_{ij}(\bm{x}_{ij},-1)\leq
0)|\Delta_{ij}=1,Y^{t^{\prime}}_{ij}=-1\big{)},$ (23)

which yields that $\bar{\phi}^{t^{\prime}t}_{ij}(\bm{x}_{ij},y^{t}_{ij})$ is a
global minimizer of (Technical Proofs), since
$\bar{\phi}^{t^{\prime}t}_{ij}(\bm{x}_{ij},1)=\mathop{\rm
sign}\big{(}\mathbb{P}_{ij}(Y^{t}_{ij}=1|Y^{t^{\prime}}_{ij}=1,\Delta_{ij}=1)-1/2\big{)}$
minimizes
$\mathbb{E}_{ij}\big{(}\mathbb{I}(Y^{t}_{ij}{\phi}^{t^{\prime}t}_{ij}(\bm{x}_{ij},1)\leq
0)|Y^{t^{\prime}}_{ij}=1,\Delta_{ij}=1\big{)}$ and
$\bar{\phi}^{t^{\prime}t}_{ij}(\bm{x}_{ij},-1)=-1$ minimizes
$\mathbb{E}_{ij}\big{(}\mathbb{I}(-\phi_{ij}^{t^{\prime}t}(\bm{x}_{ij},-1)\leq
0)|\Delta_{ij}=1,Y^{t^{\prime}}_{ij}=-1\big{)}$.

Next, we verify the two-level monotonic property (1) of
$\bar{f}^{t^{\prime}t}_{ij}$. Note that

$\displaystyle\mathbb{P}_{ij}(Y^{t}_{ij}=1|\Delta_{ij}=1)$
$\displaystyle=\mathbb{P}_{ij}(Y^{t}_{ij}=1,Y^{t-1}_{ij}=1|\Delta_{ij}=1)+\mathbb{P}_{ij}(Y^{t}_{ij}=1,Y^{t-1}_{ij}=-1|\Delta_{ij}=1)$
$\displaystyle=\mathbb{P}_{ij}(Y^{t}_{ij}=1|Y^{t-1}_{ij}=1,\Delta_{ij}=1)\mathbb{P}_{ij}(Y^{t-1}_{ij}=1|\Delta_{ij}=1)$
$\displaystyle\leq\mathbb{P}_{ij}(Y^{t-1}_{ij}=1|\Delta_{ij}=1),$ (24)

which implies that $\mathbb{P}_{ij}(Y^{t}_{ij}=1|\Delta_{ij}=1)$ is
nondecreasing over stage $t$. Moreover,

$\displaystyle\mathbb{P}_{ij}(Y^{t}_{ij}=1|Y^{t^{\prime}}_{ij}=1,\Delta_{ij}=1)=\frac{\mathbb{P}_{ij}(Y^{t}_{ij}=1,Y^{t^{\prime}}_{ij}=1|\Delta_{ij}=1)}{\mathbb{P}_{ij}(Y^{t^{\prime}}_{ij}=1|\Delta_{ij}=1)}=\frac{\mathbb{P}_{ij}(Y^{t}_{ij}=1|\Delta_{ij}=1)}{\mathbb{P}_{ij}(Y^{t^{\prime}}_{ij}=1|\Delta_{ij}=1)},$

provided that $\mathbb{P}_{ij}(Y^{t^{\prime}}_{ij}=1|\Delta_{ij}=1)\neq 0$,
which is non-decreasing when $t$ increases while $t^{\prime}$ is fixed because
$\mathbb{P}_{ij}(Y^{t}_{ij}=1|\Delta_{ij}=1)$ is non-decreasing in $t$, and is
non-increasing when $t^{\prime}$ increases while $t$ is fixed because
$\frac{1}{\mathbb{P}_{ij}(Y^{t^{\prime}}_{ij}=1|\Delta_{ij}=1)}$ is non-
increasing in $t^{\prime}$. This implies the two-level monotonic property.

Next, we construct an alternative form of the Bayes decision function
$\bar{f}_{ij}^{t^{\prime}t}(\bm{x}_{ij})$, which is additive in stage $t$:

$\bar{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})=\bar{h}^{0}_{ij}(\bm{x}_{ij})-\sum_{r=t^{\prime}+1}^{t}\bar{h}^{r}_{ij}(\bm{x}_{ij}),$

where
$\bar{h}^{0}_{ij}(\bm{x}_{ij})=c_{ij}(\bm{x}_{ij})\log_{\alpha_{ij}(\bm{x}_{ij})}(2)>0$,
$\bar{h}^{r}_{ij}(\bm{x}_{ij})=c_{ij}(\bm{x}_{ij})\log_{\alpha_{ij}(\bm{x}_{ij})}\big{(}\mathbb{P}_{ij}(Y_{ij}^{t-1}=1|\Delta_{ij}=1)/\mathbb{P}_{ij}(Y_{ij}^{t}=1|\Delta_{ij}=1)\big{)}>0$,
and $c_{ij}(\bm{x}_{ij})>0$ and $\alpha_{ij}(\bm{x}_{ij})>0$ are two arbitrary
positive functions. Clearly, $\bar{\bar{f}}^{t^{\prime}t}_{ij}$ automatically
satisfies the two-level monotonic property. Now, we need to prove that
$\bar{\bar{f}}^{t^{\prime}t}_{ij}(\bm{x}_{ij})$ has the same sign of
$\bar{f}_{ij}^{t^{\prime}t}(\bm{x}_{ij})$, that is, $\mathop{\rm
sign}\big{(}\bar{\bar{f}}^{t^{\prime}t}_{ij}(\bm{x}_{ij})\big{)}=\mathop{\rm
sign}\big{(}c_{ij}(\bm{x}_{ij})\log_{\alpha_{ij}(\bm{x}_{ij})}\big{(}2\mathbb{P}_{ij}(Y^{t}_{ij}=1|Y_{ij}^{t^{\prime}}=1,\Delta_{ij}=1)\big{)}\big{)}=1$,
when $\mathbb{P}_{ij}(Y^{t}_{ij}=1|Y_{ij}^{t^{\prime}}=1,\Delta_{ij}=1)>1/2$
and equals to $-1$ otherwise. The desired result then follows. $\square$

Proof of Lemma 2. The result follows from the fact that
$\bar{\mathbf{f}}_{ij}^{t^{\prime}t}$ is a global minimizer of
$\mathbb{E}_{ij}\big{(}L(\mathbf{f}_{ij},\bm{Z}_{ij})\big{)}$, when $V(\cdot)$
is Fisher consistency for binary classification.

Proof of Lemma 3. Note that (3), (3), (3), (3), and (3) are convex
minimization problems. Then convergence of Algorithm 1 follows from (Tseng,
2001). $\square$

Proof of Theorem 4. Our treatment of bounding
$\mathbb{P}\big{(}e(\hat{\bm{\phi}})-e(\bar{\bm{\phi}})\geq\epsilon^{2}_{nm}\big{)}$
relies on a chain argument of empirical process over a suitable partition of
$\mathcal{F}$ induced by $e(\hat{\bm{\phi}})-e(\bar{\bm{\phi}})$, as in (Dai
et al., 2019a; Wong et al., 1995). For $u\geq 1$ and $v\geq 0$, define
$A_{uv}=\\{\mathbf{f}\in\mathcal{F}:2^{u-1}\epsilon^{2}_{nm}\leq
l_{B}(\mathbf{f})-l_{B}(\bar{\mathbf{f}})\leq
2^{u}\epsilon^{2}_{nm},2^{v-1}\bar{J}\leq J(f)\leq 2^{v}\bar{J}\\}$ and
$A_{u0}=\\{\mathbf{f}\in\mathcal{F}:2^{u-1}\epsilon^{2}_{nm}\leq
l_{B}(\mathbf{f})-l_{B}(\bar{\mathbf{f}})\leq
2^{u}\epsilon^{2}_{nm},J(f)\leq\bar{J}\\}$. By Assumption A,

$\displaystyle\mathbb{P}\big{(}e(\hat{\mathbf{f}})-e(\bar{\mathbf{f}})$
$\displaystyle\geq
c_{2}\epsilon_{nm}^{2\alpha}\big{)}\leq\mathbb{P}\big{(}l_{B}(\hat{\mathbf{f}})-l_{B}(\bar{\mathbf{f}})\geq\epsilon_{nm}^{2}\big{)}$
$\displaystyle\leq\mathbb{P}\big{(}\sup_{\\{\mathbf{f}\in\mathcal{F};l_{B}(\mathbf{f})-l_{B}(\bar{\mathbf{f}})\geq\epsilon^{2}_{nm}\\}}(nm)^{-1}\sum_{i,j}\big{(}\tilde{L}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-\tilde{L}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{)}\geq
0\big{)}$
$\displaystyle\leq\mathbb{P}\big{(}\sup_{\cup_{u=1}^{\infty}\cup_{v=0}^{\infty}A_{uv}}(nm)^{-1}\sum_{i,j}\big{(}\tilde{L}_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-\tilde{L}_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{)}\geq
0\big{)}$
$\displaystyle\leq\sum_{u=1}^{\infty}\sum_{v=1}^{\infty}\mathbb{P}\Big{(}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}(nm)^{-1}\sum_{i,j}g(\mathbf{f}_{ij},\bm{Z}_{ij})\geq\delta_{uv}\Big{)}$
$\displaystyle\quad+\sum_{u=1}^{\infty}\mathbb{P}\Big{(}\sup_{\\{\mathbf{f}\in
A_{u0}\\}}(nm)^{-1}\sum_{i,j}g(\mathbf{f}_{ij},\bm{Z}_{ij})\geq\delta_{u0}\Big{)}=\sum_{u=1}^{\infty}\sum_{v=1}^{\infty}I_{uv}+\sum_{u=1}^{\infty}I_{u0},$
(25)

where
$g(\mathbf{f}_{ij},\bm{Z}_{ij})=L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})-\big{(}L_{B}(\bar{\mathbf{f}}_{ij})-L_{B}({\mathbf{f}}_{ij})\big{)}$,
$\delta_{uv}=2^{u-1}\epsilon^{2}_{nm}+\lambda(2^{v-1}-1)\bar{J}$,
$\delta_{u0}=2^{u-2}\epsilon^{2}_{nm}$,
$I_{uv}=\mathbb{P}\big{(}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}(nm)^{-1}\sum_{i,j}g(\mathbf{f}_{ij},\bm{Z}_{ij})\geq\delta_{uv}\big{)}$,
and $\lambda
J(\mathbf{f})=\lambda_{1}\|\bm{\beta}\|_{2}^{2}+\sum_{k=1}^{K}\lambda_{k}\big{(}\|\bm{a}^{k}\|_{2}^{2}+\|\bm{b}^{k}\|_{2}^{2}\big{)}$.

To apply Talagrand’s inequality (Giné et al., 2006) to $I_{uv}$, let

$\sigma^{2}_{uv}=\sup_{\\{\mathbf{f}\in
A_{uv}\\}}(nm)^{-1}\sum_{i,j}\mathop{\rm
Var}_{ij}\big{(}g(\mathbf{f}_{ij},\bm{Z}_{ij})\big{)},\quad u\geq 1,v\geq 0.$

Then

$\displaystyle I_{uv}$
$\displaystyle\leq\mathbb{P}\Big{(}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}(nm)^{-1}\sum_{i,j}g(\mathbf{f}_{ij},\bm{Z}_{ij})\geq\delta_{uv}\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}(nm)^{-1}\big{|}\sum_{i,j}g(\mathbf{f}_{ij},\bm{Z}_{ij})\big{|}-\mathbb{E}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}(nm)^{-1}\big{|}\sum_{i,j}g(\mathbf{f}_{ij},\bm{Z}_{ij})\big{|}\leq\delta_{uv}-2\mathbb{E}\big{(}\mathcal{R}_{uv}\big{)}\Big{)}$
$\displaystyle\leq
a_{1}\exp(-\frac{nm\delta^{2}_{uv}}{2a_{1}(2\sigma_{uv}^{2}+16B\mathbb{E}(\mathcal{R}_{uv})+B\delta_{uv})})\leq
a_{1}\exp(-\frac{nm\delta^{2}_{uv}}{2a_{1}(4^{\mu}c_{2}\delta^{\mu}_{uv}+5B\delta_{uv})})$
$\displaystyle\leq
a_{1}\exp(-\frac{nm\delta^{2-\min(1,\mu)}_{uv}}{2a_{1}a_{2}}),$ (26)

where $a_{2}=32c_{2}B^{\mu-\min(1,\mu)}+9B^{2-\min(1,\mu)}$ and the second
inequality follows from the fact that

$\displaystyle\mathbb{E}\Big{(}(nm)^{-1}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}\big{|}\sum_{i,j}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})-\mathbb{E}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{)}\big{)}\big{|}\Big{)}$
$\displaystyle\leq\mathbb{E}\Big{(}(nm)^{-1}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}\big{|}\mathbb{E}\big{(}\sum_{i,j}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})-(L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}^{\prime}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}^{\prime}_{ij}))\big{|}\bm{Z}\big{)}\big{|}\Big{)}$
$\displaystyle\leq\mathbb{E}\Big{(}(nm)^{-1}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}\big{|}\sum_{i,j}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})-(L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}^{\prime}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}^{\prime}_{ij}))\big{|}\Big{)}$
$\displaystyle\leq\mathbb{E}\Big{(}(nm)^{-1}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}\big{|}\sum_{i,j}\tau_{i,j}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})-(L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}^{\prime}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}^{\prime}_{ij}))\big{|}\Big{)}$
$\displaystyle\leq 2\mathbb{E}\Big{(}(nm)^{-1}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}\big{|}\sum_{i,j}\tau_{i,j}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{|}\Big{)}\leq
2\mathbb{E}(\mathcal{R}_{uv}),$

where $\bm{Z}^{\prime}_{ij}$ is an independent copy of $\bm{Z}_{ij}$,
$\mathcal{R}_{uv}=\mathbb{E}_{\tau}\Big{(}(nm)^{-1}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}\big{|}\sum_{i,j}\tau_{i,j}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{|}\Big{)}$
is the Rademacher complexity of $A_{uv}$, and $\tau_{ij}$ are independent
random variables drawn from the Rademacher distribution. The third inequality
follows from the Talagrand’s inequality and the last inequality follows from
the fact that,

$\displaystyle\sigma^{2}_{uv}=\sup_{\\{\mathbf{f}\in
A_{uv}\\}}(nm)^{-1}\sum_{i,j}\text{Var}_{ij}\big{(}L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})-L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})\big{)}\leq
c_{2}(2^{u}\epsilon^{2}_{nm})^{\mu}\leq 16c_{2}\delta_{uv}^{\mu}.$ (27)

A combination of (Technical Proofs) and (Technical Proofs) yield

$\displaystyle\mathbb{P}\big{(}e(\hat{\mathbf{f}})-e(\bar{\mathbf{f}})\leq
c_{2}\epsilon^{2\alpha}_{n,m}\big{)}\leq\sum_{u=1}^{\infty}\sum_{v=1}^{\infty}I_{uv}+\sum_{u=1}^{\infty}I_{u0}$
$\displaystyle\leq\sum_{i=1}^{\infty}\sum_{v=1}^{\infty}a_{1}\exp\big{(}-\frac{nm\delta^{2-\min(1,\mu)}_{uv}}{2a_{1}a_{2}}\big{)}+\sum_{u=1}^{\infty}a_{1}\exp\big{(}-\frac{nm\delta^{2-\min(1,\mu)}_{u0}}{2a_{1}a_{2}}\big{)}$
$\displaystyle\leq
4a_{1}\exp\big{(}-\frac{nm(\lambda\bar{J})^{2-\min(1,\mu)}}{4a_{1}a_{2}}\big{)}.$

The desired result follows immediately. $\square$

###### Lemma 5

Let $A_{uv}=\\{\mathbf{f}\in\mathcal{F}:2^{u-1}\epsilon^{2}_{nm}\leq
l_{B}(\mathbf{f})-l_{B}(\bar{\mathbf{f}})\leq
2^{u}\epsilon^{2}_{nm},2^{v-1}\bar{J}\leq J(f)\leq 2^{v}\bar{J}\\}$,
$A_{u0}=\\{\mathbf{f}\in\mathcal{F}:2^{u-1}\epsilon^{2}_{nm}\leq
l_{B}(\mathbf{f})-l_{B}(\bar{\mathbf{f}})\leq
2^{u}\epsilon^{2}_{nm},J(f)\leq\bar{J}\\}$, and
$\delta_{uv}=2^{u-1}\epsilon^{2}_{nm}+\lambda(2^{v-1}-1)\bar{J}$, and
$\delta_{u0}=2^{u-2}\epsilon^{2}_{nm}$, where

$\epsilon^{2}_{nm}\sim\Big{(}\frac{(\Lambda+T)K}{nm}\log\big{(}\frac{nm\bar{w}B\bar{J}\kappa(6T^{3}+T^{2}K)}{(\Lambda+T)K}\big{)}\Big{)}^{\frac{1}{2-\min(1,\mu)}}.$

Then

$4\mathbb{E}(\mathcal{R}_{uv})\leq\delta_{uv},$ (28)

where
$\mathcal{R}_{uv}=\mathbb{E}_{\tau}\Big{(}(nm)^{-1}\sup_{\\{\mathbf{f}\in
A_{uv}\\}}\big{|}\sum_{i,j}\tau_{i,j}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{|}\Big{)}$
is the Rademacher complexity of $A_{uv}$.

Proof of Lemma 5. To bound $\mathbb{E}(\mathcal{R}_{uv})$, we compute the
metric entropy for

$\mathcal{L}_{uv}=\big{\\{}(nm)^{-1}\sum_{i,j}\big{(}L_{B}(\bar{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{)}:\mathbf{f}\in
A_{uv}\big{\\}}.$

For any functions $\tilde{\mathbf{f}}$ and $\mathbf{f}$,

$\displaystyle(nm)^{-1}\sum_{i,j}\big{(}L_{B}(\tilde{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{)}^{2}$
$\displaystyle=\Big{(}\sum_{t^{\prime}<t}w_{t^{\prime}t}\big{(}V_{B}(Y^{t}_{ij}\tilde{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij}))-V_{B}(Y^{t}_{ij}f^{t^{\prime}t}_{ij}(\bm{x}_{ij}))\big{)}\Big{)}^{2}$
$\displaystyle\leq\bar{w}\Big{(}\sum_{t^{\prime}<t}\big{(}\tilde{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})-f^{t^{\prime}t}_{ij}(\bm{x}_{ij})\big{)}^{2}\Big{)},$
(29)

where $\bar{w}=\sum_{0\leq t^{\prime}<t\leq T}w^{2}_{t^{\prime}t}$, and the
first inequality follows from the Cauchy-Schwarz inequality. Next, we turn to
bound
$\tilde{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})-f^{t^{\prime}t}_{ij}(\bm{x}_{ij})$,

$\displaystyle\big{(}\tilde{f}^{t^{\prime}t}_{ij}(\bm{x}_{ij})-f^{t^{\prime}t}_{ij}(\bm{x}_{ij})\big{)}^{2}=\Big{(}\big{(}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ\tilde{b}(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}\tilde{q}(t^{\prime},t)-\big{(}a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}q(t^{\prime},t)\Big{)}^{2}$
$\displaystyle\leq
2\Big{(}\big{(}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ\tilde{b}(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}\big{(}\sum_{r=t^{\prime}+1}^{t}(\bm{q}_{r}-\tilde{\bm{q}}_{r})\big{)}\Big{)}^{2}$
$\displaystyle\hskip
85.35826pt+2\Big{(}\big{(}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ\tilde{b}(\bm{v}_{j},\bm{o}_{j})-a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}q(t^{\prime},t)\Big{)}^{2}=I_{1}+I_{2},$

where $q(t^{\prime},t)=\bm{1}-\sum_{r=t^{\prime}+1}^{t}\bm{q}_{r}$. Then we
bound $I_{1}$ and $I_{2}$ separately. For $I_{1}$,

$\displaystyle I_{1}$ $\displaystyle\leq
2\big{\|}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ\tilde{b}(\bm{v}_{j},\bm{o}_{j})\big{\|}_{2}^{2}\big{\|}\sum_{r=t^{\prime}+1}^{t}\bm{q}_{r}-\tilde{\bm{q}}_{r}\big{\|}_{2}^{2}\leq
2\big{\|}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\big{\|}_{2}^{2}\big{\|}\tilde{b}(\bm{v}_{j},\bm{o}_{j})\big{\|}_{2}^{2}\big{\|}\sum_{r=t^{\prime}+1}^{t}\bm{q}_{r}-\tilde{\bm{q}}_{r}\big{\|}_{2}^{2}$
$\displaystyle\leq
8(t-t^{\prime})(p_{1}+d_{1})(p_{2}+d_{2})\big{(}\|\bm{A}\|_{F}^{2}+\|\bm{a}\|_{2}^{2}\big{)}\big{(}\|\bm{B}\|_{F}^{2}+\|\bm{b}\|_{2}^{2}\big{)}\big{(}\sum_{r=t^{\prime}+1}^{t}\big{\|}\bm{q}_{r}-\tilde{\bm{q}}_{r}\big{\|}_{2}^{2}\big{)}$
$\displaystyle\leq
2^{2v+3}\bar{J}^{2}(t-t^{\prime})(p_{1}+d_{1})(p_{2}+d_{2})\big{(}\sum_{r=t^{\prime}+1}^{t}\big{\|}\bm{q}_{r}-\tilde{\bm{q}}_{r}\big{\|}_{2}^{2}\big{)},$
(30)

where the last inequality follows from the fact that

$\big{\|}a(\bm{u}_{i},\bm{s}_{i})\big{\|}_{2}^{2}\leq
2\big{(}\|\bm{A}\|_{F}^{2}\|\bm{u}_{i}\|_{2}^{2}+d_{1}\sum_{l=1}^{d_{1}}\max_{i=1,\cdots,n_{l}}\|\bm{a}_{li}\|_{2}^{2}\big{)}\leq
2\big{(}p_{1}\|\bm{A}\|_{F}^{2}+d_{1}\|\bm{a}\|_{2}^{2}\big{)},$

and

$\big{\|}b(\bm{v}_{j},\bm{o}_{j})\big{\|}_{2}^{2}\leq
2\big{(}\|\bm{B}\|_{F}^{2}\|\bm{v}_{j}\|_{2}^{2}+d_{2}\sum_{l=1}^{d_{2}}\max_{j=1,\cdots,m_{l}}\|\bm{b}_{lj}\|_{2}^{2}\big{)}\leq
2\big{(}p_{2}\|\bm{B}\|_{F}^{2}+d_{2}\|\bm{b}\|_{2}^{2}\big{)}.$

For $I_{2}$,

$\displaystyle I_{2}$ $\displaystyle\leq
2\Big{(}\big{(}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ\tilde{b}(\bm{v}_{j},\bm{o}_{j})-a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{)}^{\top}q(t^{\prime},t)\Big{)}^{2}$
$\displaystyle\leq
2\big{\|}\bm{1}-\sum_{r=t^{\prime}+1}^{t}\bm{q}_{r}\big{\|}_{2}^{2}\big{\|}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ\tilde{b}(\bm{v}_{j},\bm{o}_{j})-a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\big{\|}_{2}^{2}$ $\displaystyle\leq
4\big{(}K+(t-t^{\prime})\sum_{r=t^{\prime}+1}^{t}\|\bm{q}_{r}\|_{2}^{2}\big{)}\Big{\|}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ\tilde{b}(\bm{v}_{j},\bm{o}_{j})-\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\Big{\|}_{2}^{2}$
$\displaystyle\quad+4\big{(}K+(t-t^{\prime})\sum_{r=t^{\prime}+1}^{t}\|\bm{q}_{r}\|_{2}^{2}\big{)}\Big{\|}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})-a(\bm{u}_{i},\bm{s}_{i})\circ
b(\bm{v}_{j},\bm{o}_{j})\Big{\|}_{2}^{2}$ $\displaystyle\leq
4\big{(}K+(t-t^{\prime})\sum_{r=t^{\prime}+1}^{t}\|\bm{q}_{r}\|_{2}^{2}\big{)}\big{\|}\tilde{a}(\bm{u}_{i},\bm{s}_{i})\big{\|}_{2}^{2}\big{\|}(\tilde{\bm{B}}-\bm{B})\bm{v}_{j}-\sum_{l=1}^{d_{2}}(\tilde{\bm{b}}_{lo_{jl}}-\bm{b}_{lo_{jl}})\big{\|}_{2}^{2}$
$\displaystyle\quad+4\big{(}K+(t-t^{\prime})\sum_{r=t^{\prime}+1}^{t}\|\bm{q}_{r}\|_{2}^{2}\big{)}\big{\|}b(\bm{v}_{j},\bm{o}_{j})\big{\|}_{2}^{2}\big{\|}(\tilde{\bm{A}}-\bm{A})\bm{u}_{i}-\sum_{l=1}^{d_{1}}(\tilde{\bm{a}}_{ls_{il}}-\bm{a}_{ls_{il}})\big{\|}_{2}^{2}$
$\displaystyle\leq
2^{v+4}\kappa\bar{J}\big{(}K+(t-t^{\prime})\sum_{r=t^{\prime}+1}^{t}\|\bm{q}_{r}\|_{2}^{2}\big{)}$
$\displaystyle\hskip
113.81102pt\Big{(}p_{1}\|\tilde{\bm{A}}-\bm{A}\|_{F}^{2}+p_{2}\|\tilde{\bm{B}}-\bm{B}\|_{F}^{2}+d_{1}\|\tilde{\bm{a}}-\bm{a}\|_{2}^{2}+d_{2}\|\tilde{\bm{b}}-\bm{b}\|_{2}^{2}\Big{)},$
(31)

where $\kappa=p_{1}+p_{2}+d_{1}+d_{2}$. Combining (Technical Proofs), (30) and
(31),

$\displaystyle(nm)^{-1}\sum_{i,j}\big{(}L_{B}(\tilde{\mathbf{f}}_{ij},\bm{Z}_{ij})-L_{B}(\mathbf{f}_{ij},\bm{Z}_{ij})\big{)}^{2}$
$\displaystyle\leq\bar{w}\Big{(}2^{2v+6}\bar{J}^{2}(p_{1}+d_{1})(p_{2}+d_{2})\sum_{t^{\prime}<t}(t-t^{\prime})\sum_{r=t^{\prime}+1}^{t}\big{\|}\bm{q}_{r}-\tilde{\bm{q}}_{r}\big{\|}_{2}^{2}$
$\displaystyle\
+\kappa^{2}2^{v+4}\bar{J}\big{(}\|\tilde{\bm{A}}-\bm{A}\|_{F}^{2}+\|\tilde{\bm{B}}-\bm{B}\|_{F}^{2}+\|\tilde{\bm{a}}-\bm{a}\|_{2}^{2}+\|\tilde{\bm{b}}-\bm{b}\|_{2}^{2}\big{)}\big{(}\sum_{t^{\prime}<t}\big{(}K+(t-t^{\prime})\sum_{r=t^{\prime}+1}^{t}\|\bm{q}_{r}\|_{2}^{2}\big{)}\big{)}\Big{)}$
$\displaystyle\leq\bar{w}\Big{(}2^{2v+2}\bar{J}^{2}T^{3}(p_{1}+d_{1})(p_{2}+d_{2})\big{\|}\bm{q}-\tilde{\bm{q}}\big{\|}_{2}^{2}$
$\displaystyle\
+\kappa^{2}2^{v}\bar{J}\big{(}(T+1)TK+T^{3}\|\bm{q}\|_{2}^{2}\big{)}\big{(}\|\tilde{\bm{A}}-\bm{A}\|_{F}^{2}+\|\tilde{\bm{B}}-\bm{B}\|_{F}^{2}+\|\tilde{\bm{a}}-\bm{a}\|_{2}^{2}+\|\tilde{\bm{b}}-\bm{b}\|_{2}^{2}\big{)}\Big{)}.$

Therefore, the metric entropy for $\mathcal{L}_{uv}$ for $\nu\leq 1$ is upper
bounded:

$\displaystyle\log\mathcal{N}(\nu,\mathcal{L}_{uv})$
$\displaystyle\leq\log\mathcal{N}\Big{(}\frac{\nu}{2^{v}\bar{J}Ba_{3}},\mathcal{B}_{+}\big{(}(\Lambda+T)K\big{)}\Big{)}\leq\gamma\log\big{(}\frac{2^{v+1}\bar{J}Ba_{3}}{\nu}\big{)},$
(32)

where $a_{3}=\kappa\sqrt{\bar{w}(6T^{3}+(T+1)TK)}$,
$\mathcal{B}_{+}((\Lambda+T)K)$ is the unit $l_{2}$-ball in
$\mathbb{R}_{+}^{(\Lambda+T)K}$. Therefore, the Rademacher complexity
$\mathbb{E}(\mathcal{R}_{uv})$ can be computed based on Theorem 3.12 of
(Koltchinskii, 2011). Specifically, there exists a constant $a_{4}$, such that

$\mathbb{E}\big{(}\mathcal{R}_{uv}\big{)}\leq
a_{4}\max\Big{(}\frac{\delta^{\mu/2}_{uv}\sqrt{(\Lambda+T)K}}{\sqrt{nm}}\sqrt{\log\big{(}\frac{2^{v+1}\bar{J}Ba_{3}}{\delta^{\mu/2}_{uv}}\big{)}},\frac{B(\Lambda+T)K}{nm}\log\big{(}\frac{2^{v+1}\bar{J}Ba_{3}}{\delta^{\mu/2}_{uv}}\big{)}\Big{)}.$

Then there exists a constant $a_{5}>0$, such that
$\epsilon^{2}_{nm}=a_{5}\Big{(}\frac{(\Lambda+T)K}{nm}\log\big{(}\frac{nm\bar{w}B\bar{J}\kappa(6T^{3}+T^{2}K)}{(\Lambda+T)K}\big{)}\Big{)}^{\frac{1}{2-\min(1,\mu)}}$,
and

$\displaystyle\frac{\delta_{uv}}{4\mathbb{E}\big{(}\mathcal{R}_{uv}\big{)}}$
$\displaystyle\geq\frac{1}{4a_{4}}\min\Big{(}\frac{\sqrt{nm}\delta^{1-\mu/2}_{uv}}{\sqrt{(\Lambda+T)K}}\big{(}\log\big{(}\frac{2^{v+1}\bar{J}Ba_{3}}{\delta^{\mu/2}_{uv}}\big{)}\big{)}^{-1/2},\frac{nm\delta_{uv}}{B(\Lambda+T)K}\big{(}\log\big{(}\frac{2^{v+1}\bar{J}Ba_{3}}{\delta^{\mu/2}_{uv}}\big{)}\big{)}^{-1}\Big{)}$
$\displaystyle\geq\frac{1}{4a_{4}}\min\Big{(}\frac{\sqrt{nm}\delta^{1-\mu/2}_{u0}}{\sqrt{(\Lambda+T)K}}\big{(}\log\big{(}\frac{2\bar{J}Ba_{3}}{\delta^{\mu/2}_{u0}}\big{)}\big{)}^{-1/2},\frac{nm\delta_{u0}}{B(\Lambda+T)K}\big{(}\log\big{(}\frac{2\bar{J}Ba_{3}}{\delta^{\mu/2}_{u0}}\big{)}\big{)}^{-1}\Big{)}$
$\displaystyle\geq\frac{1}{4a_{4}}\min\Big{(}\frac{\sqrt{nm}\epsilon^{2-\mu}_{nm}}{\sqrt{(\Lambda+T)K}}\big{(}\log\big{(}\frac{2\bar{J}Ba_{3}}{\epsilon^{\mu}_{nm}}\big{)}\big{)}^{-1/2},\frac{nm\epsilon^{2}_{nm}}{B(\Lambda+T)K}\big{(}\log\big{(}\frac{2\bar{J}Ba_{3}}{\epsilon^{\mu}_{nm}}\big{)}\big{)}^{-1}\Big{)}\geq
1.$

where the first two inequalities follow from non-decreasing of $(u,v)$. This
completes the proof. $\square$

## References

  * Bhaskar (2016) S. A. Bhaskar.  Probabilistic low-rank matrix completion from quantized measurements.  _The Journal of Machine Learning Research_ , 17(1):2131–2164, 2016.
  * Bi et al. (2018) X. Bi, A. Qu, and X. Shen.  Multilayer tensor factorization with applications to recommender systems.  _The Annals of Statistics_ , 46(6B):3308 – 3333, 2018.  doi: 10.1214/17-AOS1659.  URL https://doi.org/10.1214/17-AOS1659.
  * Cawley and Talbot (2010) G. C. Cawley and N. L. Talbot.  On over-fitting in model selection and subsequent selection bias in performance evaluation.  _Journal of Machine Learning Research_ , 11(Jul):2079–2107, 2010.
  * Chen et al. (2012) B. Chen, S. He, Z. Li, and S. Zhang.  Maximum block improvement and polynomial optimization.  _SIAM Journal on Optimization_ , 22(1):87–107, 2012.
  * Cheng et al. (2014) C. Cheng, F. Xia, T. Zhang, I. King, and M. R. Lyu.  Gradient boosting factorization machines.  In _Proceedings of the 8th ACM Conference on Recommender systems_ , pages 265–272, 2014.
  * Cortes and Vapnik (1995) C. Cortes and V. Vapnik.  Support vector networks.  _Machine Learning_ , 20(3):273–297, 1995.
  * Crammer and Singer (2002) K. Crammer and Y. Singer.  Pranking with ranking.  In _Advances in Neural Information Processing Systems_ , pages 641–647, 2002.
  * Dai et al. (2019a) B. Dai, X. Shen, J. Wang, and A. Qu.  Scalable collaborative ranking for personalized prediction.  _Journal of the American Statistical Association_ , Jan. 2019a.  ISSN 0162-1459.  doi: 10.1080/01621459.2019.1691562.
  * Dai et al. (2019b) B. Dai, J. Wang, X. Shen, and A. Qu.  Smooth neighborhood recommender systems.  _Journal of Machine Learning Research_ , 20(16):1–24, 2019b.  URL http://jmlr.org/papers/v20/17-629.html.
  * Edwards (2012) D. Edwards.  _Introduction to graphical modelling_.  Springer Science & Business Media, 2012.
  * Friedman (2001) J. H. Friedman.  Greedy function approximation: a gradient boosting machine.  _Annals of statistics_ , pages 1189–1232, 2001.
  * Giné et al. (2006) E. Giné, V. Koltchinskii, et al.  Concentration inequalities and asymptotic results for ratio type empirical processes.  _Annals of Probability_ , 34(3):1143–1216, 2006\.
  * Haeffele and Vidal (2015) B. D. Haeffele and R. Vidal.  Global optimality in tensor factorization, deep learning, and beyond.  _arXiv preprint arXiv:1506.07540_ , 2015.
  * He et al. (2017) X. He, L. Liao, H. Zhang, L. Nie, X. Hu, and T.-S. Chua.  Neural collaborative filtering.  In _Proceedings of the 26th International Conference on World Wide Web_ , pages 173–182, 2017.
  * Johnson (2014) C. C. Johnson.  Logistic matrix factorization for implicit feedback data.  _Advances in Neural Information Processing Systems_ , 27, 2014.
  * Koltchinskii (2011) V. Koltchinskii.  _Oracle Inequalities in Empirical Risk Minimization and Sparse Recovery Problems: Ecole d’Eté de Probabilités de Saint-Flour XXXVIII-2008_ , volume 2033.  Springer Science & Business Media, 2011.
  * Koren (2008) Y. Koren.  Factorization meets the neighborhood: a multifaceted collaborative filtering model.  In _Proceedings of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 426–434, 2008.
  * Le and Mikolov (2014) Q. Le and T. Mikolov.  Distributed representations of sentences and documents.  In _International Conference on Machine Learning_ , pages 1188–1196, 2014.
  * McCullagh and Nelder (2019) P. McCullagh and J. A. Nelder.  _Generalized linear models_.  Routledge, 2019.
  * Schmidhuber (2015) J. Schmidhuber.  Deep learning in neural networks: An overview.  _Neural Networks_ , 61:85–117, 2015.
  * Shen et al. (2003) X. Shen, G. C. Tseng, X. Zhang, and W. H. Wong.  On $\psi$-learning.  _Journal of the American Statistical Association_ , 98(463):724–734, 2003.
  * Suphavilai et al. (2018) C. Suphavilai, D. Bertrand, and N. Nagarajan.  Predicting cancer drug response using a recommender system.  _Bioinformatics_ , 34(22):3907–3914, 2018.
  * Tran et al. (2012) T. Tran, D. Q. Phung, and S. Venkatesh.  A sequential decision approach to ordinal preferences in recommender systems.  In _Twenty-Sixth AAAI Conference on Artificial Intelligence_ , 2012\.
  * Tseng (2001) P. Tseng.  Convergence of a block coordinate descent method for nondifferentiable minimization.  _Journal of optimization theory and applications_ , 109(3):475–494, 2001.
  * Tsybakov et al. (2004) A. B. Tsybakov et al.  Optimal aggregation of classifiers in statistical learning.  _Annals of Statistics_ , 32(1):135–166, 2004\.
  * Wan and McAuley (2018) M. Wan and J. McAuley.  Item recommendation on monotonic behavior chains.  In _Proceedings of the 12th ACM Conference on Recommender Systems_ , pages 86–94. ACM, 2018.
  * Weinberger et al. (2009) K. Weinberger, A. Dasgupta, J. Attenberg, J. Langford, and A. Smola.  Feature hashing for large scale multitask learning.  _arXiv preprint arXiv:0902.2206_ , 2009.
  * Wong et al. (1995) W. H. Wong, Shen, and Xiaotong.  Probability inequalities for likelihood ratios and convergence rates of sieve mles.  _Annals of Statistics_ , 23(2):339–362, 1995\.
  * Wright (2015) S. J. Wright.  Coordinate descent algorithms.  _Mathematical Programming_ , 151(1):3–34, 2015\.
  * Zhang and Liu (2014) C. Zhang and Y. Liu.  Multicategory angle-based large-margin classification.  _Biometrika_ , 101(3):625–640, 2014.
  * Zhu and Hastie (2002) J. Zhu and T. Hastie.  Kernel logistic regression and the import vector machine.  In _Advances in Neural Information Processing Systems_ , pages 1081–1088, 2002.