5
$\begingroup$

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
$\endgroup$
3
  • $\begingroup$ 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}]$. $\endgroup$ Nov 24 '20 at 14:48
  • $\begingroup$ Isn't that identical to what I wrote? My $\Delta x$ is a column vector. $\endgroup$
    – cfp
    Nov 24 '20 at 15:06
  • $\begingroup$ You are right never mind. $\endgroup$ Nov 24 '20 at 15:07
5
+50
$\begingroup$

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:

Geometric version of the problem

And its solution:

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

First exampleSecond example Third exampleFourth example

$\endgroup$
9
  • 1
    $\begingroup$ Could you explain how OPs formulation can be translated to your geometrical problem ? Looks interesting! $\endgroup$
    – Kuifje
    Nov 25 '20 at 19:38
  • 1
    $\begingroup$ 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. $\endgroup$ Nov 25 '20 at 20:59
  • 2
    $\begingroup$ 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. $\endgroup$ Nov 26 '20 at 11:33
  • 2
    $\begingroup$ 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! $\endgroup$ Nov 28 '20 at 9:29
  • 1
    $\begingroup$ Both are equivalent here. <= would report one more bend while it's just a straight line, but would be correct. $\endgroup$ Dec 2 '20 at 15:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.