matriz sparse python

LCP con matriz dispersa

Question

Indico matrices con mayúsculas y vectores con minúsculas.

Necesito resolver el siguiente sistema de desigualdades lineales para el vector v :

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

donde 0 es un vector de ceros.

donde r es un escalar, u y s son vectores, y A es una matriz.

Al definir z = vs , B=rI - A , q=-u + Bs , puedo reescribir el problema anterior como un problema de complementariedad lineal y espero usar un solucionador LCP, por ejemplo de openopt :

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

o, en notación matricial:

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

El problema es que mi sistema de ecuaciones es enorme. Para crear A , yo

• Cree cuatro matrices A11 , A12 , A21 , A22 usando scipy.sparse.diags
• Y A = scipy.sparse.bmat([[A11, A12], [A21, A22]]) juntos como A = scipy.sparse.bmat([[A11, A12], [A21, A22]])
• (Esto también significa que A no es simétrico y, por lo tanto, algunas traducciones eficientes a problemas de QP no funcionarán)

openopt.LCP aparentemente no puede manejar matrices dispersas: cuando ejecuté esto, mi computadora se bloqueó. En general, A.todense() provocará un error de memoria. Del mismo modo, compecon-python no puede resolver problemas de LCP con matrices dispersas.

¿Qué implementaciones alternativas de LCP son adecuadas para este problema?

Realmente no pensé que se necesitaran datos de muestra para una pregunta general sobre "qué herramientas para resolver el LCP", pero de todos modos, aquí vamos

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>

Algunos consejos más:

• openopt.LCP bloquea con mis matrices, supongo que convierte las matrices en densas antes de continuar
• compecon-python falla por completo con algún error, aparentemente requiere matrices densas y no tiene respaldo para la escasez
• B no es positivo semi-definido, por lo que no puedo reformular el problema de complementariedad lineal (LCP) como un problema cuadrático convexo (QP)
• Todos los solucionadores dispersos de QP de esta exposición requieren que el problema sea convexo, que el mío no es
• En Julia, PATHSolver puede resolver mi problema (con una licencia). Sin embargo, hay problemas para llamarlo desde Python con PATHSolver ( mi informe de problema aquí )
• También Matlab tiene un solucionador LCP que aparentemente puede manejar matrices dispersas, pero su implementación es aún más extraña (y realmente no quiero recurrir a Matlab para esto)

Answer 1

Solucionador LCP para Python basado en AMPLPY

Como señaló @denfromufa, hay una interfaz AMPL para el solucionador PATH . Sugirió Pyomo ya que es de código abierto y capaz de procesar AMPL . Sin embargo, Pyomo resultó ser lento y tedioso para trabajar. Finalmente, he escrito mi propia interfaz para el solucionador PATH en cython y espero lanzar eso en algún momento, pero en este momento no tiene entrada de saneamiento, es rápido y sucio y no veo el momento de trabajar en ello.

Por ahora, puedo compartir una respuesta que usa la extensión python de AMPL . No es tan rápido como una interfaz directa a PATH : para cada LCP que queremos resolver, crea un archivo de modelo (temporal), ejecuta AMPL y recopila la salida. Es algo rápido y sucio, pero sentí que al menos debería informar partes de los resultados de los últimos meses desde que hice esta pregunta.

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

Answer 2

Este problema tiene una solución muy eficiente (tiempo lineal), aunque requiere un poco de discusión ...

Zeroth: aclarando el problema / LCP

Según aclaraciones en los comentarios, @FooBar dice que el problema original es elementwise min ; necesitamos encontrar una z (o v ) tal que

o el argumento izquierdo es cero y el argumento derecho no es negativo o el argumento izquierdo no es negativo y el argumento derecho es cero

Como @FooBar señala correctamente, una reparametrización válida conduce a un LCP. Sin embargo, a continuación muestro que se puede lograr una solución más fácil y más eficiente para el problema original (explotando la estructura en el problema original) sin necesitar el LCP. ¿Por qué debería ser esto más fácil? Bueno, observe que un LCP tiene un término cuadrático en z (Bz + q) ''z, pero que el problema original no lo tiene (solo términos lineales Bz + q y z). Explotaré ese hecho a continuación.

Primero: simplificar

Hay un detalle importante pero clave que simplifica este problema de manera importante:

• Cree cuatro matrices A11, A12, A21, A22 usando scipy.sparse.diags
• Y apilarlos juntos como A = scipy.sparse.bmat ([[A11, A12], [A21, A22]])

Esto tiene enormes implicaciones. Específicamente, este no es un gran problema, sino una gran cantidad de problemas realmente pequeños (2D, para ser precisos). Observe que la estructura diagonal de bloque de esta matriz A se conserva en todas las operaciones posteriores. Y cada operación posterior es una multiplicación matriz-vector o un producto interno. Esto significa que realmente este programa es separable en pares de variables z (o v ).

Para ser específicos, supongamos que cada bloque A11,... tiene un tamaño n por n . Luego observe críticamente que z_1 y z_{n+1} aparecen solo en ecuaciones y términos con ellos mismos, y nunca en otra parte. Entonces el problema es separable en n problemas, cada uno de los cuales es bidimensional. Por lo tanto, en lo sucesivo resolveré el problema 2D, y puede serializar o paralelizar sobre n como mejor le parezca, sin matrices dispersas o paquetes de gran opción.

Segundo: la geometría del problema 2D

Supongamos que tenemos el problema de elementos en 2D, a saber:

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

Debido a que en nuestra configuración B ahora es 2x2, la geometría de este problema corresponde a cuatro desigualdades escalares que podemos enumerar manualmente (las he llamado a1, a2, z1, z2):

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

Esto representa un poliedro (posiblemente vacío), también conocido como la intersección de cuatro medios espacios en el espacio 2d.

Tercero: resolver el problema 2D

(Editar: actualizado este bit para mayor claridad)

¿Cuál es específicamente el problema 2D entonces? Queremos encontrar una z donde se logre una de las siguientes soluciones (que no son todas distintas, pero eso no importará):

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 se logra uno de estos, tenemos una solución: az, donde el elemento mínimo mínimo de z y Bz + q es el vector 0. Si no se puede lograr ninguno de los cuatro, el problema es inviable y declararemos que no existe tal z .

El primer caso ocurre cuando el punto de intersección de a1, a2 está en el lado positivo. Simplemente elija la solución z = B ^ -1q, y verifique si es no negativa por elementos. Si es así, devuelva B^-1q (tenga en cuenta que, aunque B no sea psd, supongo que es invertible / rango completo). Asi que:

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

El segundo caso (no completamente distinto del primero) ocurre cuando el punto de intersección de a1, a2 está en cualquier lugar pero contiene el origen. En otras palabras, siempre que Bz + q> 0 para z = 0. Esto ocurre cuando q es positivo por elementos. Asi que:

elif q is elementwise nonnegative, return z = 0.

El tercer caso tiene z1 = 0, que puede sustituirse en a2 para mostrar que a2 = 0 cuando z2 = -q2 / B22. z2 debe ser> = 0, entonces -q2 / B22> = 0. a1 también debe ser> = 0, que sustituyendo en los valores de z1 y z2, se obtiene -B22 / B12 * q2 + q1> = 0. Asi que:

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

El cuarto paso es el mismo que el tercero, pero intercambiando 1 y 2. Entonces:

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

El quinto caso final es cuando el problema no es factible. Eso ocurre cuando no se cumple ninguna de las condiciones anteriores. Asi que:

else return None 

Finalmente: junta las piezas

Resolver cada problema 2D es un par de soluciones lineales simples / eficientes / triviales y se compara. Cada uno devolverá un par de números o None . Luego haga lo mismo en todos los problemas n 2D y concatene el vector. Si alguno es Ninguno, el problema es inviable (todos Ninguno). De lo contrario, tienes tu respuesta.