Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
13,797 characters · 6 sections · 22 citation commands
Bivariate Isotonic Regression by Dynamic Programming
\selectlanguage{english}
Monotone or isotonic regression is a non-parametric method introduced by Brunk1955 and extensively covered by Barlow1972. It has a wide range of applications, namely in economics, epidemiology, and biometrics. The univariate problem ($d=1$) is relatively straightforward: approximate a given sequence $y = (y_1, \cdots, y_n)$ partially ordered accordingly with another sequence $x = (x_1, \cdots, x_n)$ by an increasing sequence
minimising some risk function $R(z)$ that reflects the average error or loss in predicting $y$ with $z = (z_1, \cdots, z_n)$. If the input sequence $y$ has decreasing sections or blocks, the optimal values $z_k$ are the median or average of the corresponding elements $y_k$ with loss $\ell_1$ or $\ell_2$, respectively. Thus, isotonic regression typically produces a piecewise constant estimator Lim2025. Efficient algorithms are readily available to estimate the univariate isotonic regression, namely the pool-adjacent violators algorithm (PAVA).
The multivariate case ($d>1$) is more challenging, because the monotonicity of $z$ must be ensured for all non-decreasing sequences of the independent variables, that is,
where the partial order $\preceq$ is the standard Euclidean one, that is, $x_k \preceq x_l$ if and only if $x_{k1} \le x_{l1}, x_{k2} \le x_{l2}, \dots, x_{kd} \le x_{ld}$ Luss2024. Thus, the isotonic solution must satisfy a set of isotonicity constraints which are indexed by the set $\Gamma = \{ (k,l): x_k \preceq x_l \}$.
Most of the literature deals with bivariate isotonic regression on a grid previously ordered in both independent variables. For example, the classic algorithm of Dykstra1982, implemented in FORTRAN by Bril1984, applies PAVA successively to rows and columns until convergence.
General algorithms for multivariate isotonic regression is an open problem for future research Lim2025, with limited methods available. For example, the algorithm of Cheung2023 estimates the underlying mean response for multivariate features using an inverse projective Bayes approach. Their estimator is obtained by inverting the projective Bayes classifier across a grid of thresholds, yielding a monotone non-parametric estimate of the mean response that is non-decreasing in each feature, and a Normal-Inverse-Chi-Squared conjugate model is used to compute posterior probabilities at each threshold. Another general method is the isotonic recursive partitioning proposed by Luss2024 that follows a logic similar to classification and regression trees Breiman1984.
In this article, we extend the dynamic programming framework introduced by Rote2019 from the univariate to the bivariate isotonic problem ($d=2$), using an anti-diagonal traversal procedure. An implementation in R of the new algorithm is available on GitHub.
The isotonic problem can be defined sequentially as
where $L(z_k)$ is a loss function that measures the distance between $y_k$ and $z_k$ and $\Gamma$ is the set of isotonicity constraints. As suggested by Rote2019, this problem has a recursive functional representation of the form
where $\tilde{z}$ is the control value to be estimated in step $k$, given the state value $z$ previously estimated in $k+1$ for $k = n-1, n-2, \dots, 1$. $J(z)$ is the cost-to-go function of the state $z$, using the terminology and notation of Bertsekas1996.
Rote2019 proposed a simple algorithm with an overall run time of $O(n \log n)$ to solve this problem in the univariate case ($d=1$) using an unweighted $\ell_1$ loss function \footnote{In fact, Rote2019 also proposed a more complex algorithm for weighted loss $L(z) = w \left| z - y \right |$ with $w > 0$.}
His algorithm is based on the convexity of the cost-to-go function $J$ that first decreases monotonically to a minimum (not necessarily unique) and then increases monotonically. Thus, the goal of the first part of the algorithm ((ref)) is to find each minimum $p_k$ given $y_k$ for $k = 1, \dots, n$. For that purpose, it uses a breakpoint queue where each priority is determined by position $y_k$ and the value is 2 for new breakpoints and 1 elsewhere. The minimum is found by removing (eventual) increasing pieces from the top of the queue. The algorithm then computes the optimal solution by backward induction, noting that
The algorithm ((ref)) was made to deal with a dependent variable $y$ ordered on an independent variable $x$, but it can easily be extended to the case where $y$ is ordered on two independent variables. Thus, $y$ will live on a 2D grid ordered on both variables with $m$ rows and $n$ columns, and a generic cell $y_{ij}$ for $i = 1, \dots, m$ and $j = 1, \dots, n$.
As before, the aim of the first part of algorithm ((ref)) is to find the minima $p_{ij}$. The priority queue $Q$ is fed by an anti-diagonal traversal procedure, starting from the first anti-diagonal ($d=0$) which corresponds to the cell $y_{11}$, then to the second anti-diagonal ($d=1$) corresponding to $y_{12}$ and $y_{21}$, and so on, until the last anti-diagonal ($d=m+n-2)$ with $y_{mn}$. Cells on the same anti-diagonal share the same value of $i + j$, so $j = d - i$.
The second part of the algorithm is concerned with the computation of the optimal solution, noting that
It starts with the terminal condition $z_{mn} = p_{mn}$ and then ensures the monotonicity of the last column or the last row, following the rule ((ref)). Finally, it applies the rule ((ref)) using a backward anti-diagnonal traversal from $d = m+n-4$ to $0$.
The algorithm ((ref)) was implemented in R, the bivdp(y) function from the library libiso.R (v01.03) and is available in \url{https://github.com/pedroafonso1970/rlibs}. That library also includes the function isodp(x,y) with Rote's algorithm ((ref)) for the univariate case.
The algorithm ((ref)) was applied to the well-known baseball data set provided by Watnik1998 that describes the association of salary with a collection of player properties, including the number of runs batted and hits. The simple increasing linear ordering of the average salary in both variables is presented in table (ref).
The bivariate isotonic $\ell_1$ regression by dynamic programming produced a stepwise function of the salary on runs and hits very similar to the one derived with the classic $\ell_2$ algorithm of Dykstra1982 implemented by Bril1984, see figure (ref). The mean absolute error of the former (475.5) is also close to the last (443.0).
In many problems, the exact relationship between variables is unknown and difficult to specify parametrically Lim2025. Nevertheless, the structure of the relationship may be known or provided by some theory. In particular, this is a common situation in economics where demand (supply) typically decreases (increases) with prices, independently of the functional form. In these cases, the provided isotonic algorithm can be an elegant and flexible alternative to bivariate linear regression. It is also robust in the sense that it uses a $\ell_1$ norm.
The recursive nature of the new algorithm is also relevant, noting that dynamic programming has a wide range of applications in economics, including the savings problem, economic growth, job search, business cycles, oligopoly equilibrium, recursive contracts, and forecasting. Ljungqvist2018.
This work was financed by Fundação para a Ciência e a Tecnologia (FCT) under a doctorate auxiliary researcher grant at Universidade Católica Portuguesa (UCP) - Católica Lisbon Research Unit in Business & Economics (CUBE) with the reference 2023.15056.TENURE.067.