Return to Answer

The idea of the algorithm is to find where the next bend is going to be, and restart from there. Knowing the last 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.

def find_next_bend(n, l, u, x, last_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[last_bend]) / (i-last_bend)
        new_max_slope = (u[i] - x[last_bend]) / (i-last_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 breakpoint, 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 :)

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.

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

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

A geometric formulation

If we focus on these optimality conditions, we can see it as a geometric problem. 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$.

A geometric formulation

• $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$.

We can see it as a geometric problem. 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.

There are a few cornercases, but the idea is to find where the next bend is going to be, and restart from there. 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 ;)

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} \lt x_{i+1} - x_i~$ if $x_i = u_i$
• $x_i - x_{i-1} \gt x_{i+1} - x_i~$ if $x_i = l_i$

A geometric formulation

If we focus on these optimality conditions, we can see it as a geometric problem. 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$.

The algorithm

The idea of the algorithm is to find where the next bend is going to be, and restart from there. Knowing the last 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, last_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[last_bend]) / (i-last_bend)
        new_max_slope = (u[i] - x[last_bend]) / (i-last_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 breakpoint, 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 ;)