LCP avec matrice éparse

Question

J'indique les matrices par des lettres majuscules, et les vecteurs par des lettres minuscules.

Je dois résoudre le système d'inégalités linéaires suivant pour le vecteur v :

min(rv - (u + Av), v - s) = 0

donde 0 est un vecteur de zéros.

donde r est un scalaire, u y s sont des vecteurs, et A est une matrice.

Définition de z = v-s , B=rI - A , q=-u + Bs Je peux réécrire le problème précédent comme suit . problème de complémentarité linéaire et espérer utiliser un solveur LCP, par exemple à partir de openopt :

LCP(M, z): min(Bz+q, z) = 0

ou, en notation matricielle :

z'(Bz+q) = 0
z >= 0
Bz + q >= 0

---

Le problème est que mon système d'équations est énorme. Pour créer A , I

• Créer quatre matrices A11 , A12 , A21 , A22 en utilisant scipy.sparse.diags
• Et les empiler ensemble comme A = scipy.sparse.bmat([[A11, A12], [A21, A22]])
• (Cela signifie également que A n'est pas symétrique, et donc certaines traductions efficaces en QP les problèmes ne fonctionnent pas)

openopt.LCP ne peut apparemment pas traiter les matrices éparses : Lorsque je l'ai exécuté, mon ordinateur a planté. En général, A.todense() entraînera une erreur de mémoire. De même, compecon-python n'est pas en mesure de résoudre LCP problèmes avec des matrices éparses.

Quelle alternative LCP sont adaptées à ce problème ?

---

Je ne pensais vraiment pas données de l'échantillon était nécessaire pour une question générale "quels outils pour résoudre les LCP" était nécessaire, mais quoi qu'il en soit, nous y voici

from numpy.random import rand
from scipy import sparse

n = 3000
r = 0.03

A = sparse.diags([-rand(n)], [0])
s = rand(n,).reshape((-1, 1))
u = rand(n,).reshape((-1, 1))

B = sparse.eye(n)*r - A
q = -u + B.dot(s)

q.shape
Out[37]: (3000, 1)
B
Out[38]: 
<3000x3000 sparse matrix of type '<class 'numpy.float64'>'
    with 3000 stored elements in Compressed Sparse Row format>

---

Quelques conseils supplémentaires :

• openopt.LCP se plante avec mes matrices, je suppose qu'il convertit les matrices en dense avant de continuer
• compecon-python échoue carrément avec une erreur, il requiert apparemment des matrices denses et n'a pas de solution de rechange pour la sparsité.
• B n'est pas semi-définie positive, donc je ne peux pas reformuler le problème de complémentarité linéaire (LCP) comme un problème quadratique convexe (QP)
• Tous les solveurs de QP sparse de cette exposition exigent que le problème soit convexe, ce que le mien n'est pas.
• A Julia, PATHSolver peut résoudre mon problème (avec une licence). Cependant, il y a des problèmes pour l'appeler depuis Python avec PyJulia ( mon rapport sur les problèmes ici )
• Matlab dispose également d'un solveur LCP qui peut apparemment gérer les matrices éparses, mais son implémentation est encore plus farfelue (et je ne veux vraiment pas me rabattre sur Matlab pour cela).

Answer 1

Ce problème a une solution très efficace (temps linéaire), bien qu'elle nécessite un peu de discussion...

Zeroth : clarifier le problème / LCP

D'après les clarifications apportées dans les commentaires, @FooBar dit que le problème original est lié aux éléments. min ; nous devons trouver un z (o v ) tel que

soit l'argument de gauche est nul et l'argument de droite est non négatif, soit l'argument de gauche est non négatif et l'argument de droite est nul.

Comme le souligne correctement @FooBar, une reparamétrisation valide conduit à un LCP. Cependant, je montre ci-dessous qu'une solution plus facile et plus efficace au problème original peut être obtenue (en exploitant la structure du problème original) sans avoir besoin de la LCP. Pourquoi cela serait-il plus facile ? Eh bien, remarquez qu'une LCP a un terme quadratique dans z (Bz+q)'z, mais que le problème original ne le fait pas (seulement des termes linéaires Bz+q et z). Je vais exploiter ce fait ci-dessous.

Premièrement : simplifier

Il existe un détail important mais essentiel qui simplifie ce problème de façon majeure :

• Créer quatre matrices A11, A12, A21, A22 en utilisant scipy.sparse.diags
• Et les empiler ensemble comme A = scipy.sparse.bmat([[A11, A12], [A21, A22]])

Cela a d'énormes implications. Plus précisément, ce n'est pas simple grand site problème, mais plutôt un grand nombre de vraiment petit (2D, pour être précis). Remarquez que la structure diagonale en bloc de ce A est préservée dans toutes les opérations ultérieures. Et chaque opération ultérieure est une multiplication matrice-vecteur ou un produit interne. Cela signifie qu'en réalité ce programme est séparable par paires z (o v ) variables.

Pour être précis, supposons que chaque bloc A11,... est la taille n par n . On remarque alors de façon critique que z_1 y z_{n+1} apparaître uniquement en équations et en termes avec eux-mêmes, et jamais ailleurs. Le problème est donc séparable en n des problèmes, chacun d'entre eux étant bidimensionnel. Ainsi, je résoudrai ci-après le problème 2D, et vous pourrez sérialiser ou paralléliser sur n comme bon vous semble, sans matrices éparses ni gros paquets d'options.

Deuxièmement : la géométrie du problème 2D

Supposons que nous ayons le problème des éléments en 2D, à savoir :

find z such that (elementwise) min( Bz + q , z ) = 0, or declare that no such `z` exists.

Parce que dans notre configuration B est maintenant 2x2, cette géométrie du problème correspond à quatre inégalités scalaires que nous pouvons énumérer manuellement (je les ai nommées a1, a2, z1, z2) :

“a1”: B11*z1 + B12*z2 + q1 >=0
“a2”: B21*z1 + B22*z2 + q2 >=0
“z1”: z1 >= 0
“z2:” z2 >= 0

Cela représente un polyèdre (éventuellement vide), c'est-à-dire l'intersection de quatre demi-espaces dans l'espace 2D.

Troisièmement : résoudre le problème de la 2D

(Edition : mise à jour de cette partie pour plus de clarté)

Quel est donc le problème spécifique de la 2D ? Nous voulons trouver un z où l'on obtient l'une des solutions suivantes (qui ne sont pas toutes distinctes, mais cela n'a pas d'importance) :

1. a1>=0, z1=0, a2>=0, z2=0
2. a1=0, z1>=0, a2=0, z2>=0
3. a1>=0, z1=0, a2=0, z2>=0
4. a1=0, z1>=0, a2>=0, z2=0

Si l'un d'entre eux est atteint, nous avons une solution : a z où le minimum par élément de z et Bz+q est le vecteur 0. Si aucun des quatre n'est atteint, le problème est infaisable et nous déclarons qu'aucune solution n'est possible. z existe.

Le premier cas se produit lorsque le point d'intersection de a1,a2 est dans l'orthant positif. Il suffit de choisir la solution z = B^-1q, et de vérifier si elle est élémentairement non négative. Si c'est le cas, on retourne B^-1q (notez que, même si B n'est pas psd, je suppose qu'il est inversible/de rang complet). Donc :

if B^-1q is elementwise nonnegative, return z = B^-1q.

Le deuxième cas (qui n'est pas entièrement distinct du premier) se produit lorsque le point d'intersection de a1,a2 est n'importe où mais contient l'origine. En d'autres termes, lorsque Bz+q >0 pour z = 0. Cela se produit lorsque q est positif par élément. Donc :

elif q is elementwise nonnegative, return z = 0.

Dans le troisième cas, z1=0, ce qui peut être substitué à a2 pour montrer que a2=0 lorsque z2 = -q2/B22. z2 doit être >=0, donc -q2/B22 >=0. a1 doit également être >=0, ce qui, en substituant les valeurs de z1 et z2, donne -B22/B12*q2 + q1 >=0. Ainsi :

elif -q2/B22 >=0 and  -B22/B12*q2 + q1 >=0, return z1= 0, z2 = -q2/B22.

La quatrième étape est la même que la troisième, mais en intervertissant 1 et 2. Donc :

elif -q1/B11 >=0 and -B21/B11*q1 + q2 >=0, return z1 = -q1/B11, z2 =0.

Le cinquième et dernier cas est celui où le problème est infaisable. Cela se produit lorsqu'aucune des conditions ci-dessus n'est remplie. Donc :

else return None 

Enfin : assembler les pièces

La résolution de chaque problème 2D consiste en un couple de solutions linéaires et de comparaisons simples/effectives/triviales. Chacune renvoie une paire de nombres ou None . Ensuite, il suffit de faire la même chose pour tout le monde. n problèmes 2D, et concaténer le vecteur. Si l'un d'entre eux est None, le problème est infaisable (tous None). Sinon, vous avez votre réponse.

Answer 2

Solveur LCP pour Python basé sur AMPLPY

comme @denfromufa l'a souligné, il y a une AMPL à l'interface PATH résolveur. Il a suggéré Pyomo puisqu'il est open source et capable de traiter AMPL . Cependant, Pyomo s'est avéré être lent et fastidieux à travailler. J'ai finalement écrit ma propre interface avec le PATH J'ai un solveur en cython et j'espère le publier à un moment donné, mais pour le moment, il n'a pas d'assainissement d'entrée, il est rapide et sale et je ne vois pas le temps de travailler dessus.

Pour l'instant, je peux partager une réponse qui utilise l'extension python de AMPL . Ce n'est pas aussi rapide qu'une interface directe avec PATH : Pour chaque LCP que nous voulons résoudre, il crée un fichier modèle (temporaire), exécute AMPL et recueille le résultat. C'est un peu rapide et sale, mais j'ai pensé que je devais au moins rapporter une partie des résultats obtenus au cours des derniers mois depuis que j'ai posé cette question.

import os
# PATH license string needs to be updated
os.environ['PATH_LICENSE_STRING'] = "3413119131&Courtesy&&&USR&54784&12_1_2016&1000&PATH&GEN&31_12_2017&0_0_0&5000&0_0"

from amplpy import AMPL, Environment, dataframe
import numpy as np
from scipy import sparse
from tempfile import mkstemp
import os

import sys
import contextlib

class DummyFile(object):
    def write(self, x): pass

@contextlib.contextmanager
def nostdout():
    save_stdout = sys.stdout
    sys.stdout = DummyFile()
    yield
    sys.stdout = save_stdout

class modFile:
    # context manager: create temporary file, insert model data, and supply file name
    # apparently, amplpy coders are inable to support StringIO

    content = """
        set Rn;

        param B {Rn,Rn} default 0;
        param q {Rn} default 0;

        var x {j in Rn};

        s.t. f {i in Rn}:
                0 <= x[i]
             complements
                sum {j in Rn} B[i,j]*x[j]
                 >= -q[i];
    """

    def __init__(self):
        self.fd = None
        self.temp_path = None

    def __enter__(self):
        fd, temp_path = mkstemp()
        file = open(temp_path, 'r+')
        file.write(self.content)
        file.close()

        self.fd = fd
        self.temp_path = temp_path
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        os.close(self.fd)
        os.remove(self.temp_path)

def solveLCP(B, q, x=None, env=None, binaryDirectory=None, pathOptions={'logfile':'logpath.tmp' }, verbose=False):
    # x: initial guess
    if binaryDirectory is not None:
        env = Environment(binaryDirectory='/home/foo/amplide.linux64/')
    if verbose:
        pathOptions['output'] = 'yes'
    ampl = AMPL(environment=env)

    # read model
    with modFile() as mod:
        ampl.read(mod.temp_path)

    n = len(q)
    # prepare dataframes
    dfQ = dataframe.DataFrame('Rn', 'c')
    for i in np.arange(n):
        # shitty amplpy does not support np.float
        dfQ.addRow(int(i)+1, np.float(q[i]))

    dfB = dataframe.DataFrame(('RnRow', 'RnCol'), 'val')

    if sparse.issparse(B):
        if not isinstance(B, sparse.lil_matrix):
            B = B.tolil()
        dfB.setValues({
            (i+1, j+1): B.data[i][jPos]
            for i, row in enumerate(B.rows)
            for jPos, j in enumerate(row)
            })

    else:
        r = np.arange(n) + 1
        Rrow, Rcol = np.meshgrid(r, r, indexing='ij')
        #dfB.setColumn('RnRow', [np.float(x) for x in Rrow.reshape((-1), order='C')])
        dfB.setColumn('RnRow', list(Rrow.reshape((-1), order='C').astype(float)))
        dfB.setColumn('RnCol', list(Rcol.reshape((-1), order='C').astype(float)))
        dfB.setColumn('val', list(B.reshape((-1), order='C').astype(float)))

    # set values
    ampl.getSet('Rn').setValues([int(x) for x in np.arange(n, dtype=int)+1])
    if x is not None:
        dfX = dataframe.DataFrame('Rn', 'x')
        for i in np.arange(n):
            # shitty amplpy does not support np.float
            dfX.addRow(int(i)+1, np.float(x[i]))
        ampl.getVariable('x').setValues(dfX)

    ampl.getParameter('q').setValues(dfQ)
    ampl.getParameter('B').setValues(dfB)

    # solve
    ampl.setOption('solver', 'pathampl')

    pathOptions = ['{}={}'.format(key, val) for key, val in pathOptions.items()]
    ampl.setOption('path_options', ' '.join(pathOptions))
    if verbose:
        ampl.solve()
    else:
        with nostdout():
            ampl.solve()

    if False:
        bD = ampl.getParameter('B').getValues().toDict()
        qD = ampl.getParameter('q').getValues().toDict()
        xD = ampl.getVariable('x').getValues().toDict()
        BB = ampl.getParameter('B').getValues().toPandas().values.reshape((n, n,), order='C')
        qq = ampl.getParameter('q').getValues().toPandas().values[:, 0]
        xx = ampl.getVariable('x').getValues().toPandas().values[:, 0]
        ineq2 = BB.dot(xx) + qq
        print((xx * ineq2).min(), (xx * ineq2).max() )
    return ampl.getVariable('x').getValues().toPandas().values[:, 0]

if __name__ == '__main__':

    # solve problem from the Julia port at https://github.com/chkwon/PATHSolver.jl
    n = 4
    B = np.array([[0, 0, -1, -1], [0, 0, 1, -2], [1, -1, 2, -2], [1, 2, -2, 4]])
    q = np.array([2, 2, -2, -6])

    BSparse = sparse.lil_matrix(B)

    env = Environment(binaryDirectory='/home/foo/amplide.linux64/')
    print(solveLCP(B, q, env=env))
    print(solveLCP(BSparse, q, env=env))

    # to test licensing
    from numpy import random
    n = 1000
    B = np.diag((random.randn(n)))
    q = np.ones((n,))
    print(solveLCP(B, q, env=env))