# PAVA-like solution to simple QP

Let $$l,u\in\mathbb{R}^n$$, and consider the QP:

$$\min_{l\le x\le u} {(\Delta x)^\top (\Delta x)}$$

where $$\Delta x=[x_2-x_1,\,x_3-x_2,\,\dots,\,x_n-x_{n-1}]^\top$$.

I.e. we want to minimize the squared change in the elements of $$x$$ subject to $$x$$ being above the lower bound $$l$$ and below the upper bound $$u$$.

My hunch is that this is simple enough that it ought to have an exact "pooled adjacent violators algorithm (PAVA)" style solution. Is this correct? Has this problem been studied in the prior literature?

Aside: Of course as stated here the problem may have multiple solutions. I do not care which is returned.

Further aside: Here's MATLAB type code for an inefficient solution procedure. I expect there's a much more efficient algorithm!

while true

xo = x;

x( 1 ) = max( l( 1 ), min( u( 1 ), x( 2 ) ) );
for i = 2 : ( n - 1 )
x( i ) = max( l( i ), min( u( i ), 0.5 * ( x( i - 1 ) + x( i + 1 ) ) ) );
end
x( n ) = max( l( n ), min( u( n ), x( n - 1 ) ) );

if all( abs( x - xo ) < 1e-12 )
break
end

end

• I think it should be: $$\min_{l\le x\le u} {(\Delta x) (\Delta x)^\top}$$ where $\Delta x=[x_2-x_1,\,x_3-x_2,\,\dots,\,x_n-x_{n-1}]$. Nov 24 '20 at 14:48
• Isn't that identical to what I wrote? My $\Delta x$ is a column vector.
– cfp
Nov 24 '20 at 15:06
• You are right never mind. Nov 24 '20 at 15:07

I found the problem interesting, so I designed an algorithm to solve it (not PAVA-style). You can find a Python implementation here.

## A geometric formulation

The optimality conditions are, for $$1 \lt i \lt n$$:

• $$x_i - x_{i-1} = x_{i+1} - x_i~$$ if $$l_i \lt x_i \lt u_i$$
• $$x_i - x_{i-1} \leq x_{i+1} - x_i~$$ if $$x_i = u_i$$
• $$x_i - x_{i-1} \geq x_{i+1} - x_i~$$ if $$x_i = l_i$$

Using those, we can replace it by a geometric problem with the same optimal solutions. For each $$i$$, there is a vertical line with $$x=i$$ and $$l_i \leq y \leq u_i$$. The goal is to find the shortest path going through all these segments from left to right. Any bend in this path will correspond to $$x_i = l_i$$ or $$x_i = u_i$$.

An example: And its solution: ## The algorithm

The idea of the algorithm is to find where the next bend is going to be, and restart from there. Knowing the previous bend, we maintain the minimum and maximum slope for the next line, and which $$i$$ is at the limit, until we are forced to introduce a new bend.

So the core loop is:

def find_next_bend(n, l, u, x, previous_bend):
min_slope = -float("inf")
max_slope = float("inf")
min_slope_ind = -1
max_slope_ind = -1
for i in range(last_bend+1, n):
new_min_slope = (l[i] - x[previous_bend]) / (i-previous_bend)
new_max_slope = (u[i] - x[previous_bend]) / (i-previous_bend)
if new_min_slope > max_slope:
# Bend upwards
x[max_slope_ind] = u[max_slope_ind]
return max_slope_ind
if new_max_slope < min_slope:
# Bend downwards
x[min_slope_ind] = l[min_slope_ind]
return min_slope_ind
if new_max_slope < max_slope:
# Restrict max possible slope
max_slope_ind = i
max_slope = new_max_slope
if new_min_slope > min_slope:
# Restrict min possible slope
min_slope_ind = i
min_slope = new_min_slope


There are cornercases for the first and last bend, that I won't describe here: they are commented in the full code. The proof of the algorithm is left as an exercise for the reader :)

## Computational results

The algorithm scales perfectly on random data, and finds the optimal solution for $$n > 10^6$$ in seconds. Theoretically, it is only $$O(n^2)$$ as far as I can tell, but has linear complexity in practice. I am convinced that it's possible to find something that is $$O(n \log n)$$, as it bears some similarity to 2D convex hull algorithms, but that's fun for another day ;)

• Could you explain how OPs formulation can be translated to your geometrical problem ? Looks interesting! Nov 25 '20 at 19:38
• Thanks! I don't have much time right now, but I'll try detailing it tomorrow. I think the easiest way to see the equivalence of the formulations despite the different objective function is to look at the optimality conditions: in both cases, they are $x_i - x_{i-1} = x_{i+1} - x_i$ if $x_i$ is not at its bound, and $x_i - x_{i-1} \gt x_{i+1} - x_i$ or $x_i - x_{i-1} \lt x_{i+1} - x_i$ depending on the bound otherwise. Nov 25 '20 at 20:59
• I just added some code and the above explanation in the answer. For the full algorithm, you will still have to look at the full Python code: I don't want to detail the cornercases in the answer, since they are not interesting and make the algorithm quite long. Nov 26 '20 at 11:33
• Thank you :-) If only all QP were like this :D I'd like to see someone come up with an O(n log n) algorithm, too! Nov 28 '20 at 9:29
• Both are equivalent here. <= would report one more bend while it's just a straight line, but would be correct. Dec 2 '20 at 15:02