PAVA-like solution to simple QP

Question

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

Answer 1

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 ;)