LCP с разреженной матрицей

Я указываю матрицы заглавными буквами, а векторы - маленькими буквами.

Мне нужно решить следующую линейную систему уравнений для v.

rv >= u + Av
v >= s

Определяя z = v-s, B=rI - A, q=-u + BS, я могу переписать предыдущую проблему как проблему линейной взаимодополняемости и надеемся использовать решатель LCP, например, от openopt:

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

или, в матричной нотации:

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

Проблема в том, что моя система уравнений огромна. Чтобы создать A, I

  • Создайте четыре матрицы A11, A12, A21, A22 с помощью scipy.sparse.diags
  • И соедините их вместе как A = scipy.sparse.bmat([[A11, A12], [A21, A22]])
  • (Это также означает, что A не является симметричным, и, следовательно, некоторые эффективные переводы в QP проблемы не будут работать)

openopt.LCP, по-видимому, не может иметь дело с разреженными матрицами: когда я побежал, мой компьютер разбился. Как правило, A.todense() приведет к ошибке памяти. Аналогично, compecon-python не может решить проблемы LCP с разреженными матрицами.

Какие альтернативные реализации LCP подходят для этой проблемы?

Я действительно не думал, что необходимы образцы данных для общего "того, какие инструменты для решения LCP" были заданы, но в любом случае здесь мы идем

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>

Еще несколько указателей:

  • openopt.LCP падает с моими матрицами, я предполагаю, что он преобразует матрицы в плотные, прежде чем продолжить
  • compecon-python ошибка сбоя с некоторой ошибкой, она, по-видимому, требует плотных матриц и не имеет резерва для разреженности
  • B не является положительным полуопределенным, поэтому я не могу перефразировать проблему линейной комплементарности (LCP) как выпуклую квадратичную задачу (QP)
  • Все разрешенные QP-решения из это изложение требуют, чтобы проблема была выпуклой, а моя не
  • В Julia PATHSolver может решить мою проблему (с учетом лицензии). Однако есть проблемы, вызывающие его из Python с PyJulia (моя проблема отчет здесь)
  • Также у Matlab есть решатель LCP, который, по-видимому, может обрабатывать разреженные матрицы, но реализация еще более дурацкая (и я действительно не хочу отступать на Matlab для этого)

Ответы

Ответ 1

Используйте SUMT by Stephen Boyd. Я использовал его раньше для линейных систем, у которых есть <= или >= ограничения, и он работает очень хорошо.

SUMT = Последовательная техническая техника с минимальной минимизацией.