Solution de type PAVA pour un simple QP

Nov 09 2020

Laisser $l,u\in\mathbb{R}^n$, et considérez le QP:

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

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

C'est-à-dire que nous voulons minimiser le changement au carré des éléments de $x$ sujet à $x$ être au-dessus de la limite inférieure $l$ et en dessous de la limite supérieure $u$.

Mon intuition est que c'est assez simple pour qu'il devrait avoir une solution exacte de style "algorithme de violation adjacente groupée (PAVA)". Est-ce correct? Ce problème a-t-il été étudié dans la littérature antérieure?


A part: Bien sûr, comme indiqué ici, le problème peut avoir plusieurs solutions. Je me fiche de ce qui est retourné.


Plus loin: voici le code de type MATLAB pour une procédure de solution inefficace. Je pense qu'il y a un algorithme beaucoup plus efficace!

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

Réponses

5 GabrielGouvine Nov 26 2020 at 01:34

J'ai trouvé le problème intéressant, j'ai donc conçu un algorithme pour le résoudre (pas de style PAVA). Vous pouvez trouver une implémentation Python ici .

Une formulation géométrique

Les conditions d'optimalité sont, pour $1 \lt i \lt n$:

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

En utilisant ceux-ci, nous pouvons le remplacer par un problème géométrique avec les mêmes solutions optimales. Pour chaque$i$, il y a une ligne verticale avec $x=i$ et $l_i \leq y \leq u_i$. Le but est de trouver le chemin le plus court passant par tous ces segments de gauche à droite. Tout coude de ce chemin correspondra à$x_i = l_i$ ou $x_i = u_i$.

Un exemple:

Et sa solution:

L'algorithme

L'idée de l'algorithme est de trouver où sera le prochain virage et de recommencer à partir de là. Connaissant le virage précédent, nous maintenons la pente minimale et maximale pour la ligne suivante, et qui$i$ est à la limite, jusqu'à ce que nous soyons obligés d'introduire un nouveau virage.

La boucle principale est donc:

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

Il y a des majuscules pour le premier et le dernier virage, que je ne décrirai pas ici: elles sont commentées dans le code complet . La preuve de l'algorithme est laissée comme exercice au lecteur :)

Résultats de calcul

L'algorithme s'adapte parfaitement aux données aléatoires et trouve la solution optimale pour $n > 10^6$en secondes. Théoriquement, c'est seulement$O(n^2)$pour autant que je sache, mais a une complexité linéaire dans la pratique. Je suis convaincu qu'il est possible de trouver quelque chose qui$O(n \log n)$, car il présente une certaine similitude avec les algorithmes de coque convexe 2D, mais c'est amusant pour un autre jour;)