The problem of constructing the best fitted monotone regression is NP-hard problem and can be formulated in the form of a convex programming problem with linear constraints. The paper proposes a simple greedy algorithm for finding a sparse monotone regression using Frank-Wolfe-type approach. A software package for this problem is developed and implemented in R and C++. The proposed method is compared with the well-known pool-adjacent-violators algorithm (PAVA) using simulated data.
The recent years have seen an increasing interest in shape-constrained estimation
in statistics [
Using mathematical programming approach the works [
Monotone regression is widely used in mathematical statistics [
Constructing monotone regression we assume a relationship between a predictor x = (x1; : : : ; xn) and a response y = (y1; : : : ; yn). In the general case it is expected that xi+1 xi 6= const, xi < xi+1, i = 1; : : : ; n 1. ? The work was supported by RFBR (grant 16-01-00507).
A sequence z = (z1; : : : ; zn) 2 Rn is called monotone if zi zi 1
0; i = 2; : : : ; n: Denote 1n the set of all vectors from Rn, which are monotone.
The problem of constructing monotone regression can be formulated in the form of a convex programming problem with linear constraints as follows: it is necessary to find a vector z 2 Rn with the lowest mean square error of approximation to the given vector y 2 Rn under condition of monotonicity of z: n f (z) = 1 X(zi n i=1 yi)2
min ; ! z2 1n (1)
In many situations researchers have no information regarding the
mathematical specification of the true regression function. Typically, this involves
non-decreasing of yi’s with the ordered xi’s. Such a situation is called isotonic
regression. Isotonic regression (monotone regression) is a special case to the
kmonotone regression [
It is well-known that the problem (1) is NP-hard problem [
For the convenience of solving the problem (1), we move from points zi to its increments i, where i = zi+1 zi, i = 1; : : : ; n 1, 0 = z1. Then monotonicity of z corresponds non-negativity of i’s (exept 0). The proposed method is compared with the well-known pool-adjacent-violators algorithm (PAVA) using simulated data. 2 2.1
Algorithms for monotone regression
Simple iterative algorithm for solving the problem (1) is called
Pool-AdjacentViolators Algorithm (PAVA) [
PAVA computes a non-decreasing sequence of values z = (zi)in=1 such that
the problem (1) is optimized. In the simple monotone regression case we have
the measurement pairs (xi; yi). Let us assume that these pairs are ordered with
respect to the predictors. The following (Algorithm 1) is a pseudocode of PAVA
for the problem. The generalized pool-adjacent-violators algorithm (GPAVA),
which is a strict generalization of PAVA, was developed in the article [
The block values are expanded with respect to the observations i = 1; : : : ; n
such that the final result is the vector z of length n with elements zi of increasing
e
order [
Algorithm 1: Pool-Adjacent-Violators Algorithm (PAVA) begin
Let zj(0) := yi be the start point, l = 0; The index for the blocks is r = 1; : : : ; B where at step l = 0 we set B := n, i.e. each observation zr(0) forms a block; repeat (Adjacent pooling ) Merge values of z(l) into blocks if zr+1 < zr(l); (l) Solve f (z) for each block r, i.e., compute the update based on the solver which gives zr(l+1) := s(zr(l)), the solver s is conditional (weighted) mean and (weighted) quantities; If zr(l+)1 zr(l) then l := l + 1;
(l) until the z-blocks are increasing, i.e. zr+1
Return z; zr(l) for all r; end 2.2
Frank–Wolfe method (or conditional gradient method) solves conditional convex
optimization problems in vector finite-dimensional space. The method was
introduced in 1956. The original algorithm did not use a fixed step size, and has the
complexity of the linear programming. Frank–Wolfe method was developed by
Levitin and Polyak in 1966, and V.F. Demianov and A.M. Rubinov generalized
it to the case of arbitrary Banach spaces in 1970 [
As it was mentioned above, for computational convenience of the problem (1), we moved from points zi to increments i = zi+1 zi, i = 1; : : : ; n 1, 0 = z1. Then the problem (1) can be rewritten as follows:
n g( ) := 1 X n i=1 i 1 X j=0 j yi !2 ! min; 2S (2) where S denotes the set of all = ( 0; 1; : : : ; n 1) 2 Rn such that 0 2 R, ( 1; : : : ; n 1) 2 Rn+ 1 and Pjn=01 j maxi yi.
Let rg( ) denote the gradient of function g at point .
It should be noted that for larger-scale problems the solution can appear computationally quite challenging. In this regard, the present study proposes to use a greedy algorithm of Frank–Wolfe type for solving this problem.
The following (Algorithm 2) is a pseudocode of Frank–Wolfe-type algorithm for the problem (2).
The rate of convergence is estimated according to the following theorem. Algorithm 2: Greedy algorithm for sparse monotone regression begin
Let N be the number of iteration; Our function g( ) and the feasible set S were defined above.
@g @g @g
be the gradient of function g at Let rg( ) = i=k j=0 Let zero vector 0 = (0; : : : ; 0) be the start point, and let the counter t = 0; while t < N do
Calculate rg( t), the gradient of the function g at the point t; Let et be the solution of the linear optimization problem hrg( t)T ; i ! min; where hrg( t)T ; i is a scalar product of vectors; 2S (Update step) Let t+1 = t + t(et
t), t = t+22 and than t := t + 1; Recover the monotone sequence z = (z1; : : : ; zn) from the vector of increments N ; end Theorem 1. Let f tg be generated according to the Frank–Wolfe method (Algorithm 2) using the step-size rule t = t+22 : Then for all t 2 g( t) g 4 r n(n + 1)(2n + 1) (maxi yi 6n2 t + 2 mini yi)2 ; (3) where g is the optimal solution of (2).
Proof. It it is know [
Let r2g( ) := : It is well-known that if rg is differentiable then its Lipschitz constant L satisfies the inequality
L
sup kr2g( )k2:
L
vun 1 sup tuX k=0 It is easy to prove and Diam(S) := p2(maxi yi mini yi).
The disadvantage of this method is the dependence of the theoretical degree
of convergence on the dimensionality of the problem. The papers [
The algorithms have been implemented both in R and C++. We compared the performance of the greedy algorithm (Algorithm 2) with the performance of PAVA (Algorithm 1) using simulated data sets.
It should also be noted that the PAVA’s speed is significantly higher for small-scale tasks in R. But if the number of points is greater than at least 2000, the greedy algorithm spends less time searching for a solution (Fig. 1).
Tables 1, 2 present empirical results for PAVA and greedy algorithms for a simulated set of points. The simulated points are obtained as the values of logarithm function with added normally distributed noise: A = f(xi; yi); yi = ln(x0 + i4x) + 'i; 'i N (0; 1); x0 = 1; 4x = 1; i = 1; : : : ; 10000: The dimension of the problem is 10000 points. The tables contain information on errors 1 Pn n i=1(zi yi)2, elapsed time, cardinality and greedy algorithm’s iteration number.
The results show that error of greedy algorithm are getting closer to the error of PAVA with increase of number of iterations for greedy algorithm. While PAVA is better than greedy algorithm in terms of errors, the solutions of greedy algorithm have a better sparsity. Greedy algorithm’s output solution is more sparse. It should be noted that the elapsed time for PAVA implemented in C++ is smaller than for greedy algorithm. However, greedy algorithm has a better rate of convergence if number of iterations is less than 700 for the algorithms implemented in R,. Both algorithms obtain a sparse solutions, but we can control the number of nonzero elements (cardinality) in the greedy algorithm as opposed to PAVA. Generally, the greedy algorithm’s cardinality increases by one at each iteration. Consequently, we should limit the number of iterations to obtain more sparse solution.
Figure 2 shows simulated points (N = 100) with logarithm structure and isotonic regressions, where green line represents the greedy algorithm’s isotonic e m i T 1:5
1 0:5 0 1000 2000
3000 Points number 4000 regression and red line presents PAVA’s isotonic regression. Greedy algorithm gives a solution with 14 jumps, and PAVA with 16 jumps. Since the solutions of the greedy algorithm are more sparse, the greedy algorithm error (") is slightly higher than the PAVA.
The obtained empirical results for the greedy algorithm show that the degree of convergence for the considered examples is much higher than its theoretical estimates obtained in Theorem 1. Our research proposes an algorithm for solving the problem of constructing the best fitted monotone regression by using the Frank–Wolfe method. The software was implemented in R and C++. We compared the performance of the greedy algorithm with the performance of PAVA using simulated data sets. While PAVA gives a slightly smaller errors than greedy algorithm, greedy algorithm obtains significantly sparser solutions. The advantages of greedy algorithm are the simplicity of implementation, the potential for controlling cardinality and the elapsed time is lower for the implementation in R in the case of problem with large dimension.