Introduction

The recent years have seen an increasing interest in shape-constrained estimation in statistics [10]. One of such problems is the problem of constructing monotone regression. The problem is to find best fitted non-decreasing points to a given set of points on the plane. The survey of results on monotone regression can be found in the book by Robertson and Dykstra [25]. The papers of Barlow and Brunk [3], Dykstra [16], Best and Chakravarti [4], Best [5] consider the problem of finding monotone regression in quadratic and convex programming frameworks.

Using mathematical programming approach the works [1,21,31] have recently provided some new results on the topic. The papers [7,17] extend the problem to particular orders defined by the variables of a multiple regression. The paper [8] investigates a dual active-set algorithm for regularized monotonic regression.

Monotone regression is widely used in mathematical statistics [2,10]; in smoothing of empirical data [15]; in shape-preserving approximation [19], [26], [30], [6], [27], [13]; in shape-preserving dynamic programming [9].

Constructing monotone regression we assume a relationship between a predictor x = (x 1 , . . . , x n ) and a response y = (y 1 , . . . , y n ). In the general case it is expected that

x i+1 − x i = const, x i < x i+1 , i = 1, . . . , n − 1.

The work was supported by RFBR (grant 16-01-00507).

A sequencez = (z 1 , . . . , z n ) ∈ R n is called monotone if z i − z i−1 ≥ 0, i = 2, . . . , n.

Denote ∆ n

1 the set of all vectors from R n , 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 ∈ R n with the lowest mean square error of approximation to the given vector y ∈ R n under condition of monotonicity of z:

f (z) = 1 n n i=1 (z i − y i ) 2 → min z∈∆ n 1 ,(1)

In many situations researchers have no information regarding the mathematical specification of the true regression function. Typically, this involves non-decreasing of y i 's with the ordered x i 's. Such a situation is called isotonic regression. Isotonic regression (monotone regression) is a special case to the kmonotone regression [24].

It is well-known that the problem (1) is NP-hard problem [24]. In this paper we present a simple greedy algorithm which employs Frank-Wolfe-type approach for finding sparse monotone regression. A software package for this problem is developed and implemented in R and C++.

For the convenience of solving the problem (1), we move from points z i to its increments ζ i , where

ζ i = z i+1 − z i , i = 1, . . . , n − 1, ζ 0 = z 1 .

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.

Algorithms for monotone regression

PAVA

Simple iterative algorithm for solving the problem (1) is called Pool-Adjacent-Violators Algorithm (PAVA) [11,24]. The work [4] examined the generalization of this algorithm. The paper [32] studied this problem as the problem of identifying the active set and proposed a direct algorithm of the same complexity as the PAVA (the dual algorithm).

PAVA computes a non-decreasing sequence of values z = (z i ) n i=1 such that the problem (1) is optimized. In the simple monotone regression case we have the measurement pairs (x i , y i ). 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 [33].

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 z i of increasing order [24].

:= s(z (l)

r ), the solver s is conditional (weighted) mean and (weighted) quantities;

• If z (l) r+1 ≤ z (l) r then l := l + 1; until the z-blocks are increasing, i.e. z (l) r+1 ≥ z (l)

r for all r; • Return z; end

Frank-Wolfe type greedy algorithm

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 [14]. Recently Frank-Wolfe type methods have caused an increased interest related to the possibility of obtaining sparse solutions, as well as a good scaling [12,23]. In particular, [22,34] researched algorithms for solving problems with penalty functions (instead of considering the conditional optimization problems). Besides, the paper [34] uses interlacing boosting with fixed-rank local optimization.

As it was mentioned above, for computational convenience of the problem (1), we moved from points z i to increments

ζ i = z i+1 − z i , i = 1, . . . , n − 1, ζ 0 = z 1 .

Then the problem (1) can be rewritten as follows:

g(ζ) := 1 n n i=1 i−1 j=0 ζ j − y i 2 → min ζ∈S ,(2)

where S denotes the set of all

ζ = (ζ 0 , ζ 1 , . . . , ζ n−1 ) ∈ R n such that ζ 0 ∈ R, (ζ 1 , . . . , ζ n−1 ) ∈ R n−1 + and n−1 j=0 ζ j ≤ max i y i . Let ∇g(ζ) 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

(max i y i − min i y i ) 2 t + 2 , (3)

where g * is the optimal solution of (2).

Proof. It it is know [18] that for all t ≥ 2:

g(ζ t ) − g * ≤ 2L(Diam(S)) 2 t + 2 ,

where L is the Lipschitz constant and and Diam(S) is the diameter of S. Let

∇ 2 g(ζ) := ∂ 2 g ∂ζ 2 0 , ∂ 2 g ∂ζ 2 1 , . . . , ∂ 2 g ∂ζ 2 n−1 .

It is well-known that if ∇g is differentiable then its Lipschitz constant L satisfies the inequality

L ≤ sup ζ ∇ 2 g(ζ) 2 . Then L ≤ sup ζ n−1 k=0 ∂ 2 g ∂ζ 2 k 2 = = 1 n n k=1 (2(n − k + 1)) 2 = 2 n n k=1 k 2 = 2 n(n + 1)(2n + 1) 6n 2 . (4)

It is easy to prove and Diam(S)

:= √ 2(max i y i − min i y i ).

The disadvantage of this method is the dependence of the theoretical degree of convergence on the dimensionality of the problem. The papers [28], [20], [29] suggest to use the values of duality gap as the stopping criterion for Frank-Wolfe type algorithms.

Empirical Result

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 = {(x i , y i ), y i = ln(x 0 + i x) + ϕ i , ϕ i ∼ N (0, 1), x 0 = 1, x = 1, i = 1, . . . , 10000. The dimension of the problem is 10000 points. The tables contain information on errors

1 n n i=1 (z i −y i ) 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 Table 1. Comparison of algorithms PAVA and greedy algorithm (Greedy) on an example of the simulated data (implementation in language C++): 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.

A = {(xi, yi), yi = ln(x0 + i x) + ϕi, ϕi ∼ N (0, 1), x0 = 1, x = 1, i = 1, . . . ,10000

Conclusion

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.