# Parallel Algorithm for Optimal Threshold Labeling of Ordinal Regression
Methods

Ryoya Yamasaki &Toshiyuki Tanaka

###### Abstract

Ordinal regression (OR) is classification of ordinal data in which the
underlying categorical target variable has a natural ordinal relation for the
underlying explanatory variable. For $K$-class OR tasks, threshold methods
learn a one-dimensional transformation (1DT) of the explanatory variable so
that 1DT values for observations of the explanatory variable preserve the
order of label values $1,\ldots,K$ for corresponding observations of the
target variable well, and then assign a label prediction to the learned 1DT
through threshold labeling, namely, according to the rank of an interval to
which the 1DT belongs among intervals on the real line separated by $(K-1)$
threshold parameters. In this study, we propose a parallelizable algorithm to
find the optimal threshold labeling, which was developed in previous research,
and derive sufficient conditions for that algorithm to successfully output the
optimal threshold labeling. In a numerical experiment we performed, the
computation time taken for the whole learning process of a threshold method
with the optimal threshold labeling could be reduced to approximately 60 % by
using the proposed algorithm with parallel processing compared to using an
existing algorithm based on dynamic programming.

##  1 Introduction

_Ordinal regression_ (_OR_ , or called ordinal classification) is
classification of _ordinal data_ in which the underlying target variable is
categorical and considered to be equipped with a _natural ordinal relation_
for the underlying explanatory variable. Various applications, including age
estimation (Niu et al., 2016; Cao et al., 2020; Yamasaki, 2023), information
retrieval (Liu, 2011), movie rating (Yu et al., 2006), and questionnaire
survey (Bürkner and Vuorre, 2019), leverage OR techniques.

_Threshold methods_ (or called threshold models) have been actively studied in
machine learning research (Shashua and Levin, 2003; Lin and Li, 2006; Chu and
Keerthi, 2007; Lin and Li, 2012; Li and Lin, 2007; Pedregosa et al., 2017;
Yamasaki, 2023), and they are popularly used for OR tasks as a simple way to
capture the ordinal relation of ordinal data. Those methods learn a _one-
dimensional transformation_ (_1DT_) of the observation of the explanatory
variable so that an observation with a larger class label tends to have a
larger 1DT value; they predict a label value through _threshold labeling_ ,
that is, according to the rank of an interval to which the learned 1DT belongs
among intervals on the real line separated by $(K-1)$ _threshold parameters_
for a $K$-class OR task.

Yamasaki (2023) proposed to use a threshold parameter vector that minimizes
the empirical task risk, and showed experimentally that the corresponding
threshold labeling function (_optimal threshold labeling_) could lead to
better classification performance than other labeling functions applied in
previous OR methods (McCullagh, 1980; Shashua and Levin, 2002; Chu and
Keerthi, 2005; Lin and Li, 2006, 2012; Pedregosa et al., 2017; Cao et al.,
2020). The previous study (Yamasaki, 2023) employs a _dynamic-programming-
based (DP) algorithm_ (Lin and Li, 2006) for optimization of the threshold
parameter vector, but that algorithm does not suit acceleration based on
parallel processing and requires long computation time when the training data
size and number of classes are large. In order to mitigate this trouble, in
this study, we propose another parallelizable algorithm, which we call
_independent optimization_ (_IO_) _algorithm_ , to find the optimal threshold
labeling, and derive sufficient conditions for that algorithm to successfully
output the optimal threshold labeling. Moreover, we demonstrate, through a
numerical experiment, that the paralleled IO algorithm can reduce the
computation time.

The rest of this paper is organized as follows. Section 2 introduces
formulations of the OR task, threshold method, and optimal threshold labeling,
in preparation for discussing calculation algorithms for the optimal threshold
labeling. Section 3 provides a review of an existing calculation algorithm,
the DP algorithm. Section 4 presents the main results of this paper, IO
algorithm and its parallelization procedure and theoretical properties.
Section 5 gives an experimental comparison of the computational efficiency of
the DP, non-parallelized IO, and parallelized IO algorithms. Section 6
concludes this paper. Additionally, in Appendix A, we mention another
parallelization procedure of the IO algorithm, which is more complex than that
described in the main text.

##  2 Preparation

###  2.1 Formulation of Ordinal Regression Task

Suppose that we are interested in behaviors of the categorical target variable
$Y\in[K]:=\\{1,\ldots,K\\}$ given a value of the explanatory variable
${\bm{X}}\in{\mathbb{R}}^{d}$, and that we have ordinal data
$({\bm{x}}_{1},y_{1}),\ldots,({\bm{x}}_{n},y_{n})\in{\mathbb{R}}^{d}\times[K]$
that are considered to be independent and identically distributed observations
of the pair $({\bm{X}},Y)$ and have a natural ordinal relation like examples
described in Section 1, with $K,d,n\in{\mathbb{N}}$ such that $K\geq 3$.111In
this paper, we do not discuss what the natural ordinal relation is. Refer to,
for example, da Costa et al. (2008); Yamasaki (2022) for this topic. An OR
task is a classification task of the ordinal data. In this paper, we focus on
a popular formulation of the OR task, searching for a classifier
$f:{\mathbb{R}}^{d}\to[K]$ that is good in minimization of the _task risk_
$\operatorname{{\mathbb{E}}}[\ell(f({\bm{X}}),Y)]$ for a user-specified _task
loss function_ $\ell:[K]^{2}\to[0,+\infty)$, where the expectation
$\operatorname{{\mathbb{E}}}[\cdot]$ is taken for $({\bm{X}},Y)$. Popular task
loss functions include not only the zero-one task loss $\ell_{\rm
zo}(k,l):=\operatorname{{\mathbbm{1}}}(k\neq l)$ (for minimization of
misclassification rate), but also V-shaped losses (for cost-sensitive
classification tasks) reflecting the user’s preference of smaller prediction
errors over larger ones such as the absolute task loss $\ell_{\rm
ab}(k,l):=|k-l|$, and squared task loss $\ell_{\rm sq}(k,l):=(k-l)^{2}$, where
$\operatorname{{\mathbbm{1}}}(c)$ is 1 if the condition $c$ is true and 0
otherwise.

In a formulation of the OR task, one may adopt other criteria that cannot be
decomposed into a sum of losses for each data point: for example, quadratic
weighted kappa (Cohen, 1960, 1968). Our discussion in this paper does not
cover such criteria.

###  2.2 Formulation of Threshold Method

_1DT-based methods_ (Yamasaki, 2023) learn a 1DT
$a:{\mathbb{R}}^{d}\to{\mathbb{R}}$ of the explanatory variable ${\bm{X}}$ so
that learned 1DT values $\hat{a}({\bm{x}}_{i})$’s preserve the order of label
values $y_{i}$’s well (McCullagh, 1980; Shashua and Levin, 2002; Chu and
Keerthi, 2005; Lin and Li, 2006, 2012; Pedregosa et al., 2017; Cao et al.,
2020; Yamasaki, 2022, 2023). For example, _ordinal logistic regression (OLR)_
(McCullagh, 1980) models the probability of $Y=y$ conditioned on
${\bm{X}}={\bm{x}}$ by $P(y,a({\bm{x}}),{\bm{b}})=\sigma(b_{1}-a({\bm{x}}))$
(for $y=1$), $\sigma(b_{y}-a({\bm{x}}))-\sigma(b_{y-1}-a({\bm{x}}))$ (for
$y\in\\{2,\ldots,K-1\\}$), $1-\sigma(b_{K-1}-a({\bm{x}}))$ (for $y=K$) with a
1DT $a$, a bias parameter vector
${\bm{b}}=(b_{k})_{k\in[K-1]}\in{\mathbb{R}}^{K-1}$ satisfying
$b_{1}\leq\cdots\leq b_{K-1}$, and the sigmoid function
$\sigma(u):=1/(1+e^{-u})$, and it can learn $a$ and ${\bm{b}}$ via, for
example, the maximum likelihood method
$(\hat{a},\hat{{\bm{b}}})\in\operatorname{{arg\,max}}_{a,{\bm{b}}}\prod_{i=1}^{n}P(y_{i},a({\bm{x}}_{i}),{\bm{b}})$.
1DT-based methods construct a classifier $f$ as $f=h\circ\hat{a}$ with a
learned 1DT $\hat{a}$ and a _labeling function_ $h:{\mathbb{R}}\to[K]$. Many
of existing 1DT-based methods can be seen as adopting a _threshold labeling
function_ ,

$h_{\rm thr}(u;{\bm{t}}):=1+\sum_{k=1}^{K-1}\operatorname{{\mathbbm{1}}}(u\geq
t_{k})$ (1)

with a threshold parameter vector
${\bm{t}}=(t_{k})_{k\in[K-1]}\in{\mathbb{R}}^{K-1}$; we call such
classification methods threshold methods.

###  2.3 Optimal Threshold Labeling

Yamasaki (2023) proposed to use a threshold parameter vector $\hat{{\bm{t}}}$
that minimizes the empirical task risk for a learned 1DT $\hat{a}$ and
specified task loss $\ell$:

$\hat{{\bm{t}}}\in\underset{{\bm{t}}\in{\mathbb{R}}^{K-1}}{\operatorname{{arg\,min}}}\,\frac{1}{n}\sum_{i=1}^{n}\ell(h_{\rm
thr}(\hat{a}({\bm{x}}_{i});{\bm{t}}),y_{i}).$ (2)

It was also shown experimentally that, for various learning methods of the
1DT, the optimal threshold labeling function $h_{\rm thr}(u;\hat{{\bm{t}}})$
could yield smaller test task risk than other labeling functions used in
previous studies (McCullagh, 1980; Shashua and Levin, 2002; Chu and Keerthi,
2005; Lin and Li, 2006, 2012; Pedregosa et al., 2017; Cao et al., 2020), for
example, $h(u)=h_{\rm thr}(u;\hat{{\bm{b}}})$ with a learned bias parameter
vector $\hat{{\bm{b}}}$ and likelihood-based labeling function
$h(u)\in\operatorname{{arg\,min}}_{k\in[K]}(\sum_{l=1}^{K}\ell(k,l)P(l,u,\hat{{\bm{b}}}))$
for OLR reviewed in Section 2.2.

It is important to reduce the computation time required to optimize the
threshold parameter vector. For example, when a user learns the 1DT with an
iterative optimization algorithm and employs early stopping (Prechelt, 2002),
the user evaluates the empirical task risk using a holdout validation dataset
at every epoch. For this purpose it is necessary to calculate the optimal
threshold parameter vector for a 1DT available at each epoch, and the
computation time for this calculation can account for a high percentage of
that for the whole learning process for a threshold method, which will be
shown by a numerical experiment we took (see Table 1). Therefore, we study
computationally efficient algorithms for calculating the optimal threshold
parameter vector.

##  3 Existing Dynamic-Programming-based (DP) Algorithm

It was shown in Yamasaki (2023, Theorem 6) that the minimization problem (2)
can be solved by the DP algorithm (Algorithms 1 and 2), which was developed by
Lin and Li (2006), as

$\begin{split}&\hat{{\bm{t}}}\in\underset{{\bm{t}}\in{\mathbb{R}}^{K-1}}{\operatorname{{arg\,min}}}\,\frac{1}{n}\sum_{i=1}^{n}\ell(h_{\rm
thr}(\hat{a}({\bm{x}}_{i});{\bm{t}}),y_{i}),\\\
&\text{s.t.~{}}t_{1},\ldots,t_{K-1}\in\\{c_{j}\\}_{j\in[N+1]}\text{~{}and~{}}t_{1}\leq\cdots\leq
t_{K-1}\end{split}$ (3)

with the _candidate vector_ ${\bm{c}}=(c_{j})_{j\in[N+1]}$. Here, $N$ is the
total number of the sorted unique elements $\hat{a}_{j}$, $j=1,\ldots,N$ of
$\\{\hat{a}({\bm{x}}_{i})\\}_{i\in[n]}$ that satisfy
$\hat{a}_{1}<\cdots<\hat{a}_{N}$ and
$\hat{a}({\bm{x}}_{i})\in\\{\hat{a}_{j}\\}_{j\in[N]}$ for all $i\in[n]$, and
we define ${\bm{c}}$ by $c_{j}:=-\infty$ (for $j=1$),
$(\hat{a}_{j-1}+\hat{a}_{j})/2$ (for $j\in\\{2,\ldots,N\\}$), $+\infty$ (for
$j=N+1$). Note that there is a degree of freedom in definition of the
candidate vector ${\bm{c}}$, and that the value of the training task risk does
not change even if $c_{1}$, $c_{j}$ for $j\in\\{2,\ldots,N\\}$, and $c_{N+1}$
are replaced with a smaller value than $\min(\hat{a}_{j})_{j\in[N]}$, a value
in $(\hat{a}_{j-1},\hat{a}_{j})$, and a larger value than
$\max(\hat{a}_{j})_{j\in[N]}$.

Algorithm 2 does not suit parallel processing, and takes the computation time
of the order $O(N\cdot K)$. Note that Algorithm 1 takes the computation time
of the order $O(n\log n)$ in average by, for example, quick sort, but the
actual computation time for Algorithm 2 can be longer than that for Algorithm
1; see experimental results in Section 5.

Input: Training data $\\{({\bm{x}}_{i},y_{i})\\}_{i\in[n]}$, learned 1DT
$\hat{a}$, and task loss $\ell$.

1 Calculate $\hat{a}({\bm{x}}_{i})$, $i=1,\ldots,n$;

2 Calculate sorted unique elements $\hat{a}_{j}$, $j=1,\ldots,N$ of
$\\{\hat{a}({\bm{x}}_{i})\\}_{i\in[n]}$;

3 Calculate sets
${\mathcal{Y}}_{j}=\\{y_{i}\mid\hat{a}({\bm{x}}_{i})=\hat{a}_{j}\\}_{i\in[n]}$,
$j=1,\ldots,N$;

/* Calculate loss matrix $(M_{j,k})_{j\in[N],k\in[K]}$. Parallelizable with
$j\in[N]$ and $k\in[K]$. */

4 for  _$j=1,\ldots,N$ and $k=1,\ldots,K$_ do
$M_{j,k}=\sum_{y_{i}\in{\mathcal{Y}}_{j}}\ell(k,y_{i})$;

5 Define the candidate vector ${\bm{c}}$ as $c_{j}:=-\infty$ (for $j=1$),
$(\hat{a}_{j-1}+\hat{a}_{j})/2$ (for $j\in\\{2,\ldots,N\\}$), $+\infty$ (for
$j=N+1$);

Output: $\hat{a}_{j}$ and ${\mathcal{Y}}_{j}$ ($j=1,\ldots,N$), loss matrix
$M$, and candidate vector ${\bm{c}}$.

Algorithm 1 Preparation for Algorithms 2, 3, and 5

Input: $\hat{a}_{j}$ and ${\mathcal{Y}}_{j}$ ($j=1,\ldots,N$), loss matrix
$M$, and candidate vector ${\bm{c}}$ prepared by Algorithm 1.

/* Calculate $(L_{j,k})_{j\in[N],k\in[K]}$ sequentially. */

1 for  _$k=1,\ldots,K$_ do $L_{1,k}=M_{1,k}$;

2 for  _$j=2,\ldots,N$_ do

/* An efficient implementation of $L_{j,k}=\min(L_{j-1,l})_{l\in[k]}+M_{j,k}$.
*/

3 $O_{j}\leftarrow+\infty$;

4 for  _$k=1,\ldots,K$_ do

5 if  _$L_{j-1,k} <O_{j}$_ then $O_{j}\leftarrow L_{j-1,k}$;

6 $L_{j,k}=O_{j}+M_{j,k}$;

/* Calculate threshold parameters $(\hat{t}_{k})_{k\in[K-1]}$ sequentially. */

7 $I\leftarrow\min(\operatorname{{arg\,min}}(L_{N,k})_{k\in[K]})$;

8 if  _$I\neq K$_ then $\hat{t}_{k}=c_{N+1}$ for $k=I,\ldots,K-1$;

9 for  _$j=N-1,\ldots,1$_ do

10 $J\leftarrow\min(\operatorname{{arg\,min}}(L_{j,k})_{k\in[I]})$;

11 if  _$I\neq J$_ then $\hat{t}_{k}=c_{j+1}$ for $k=J,\ldots,I-1$, and
$I\leftarrow J$;

12if  _$I\neq 1$_ then $\hat{t}_{k}=c_{1}$ for $k=1,\ldots,I-1$;

Output: A threshold parameter vector $(\hat{t}_{k})_{k\in[K-1]}$.

Algorithm 2 DP algorithm to calculate the optimal threshold parameter
vector333Optimality is guaranteed regardless of which element of
$\operatorname{{arg\,min}}(L_{N,k})_{k\in[K]}$ is set to $I$ in Line 2 and
which element of $\operatorname{{arg\,min}}(L_{j,k})_{k\in[I]}$ is set to $J$
in Line 2. This is also true for Line 3 of Algorithm 3 and Line 4 of Algorithm
5 described in Appendix A.

##  4 Proposed Independent Optimization (IO) Algorithm

We here propose to learn the threshold parameters $t_{k}$, $k\in[K-1]$
independently by Algorithm 3 (IO algorithm), which is parallelizable with
$k\in[K-1]$. This algorithm is developed according to the conditional relation

$\frac{1}{n}\sum_{i=1}^{n}\ell(h_{\rm
thr}(\hat{a}({\bm{x}}_{i});{\bm{t}}),y_{i})\\\
=\sum_{k=1}^{K-1}R_{k}(t_{k})-\underbrace{\sum_{k=2}^{K-1}R_{k}(+\infty)}_{\text{${\bm{t}}$-independent}}\text{
if }t_{1}\leq\cdots\leq t_{K-1}$ (4)

with the functions $R_{k}$, $k\in[K-1]$ defined by

$R_{k}(t):=\frac{1}{n}\sum_{i=1}^{n}\ell(h_{\rm
thr}(\hat{a}({\bm{x}}_{i});(\cdots,-\infty,\underbrace{t}_{\text{$k$-th}},+\infty,\cdots)),y_{i}),$
(5)

which is the empirical task risk when it labels $\hat{a}({\bm{x}}_{i})<t$ as
$k$ and $\hat{a}({\bm{x}}_{i})\geq t$ as $(k+1)$. Figure 1 schematizes
Equation (4): the summation of the red parts implies the left-hand side term
of (4), and the summation of the blue parts implies the latter term of the
right-hand side of (4). This relation implies the equivalence between
minimizing the empirical task risk $\frac{1}{n}\sum_{i=1}^{n}\ell(h_{\rm
thr}(\hat{a}({\bm{x}}_{i});{\bm{t}}),y_{i})$ and minimizing $R_{k}(t_{k})$,
$k=1,\ldots,K-1$ independently for each $k$ when the threshold parameters
satisfy the order condition $t_{1}\leq\cdots\leq t_{K-1}$. Therefore,
Algorithm 3 solves the latter subproblems associated with each threshold
parameter independently.

Input: $\hat{a}_{j}$ and ${\mathcal{Y}}_{j}$ ($j=1,\ldots,N$), loss matrix
$M$, and candidate vector ${\bm{c}}$ prepared by Algorithm 1.

/* Parallelizable with $k\in[K-1]$. */

1 for  _$k=1,\ldots,K-1$_ do

/* Calculate $(R_{j,k})_{j\in[N+1]}$. Further parallelizable (see Appendix A).
*/

2 $R_{1,k}=0$;

3 for  _$j=1,\ldots,N$_ do $R_{j+1,k}=R_{j,k}+M_{j,k}-M_{j,k+1}$;

/* Calculate threshold parameter $\hat{t}_{k}$. */

4 $\hat{t}_{k}=c_{j_{k}}$ with
$j_{k}=\min(\operatorname{{arg\,min}}(R_{j,k})_{j\in[N+1]})$;

Output: A threshold parameter vector $(\hat{t}_{k})_{k\in[K-1]}$.

Algorithm 3 IO algorithm to calculate the optimal threshold parameter vector
Figure 1: This figure assists in understanding Equation (4). In calculation of
$R_{k}(t_{k})$, the learned 1DT value $\hat{a}({\bm{x}}_{i})$ located in the
range marked with $l$ is labeled as $l$, and the corresponding task loss is
$\ell(l,y_{i})$.

As the condition of Equation (4) also suggests, the threshold parameter vector
$\hat{{\bm{t}}}$ obtained by Algorithms 1 and 3 becomes optimal in
minimization of the empirical task risk as long as the learned vector
$\hat{{\bm{t}}}$ eventually follows the appropriate order
$\hat{t}_{1}\leq\cdots\leq\hat{t}_{K-1}$:

######  Theorem 1.

For any data $\\{({\bm{x}}_{i},y_{i})\\}_{i\in[n]}$, learned 1DT $\hat{a}$,
and task loss $\ell$, the threshold parameter vector ${\bm{t}}=\hat{{\bm{t}}}$
obtained by Algorithms 1 and 3 minimizes the empirical task risk
$\frac{1}{n}\sum_{i=1}^{n}\ell(h_{\rm
thr}(\hat{a}({\bm{x}}_{i});{\bm{t}}),y_{i})$, if it satisfies the order
condition $\hat{t}_{1}\leq\cdots\leq\hat{t}_{K-1}$.

###### Proof of Theorem 1.

In this proof, we use the notations $\hat{a}_{j}$, $R_{j,k}$, and $j_{k}$,
defined in Algorithms 1 and 3 as well. Because of the definition, it holds
that

$R_{j,k}=n\\{R_{k}(c_{j})-R_{k}(c_{1})\\}\text{~{}for
all~{}}j\in[N+1],\,k\in[K-1].$ (6)

According to Equations (4) and (6), one has that

$\begin{split}&\frac{1}{n}\sum_{i=1}^{n}\ell(h_{\rm
thr}(\hat{a}({\bm{x}}_{i});(c_{i_{k}})_{k\in[K-1]},y_{i})\\\
&=\sum_{k=1}^{K-1}R_{k}(c_{i_{k}})-\sum_{k=2}^{K-1}R_{k}(+\infty)\quad(\because\eqref{eq:IOTEQ})\\\
&=\sum_{k=1}^{K-1}\biggl{\\{}\frac{1}{n}R_{i_{k},k}+R_{k}(c_{1})\biggr{\\}}-\sum_{k=2}^{K-1}R_{k}(+\infty)\quad(\because\eqref{eq:IOTEQ2})\\\
&=\frac{1}{n}\sum_{k=1}^{K-1}R_{i_{k},k}+\text{($(i_{k})_{k\in[K-1]}$-independent
term)}\end{split}$ (7)

with indices $i_{1},\ldots,i_{K-1}\in[N+1]$ satisfying $i_{1}\leq\cdots\leq
i_{K-1}$. Minimization of $\sum_{k=1}^{K-1}R_{i_{k},k}$ regarding
$i_{1},\ldots,i_{K-1}$ amounts to minimization of the empirical task risk with
the threshold parameters $(c_{i_{k}})_{k\in[K-1]}$ regarding
$i_{1},\ldots,i_{K-1}$ as far as the solution of the former problem keeps the
ascending order. Algorithm 3 finds the index $j_{k}$ that minimizes
$(R_{j,k})_{j\in[N+1]}$ for each $k\in[K-1]$. Namely, $(j_{k})_{k\in[K-1]}$ is
a minimizer of $\sum_{k=1}^{K-1}R_{i_{k},k}$ regarding $i_{1},\ldots,i_{K-1}$.
Therefore, under the assumption of this theorem
($\hat{t}_{1}\leq\cdots\leq\hat{t}_{K-1}$ or equivalently $j_{1}\leq\cdots\leq
j_{K-1}$), the threshold parameter vector
$\hat{{\bm{t}}}=(c_{j_{k}})_{k\in[K-1]}$ minimizes the empirical task risk. ∎

Furthermore, we found that the assumption of Theorem 1 holds for the task with
a convex task loss, which is defined by being convex regarding the first
augment, for example, $\ell=\ell_{\rm ab},\ell_{\rm sq}$:

######  Theorem 2.

For any data $\\{({\bm{x}}_{i},y_{i})\\}_{i\in[n]}$ and learned 1DT $\hat{a}$,
the threshold parameter vector ${\bm{t}}=\hat{{\bm{t}}}$ obtained by
Algorithms 1 and 3 satisfies the order condition
$\hat{t}_{1}\leq\cdots\leq\hat{t}_{K-1}$ and hence minimizes the empirical
task risk $\frac{1}{n}\sum_{i=1}^{n}\ell(h_{\rm
thr}(\hat{a}({\bm{x}}_{i});{\bm{t}}),y_{i})$, if the task loss function $\ell$
satisfies

$\ell(k,l)-2\ell(k+1,l)+\ell(k+2,l)\geq 0\text{~{}for
all~{}}k\in[K-2],\,l\in[K].$ (8)

###### Proof of Theorem 2.

In this proof, we use the notations $\hat{a}_{j}$, ${\mathcal{Y}}_{j}$,
$M_{j,k}$, $R_{j,k}$, and $j_{k}$, defined in Algorithms 1 and 3 as well. We
prove this theorem using the proof by contradiction. Assume
$\hat{t}_{k}>\hat{t}_{k+1}$ (i.e., $j_{k}>j_{k+1}$) with some $k\in[K-2]$. One
has that

$\begin{split}&\\{R_{j_{k},k}+R_{j_{k+1},k+1}\\}-\\{R_{j_{k+1},k}+R_{j_{k},k+1}\\}\\\
&=\biggl{\\{}\sum_{j<j_{k}}(M_{j,k}-M_{j,k+1})+\sum_{j<j_{k+1}}(M_{j,k+1}-M_{j,k+2})\biggr{\\}}\\\
&-\biggl{\\{}\sum_{j<j_{k+1}}(M_{j,k}-M_{j,k+1})+\sum_{j<j_{k}}(M_{j,k+1}-M_{j,k+2})\biggr{\\}}\quad(\because\text{Algorithm~{}\ref{alg:IO-
algorithm}, Line~{}\ref{IO3}})\\\
&=\sum_{j<j_{k}}(M_{j,k}-2M_{j,k+1}+M_{j,k+2})-\sum_{j<j_{k+1}}(M_{j,k}-2M_{j,k+1}+M_{j,k+2})\\\
&=\sum_{j\in\\{j_{k+1},\ldots,j_{k}-1\\}}(M_{j,k}-2M_{j,k+1}+M_{j,k+2})\\\
&=\sum_{j\in\\{j_{k+1},\ldots,j_{k}-1\\}}\sum_{y_{i}\in{\mathcal{Y}}_{j}}\\{\ell(k,y_{i})-2\ell(k+1,y_{i})+\ell(k+2,y_{i})\\}\quad(\because\text{Algorithm~{}\ref{alg:Preparation},
Line~{}\ref{P4}})\\\ &\geq 0\quad(\because\eqref{eq:Convex}).\end{split}$ (9)

When the equality holds in the last inequality of (9),
$R_{j_{k},k}+R_{j_{k+1},k+1}=R_{j_{k+1},k}+R_{j_{k},k+1}$ also holds.
Furthermore, $R_{j_{k},k}\leq R_{j_{k+1},k}$ and $R_{j_{k+1},k+1}\leq
R_{j_{k},k+1}$ hold because $j_{k}$ minimizes $(R_{j,k})_{j\in[N+1]}$ and
$j_{k+1}$ minimizes $(R_{j,k+1})_{j\in[N+1]}$. From these results, it
necessarily holds that $R_{j_{k},k}=R_{j_{k+1},k}$ and
$R_{j_{k+1},k+1}=R_{j_{k},k+1}$, which contradicts the definition
$j_{k}=\min(\operatorname{{arg\,min}}(R_{j,k})_{j\in[N+1]})=\min(\\{j_{k+1},j_{k},\ldots\\})$
since $j_{k}>j_{k+1}$. When the equality does not hold in the last inequality
of (9), $R_{j_{k},k}+R_{j_{k+1},k+1}>R_{j_{k+1},k}+R_{j_{k},k+1}$ holds, and
hence the threshold parameter vector $(\check{t}_{k})_{k\in[K-1]}$ consisting
of $\check{t}_{k}=\hat{t}_{k+1}$, $\check{t}_{k+1}=\hat{t}_{k}$, and
$\check{t}_{j}=\hat{t}_{j}$ for $j\neq k,k+1$ is better than $\hat{{\bm{t}}}$
in minimization of $\sum_{k=1}^{K-1}R_{i_{k},k}$ regarding
$i_{1},\ldots,i_{K-1}$. According to the sub-optimality of $\hat{{\bm{t}}}$,
it can be seen that Algorithm 3 does not output $\hat{{\bm{t}}}$, which is a
contradiction to the setting of this proof. Therefore, the proof is concluded.
∎

On the ground of Theorems 1 and 2, it is guaranteed that one can obtain the
optimal threshold parameter vector efficiently by the _parallelized IO_
(_PIO_) _algorithm_ (Algorithm 3, or Algorithm 5 in Appendix A) for the task
with a convex task loss function.

##  5 Experiment

In order to verify that the proposed PIO algorithm can be faster than the
existing DP algorithm, we performed a practical OR experiment using real-world
datasets to compare the computational efficiency of the DP, (non-parallelized)
IO, and PIO algorithms.

We addressed an OR task of estimating age from facial image by a discrete-
value prediction, with reference to experiment settings in the previous OR
studies (Cao et al., 2020; Yamasaki, 2023). The used datasets are MORPH-2
(Ricanek and Tesafaye, 2006), CACD (Chen et al., 2014), and AFAD (Niu et al.,
2016). After pre-processing as that in Cao et al. (2020); Yamasaki (2023), we
used data of the total data size $n_{\rm tot}=55013,159402,164418$ and
$K=55,49,26$ in these datasets. For each dataset, we resized all images to
$128\times 128$ pixels of 3 RGB channels and divided the dataset into a 90 %
training set and a 10 % validation set. As a learning method of the 1DT, we
tried OLR that we reviewed in Section 2. We implemented the 1DT with ResNet18
or ResNet34 (He et al., 2016), and trained the network using the negative-log-
likelihood loss function and Adam with the mini-batch size 256 for 100 epochs.
For the task with the absolute task loss $\ell=\ell_{\rm ab}$, we optimized
the threshold parameters by the DP, IO, or PIO algorithms and calculated the
empirical task risk for the validation set, validation error (VE), at every
epoch. We used images randomly cropped to $120\times 120$ pixels of 3 channels
in ‘training of 1DT’, and images center-cropped to the same size in
‘preparation (Algorithm 1)’, ‘optimization of threshold parameters’, and
‘calculation of VE’, as inputs.

We experimented with Python 3.8.13 and PyTorch 1.10.1 under the computational
environment with 4 CPUs Intel Xeon Silver 4210 and a single GPU GeForce RTX
A6000. We performed ‘training of 1DT’, ‘Line 1 of preparation’, and
‘calculation of VE’ using the GPU, and ‘remaining lines of preparation’ and
‘optimization of threshold parameters’ using the CPUs, which was the best
allocation of the CPUs and GPU among those we tried. We employed the
parallelization with DataLoader of num_workers 8 for the calculations on the
GPU, $K$-parallelization for ‘Line 1 of preparation’ and
$(K-1)$-parallelization for ‘Lines 3–3 of Algorithm 3’ of the PIO among the
calculations on the CPUs. See https://github.com/yamasakiryoya/ECOTL for
program codes we used.

Table 1: Calculation time for 100 epochs (in minutes).

|  |  MORPH-2 |  CACD |  AFAD
---|---|---|---|---
Training of 1DT | ‘ResNet18, GPU’ or |  32.00 |  86.59 |  83.46
‘ResNet34, GPU’ |  45.27 |  125.27 |  121.54
| ‘Line 1, ResNet18, GPU’ or |  23.22 |  48.11 |  58.30
Preparation | ‘Line 1, ResNet34, GPU’, |  21.33 |  50.61 |  59.46
(Algorithm 1) | ‘Lines 1 and 1, CPUs’, and |  1.05 |  3.61 |  3.63
| ‘Lines 1 and 1, CPUs’ |  9.66 |  13.60 |  9.94
Optimization of | ‘DP, CPUs’, |  81.04 |  186.60 |  130.35
threshold parameters | ‘IO, CPUs’, or |  87.71 |  224.09 |  123.55
‘PIO, CPUs’ |  18.45 |  24.38 |  15.26
Calculation of VE | ‘ResNet18, GPU’ or |  5.20 |  8.40 |  9.68
‘ResNet34, GPU’ |  5.07 |  9.18 |  9.85

Table 1 shows the calculation time taken for each process. The total
calculation time with the DP, IO, and PIO algorithms ($t_{\rm DP}$, $t_{\rm
IO}$, and $t_{\rm IO}$) had the following relationship: the ratios $(t_{\rm
IO}/t_{\rm DP},t_{\rm PIO}/t_{\rm DP})$ are $(1.04,0.59)$ for MORPH-2,
$(1.11,0.53)$ for CACD, and $(0.98,0.61)$ for AFAD for the ResNet18-based 1DT;
$(1.04,0.62)$ for MORPH-2, $(1.10,0.58)$ for CACD, and $(0.98,0.66)$ for AFAD
for the ResNet34-based 1DT. Numerical evaluations related to the computation
time may change slightly depending on a used computational environment, but we
consider that the following facts will hold universally: the DP and IO
algorithms have comparable computational efficiency, while the PIO algorithm
can be faster than the former two algorithms.

##  6 Conclusion

In this paper, we proposed the parallelizable IO algorithm (Algorithm 3) for
acquiring the optimal threshold labeling (Yamasaki, 2023), and derived
sufficient conditions (Theorems 1 and 2) for this algorithm to successfully
output the optimal threshold labeling. Our numerical experiment showed that
the computation time for the whole learning process of a threshold method with
the optimal threshold labeling was reduced to approximately 60 % under our
computational environment, by using the parallelized IO algorithm compared to
using the existing DP algorithm. On the ground of these results, we conclude
that the proposed algorithm can be useful to accelerate the application of
threshold methods especially for cost-sensitive tasks with a convex task loss,
which are often tackled in OR research.

## References

  * Bürkner and Vuorre [2019] P.-C. Bürkner and M. Vuorre.  Ordinal regression models in psychology: A tutorial.  _Advances in Methods and Practices in Psychological Science_ , 2(1):77–101, 2019.
  * Cao et al. [2020] W. Cao, V. Mirjalili, and S. Raschka.  Rank consistent ordinal regression for neural networks with application to age estimation.  _Pattern Recognition Letters_ , 140:325–331, 2020.
  * Chen et al. [2014] B.-C. Chen, C.-S. Chen, and W. H. Hsu.  Cross-age reference coding for age-invariant face recognition and retrieval.  In _Proceedings of the European Conference on Computer Vision_ , pages 768–783, 2014.
  * Chu and Keerthi [2005] W. Chu and S. S. Keerthi.  New approaches to support vector ordinal regression.  In _Proceedings of the International Conference on Machine Learning_ , pages 145–152, 2005.
  * Chu and Keerthi [2007] W. Chu and S. S. Keerthi.  Support vector ordinal regression.  _Neural Computation_ , 19(3):792–815, 2007.
  * Cohen [1960] J. Cohen.  A coefficient of agreement for nominal scales.  _Educational and Psychological Measurement_ , 20(1):37–46, 1960.
  * Cohen [1968] J. Cohen.  Weighted kappa: nominal scale agreement provision for scaled disagreement or partial credit.  _Psychological Bulletin_ , 70(4):213, 1968.
  * da Costa et al. [2008] J. F. P. da Costa, H. Alonso, and J. S. Cardoso.  The unimodal model for the classification of ordinal data.  _Neural Networks_ , 21(1):78–91, 2008.
  * He et al. [2016] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 770–778, 2016.
  * Li and Lin [2007] L. Li and H.-T. Lin.  Ordinal regression by extended binary classification.  In _Advances in Neural Information Processing Systems_ , pages 865–872, 2007.
  * Lin and Li [2006] H.-T. Lin and L. Li.  Large-margin thresholded ensembles for ordinal regression: Theory and practice.  In _Algorithmic Learning Theory_ , pages 319–333, 2006.
  * Lin and Li [2012] H.-T. Lin and L. Li.  Reduction from cost-sensitive ordinal ranking to weighted binary classification.  _Neural Computation_ , 24(5):1329–1367, 2012\.
  * Liu [2011] T.-Y. Liu.  _Learning to Rank for Information Retrieval_.  Springer Science & Business Media, 2011.
  * McCullagh [1980] P. McCullagh.  Regression models for ordinal data.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 42(2):109–127, 1980.
  * Niu et al. [2016] Z. Niu, M. Zhou, L. Wang, X. Gao, and G. Hua.  Ordinal regression with multiple output cnn for age estimation.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 4920–4928, 2016.
  * Pedregosa et al. [2017] F. Pedregosa, F. Bach, and A. Gramfort.  On the consistency of ordinal regression methods.  _Journal of Machine Learning Research_ , 18(Jan):1769–1803, 2017.
  * Prechelt [2002] L. Prechelt.  Early stopping - but when?  In _Neural Networks: Tricks of the Trade_ , pages 55–69. Springer, 2002.
  * Ricanek and Tesafaye [2006] K. Ricanek and T. Tesafaye.  Morph: A longitudinal image database of normal adult age-progression.  In _Proceedings of the IEEE International Conference on Automatic Face and Gesture Recognition_ , pages 341–345, 2006.
  * Shashua and Levin [2002] A. Shashua and A. Levin.  Taxonomy of large margin principle algorithms for ordinal regression problems.  In _Advances in Neural Information Processing Systems_ , pages 937–944, 2002.
  * Shashua and Levin [2003] A. Shashua and A. Levin.  Ranking with large margin principle: Two approaches.  In _Advances in Neural Information Processing Systems_ , pages 961–968, 2003.
  * Yamasaki [2022] R. Yamasaki.  Unimodal likelihood models for ordinal data.  _Transactions on Machine Learning Research_ , pages 1–25, 2022.  URL https://openreview.net/forum?id=1l0sClLiPc.
  * Yamasaki [2023] R. Yamasaki.  Optimal threshold labeling for ordinal regression methods.  _Transactions on Machine Learning Research_ , pages 1–42, 2023.  URL https://openreview.net/forum?id=mHSAy1n65Z.
  * Yu et al. [2006] S. Yu, K. Yu, V. Tresp, and H.-P. Kriegel.  Collaborative ordinal regression.  In _Proceedings of the International Conference on Machine Learning_ , pages 1089–1096, 2006.

##  Appendix A Further Parallelization of IO Algorithm

The IO algorithm (Algorithm 3) can be parallelized further. Algorithm 5
describes a parallelization procedure regarding $k\in[K-1]$ and $l\in[L^{\rm
par}]$ with $L^{\rm par}=\lceil(N+1)/L\rceil$ for a user-specified integer
$L$. This algorithm calculates the translated version of $L^{\rm par}$ blocks
splitting $(R_{j,k})_{j\in[N+1]}$ (Lines 4–4) and then modifies the
translation of each block (Lines 4–4) in parallel, to calculate the whole of
$(R_{j,k})_{j\in[N+1]}$, instead of Lines 3 and 3 of Algorithm 3, for each
$k\in[K-1]$. Algorithms 3 and 5 yield the same outputs. If one has sufficient
computing resources, it would be better to set $L$ to an integer of the order
$O(\sqrt{N})$ theoretically since this parallelized IO algorithm takes
computation time of the order $O(L+L^{\rm par})$, and this implementation
improves the order of computation time from $O(N\cdot K)$ of the DP algorithm
(Algorithms 2) to $O(\sqrt{N})$. However, we note that it is not realistic to
expect as dramatic an improvement in actual computation time as the order
evaluation suggests, owing to the limitation of computing resources: under the
setting of the experiment in Section 5, Algorithm 5 with
$L=\lceil(N+1)/2\rceil,\lceil(N+1)/3\rceil,\lceil(N+1)/4\rceil$ (such that
$L^{\rm par}=2,3,4$) and CPUs takes $37.27,52.07,67.08$ for MORPH-2,
$40.38,51.98,64.95$ for CACD, and $23.33,28.92,35.43$ for AFAD (in minutes)
for ‘optimization of threshold parameters’ for 100 epochs. ‘PIO, CPUs’
following Algorithm 3 was better than Algorithm 5 under our computational
environment.

Input: $\hat{a}_{j}$ and ${\mathcal{Y}}_{j}$ ($j=1,\ldots,N$), loss matrix
$M$, candidate vector ${\bm{c}}$ prepared by Algorithm 1, $L\in{\mathbb{N}}$,
$L^{\rm par}:=\lceil(N+1)/L\rceil$, $L^{\rm quo}:=\lfloor(N+1)/L\rfloor$,
$L^{\rm rem}:=(N+1)-L\cdot L^{\rm quo}$, and $L_{l}:=(l-1)L$.

/* Parallelizable with $k\in[K-1]$. */

1 for  _$k=1,\ldots,K-1$_ do

/* Calculate $(Q_{j,k})_{j\in[N+1]}$. Parallelizable with $l\in[L^{\rm par}]$.
*/

2 for  _$l=1,\ldots,L^{\rm par}$_ do

3 $Q_{L_{l}+1,k}=0$;

4 if  _$l\leq L^{\rm quo}$_ then for  _$j=1,\ldots,L-1$_ do
$Q_{L_{l}+j+1,k}=Q_{L_{l}+j,k}+M_{L_{l}+j,k}-M_{L_{l}+j,k+1}$;

5 else for  _$j=1,\ldots,L^{\rm rem}-1$_ do
$Q_{L_{l}+j+1,k}=Q_{L_{l}+j,k}+M_{L_{l}+j,k}-M_{L_{l}+j,k+1}$;

/* Calculate $(R_{j,k})_{j\in[N+1]}$. Parallelizable with $l\in[L^{\rm par}]$.
*/

6 for  _$l=1,\ldots,L^{\rm par}$_ do

7 $S_{k,l}=\operatorname{{\mathbbm{1}}}(l\neq
1)\sum_{m\in[l-1]}(Q_{L_{m+1},k}-Q_{L_{m}+1,k}+M_{L_{m+1},k}-M_{L_{m+1},k+1})$;

8 if  _$l\leq L^{\rm quo}$_ then for  _$j=1,\ldots,L$_ do
$R_{L_{l}+j,k}=Q_{L_{l}+j,k}+S_{k,l}$;

9 else for  _$j=1,\ldots,L^{\rm rem}$_ do
$R_{L_{l}+j,k}=Q_{L_{l}+j,k}+S_{k,l}$;

/* Calculate threshold parameter $\hat{t}_{k}$. */

10 $\hat{t}_{k}=c_{j_{k}}$ with
$j_{k}=\min(\operatorname{{arg\,min}}(R_{j,k})_{j\in[N+1]})$;

Output: A threshold parameter vector $(\hat{t}_{k})_{k\in[K-1]}$.

Algorithm 4 Further parallelized IO algorithm to calculate the optimal
threshold parameter vector555 $\lceil\cdot\rceil$ is the ceiling function, and
$\lfloor\cdot\rfloor$ is the floor function.