Решить систему линейных целочисленных уравнений в Python

Я ищу метод для решения системы линейных уравнений в Python. В частности, я ищу наименьший целочисленный вектор, который больше всех нулей и решает данное уравнение. Например, у меня есть следующее уравнение:

введите описание изображения здесь

и хочу решить введите описание изображения здесь,

В этом случае наименьший целочисленный вектор, который решает это уравнение введите описание изображения здесь,

Однако как я могу определить это решение автоматически? Если я использую scipy.optimize.nnls, лайк

A = np.array([[1,-1,0],[0,2,-1],[2,0,-1]])
b = np.array([0,0,0])
nnls(A,b)

результат (array([ 0., 0., 0.]), 0.0), Что тоже правильно, но не желаемое решение...

Редактировать: я прошу прощения за неточность в определенных аспектах. Если кого-то интересуют подробности, проблема возникает из статьи "Статическое планирование программ синхронного потока данных для цифровой обработки сигналов", Эдвард А. Ли и Дэвид Г. Мессершмитт, IEEE Transactions on Computers, Vol. С-36, № 1, с. 24-35, январь 1987 г.

Теорема 2 говорит

Для связного графа SDF с s узлами и матрицей топологии A и с рангом (A)=s-2 мы можем найти положительный целочисленный вектор b!= 0 такой, что Ab = 0, где 0 - нулевой вектор.

Сразу после доказательства теоремы 2 говорят

Может быть желательно найти для наименьшего положительного целочисленного вектора в пустом пространстве. Чтобы сделать это, уменьшите каждую рациональную запись в u'так, чтобы ее числитель и знаменатель были относительно простыми. Алгоритм Евклида будет работать для этого.

Источник

6 ответов

Решение

Чтобы найти точное решение, которое вы хотите, numpy и scipy, вероятно, не самые лучшие инструменты. Их алгоритмы обычно работают с плавающей запятой и не гарантируют точного ответа.

Ты можешь использовать sympy чтобы получить точный ответ на эту проблему. Matrix класс в sympy обеспечивает метод nullspace() который возвращает список базовых векторов для пустого пространства. Вот пример:

In [20]: from sympy import Matrix, lcm

In [21]: A = Matrix([[1, -1, 0], [0, 2, -1], [2, 0, -1]])

Получить вектор в пустом пространстве. nullspace() возвращает список; этот код предполагает, что ранг A равен 2, поэтому список имеет длину один:

In [22]: v = A.nullspace()[0]

In [23]: v
Out[23]: 
Matrix([
[1/2],
[1/2],
[  1]])

Найти наименьшее общее кратное знаменателей значений в v так что мы можем масштабировать результат до наименьших целых чисел:

In [24]: m = lcm([val.q for val in v])

In [25]: m
Out[25]: 2

x содержит целочисленный ответ:

In [26]: x = m*v

In [27]: x
Out[27]: 
Matrix([
[1],
[1],
[2]])

Чтобы преобразовать этот результат в массив целых чисел, вы можете сделать что-то вроде:

In [52]: np.array([int(val) for val in x])
Out[52]: array([1, 1, 2])
Источник

На самом деле, это просто базовая линейная алгебра.

>>> A = np.array([[1,-1,0], [0,2,-1],[2,0,-1]])

Давайте вычислим собственные значения и собственный вектор для этой матрицы.

>>> e = np.linalg.eig(A)
>>> e
(array([ -4.14213562e-01,  -1.05471187e-15,   2.41421356e+00]), array([[-0.26120387, -0.40824829,  0.54691816],
       [-0.36939806, -0.40824829, -0.77345908],
       [-0.89180581, -0.81649658,  0.32037724]]))    
>>>> np.round(e[0], 10)
array([-0.41421356, -0.        ,  2.41421356])

После округления мы видим, что 0 является собственным значением нашей матрицы A. Таким образом, собственный вектор s для 0-собственного значения является хорошим кандидатом для вашей системы уравнений.

>>> s = e[1][:,1]
>>> s
array([-0.40824829, -0.40824829, -0.81649658])

Умножение собственного вектора на соответствующую матрицу приводит к получению самого собственного вектора, масштабированного на соответствующее собственное значение. Итак, в приведенном выше случае мы видим: As = 0s = 0

>>> np.round(A.dot(s), 10)
array([ 0.,  0.,  0.])

Поскольку нас интересует целочисленное решение, мы должны масштабировать вектор решения.

>>> x = s / s[1]
>>> x
array([ 1.,  1.,  2.])

Надеюсь, что этот ответ решит вашу проблему.

Источник

В Pypi есть решение для такого рода уравнения. Очевидно, он может вычислить нормальную форму матрицы Эрмита, которая, в свою очередь, может быть использована для решения вашей проблемы.

Пакет имеет версию 0.1.

Мудрец также поддерживает нормальную форму Эрмита.

Частный случай однородной системы, т. Е. B =0, немного проще, вот решатель для простейшего возможного случая матрицы, имеющей ранг n-1

import sympy
import numpy as np

def create_rd(n, defect=1, range=10):
    while True:
        res = sympy.Matrix((np.random.randint(-range+1,range,(n, n-defect))
                           @np.random.randint(0,2,(n-defect, n)))
                           .astype(object))
        if res.rank() == n-defect:
            break
    return res

def solve(M):
    ns = M.nullspace()
    ns = [n / sympy.gcd(list(n)) for n in ns]
    nsnp = np.array([[int(k) for k in n] for n in ns])
    if len(ns) == 1:
        return ns[0], nsnp[0]
    else:
        raise NotImplementedError

Образец вывода:

>>> M = create_rd(4)  # creates a rank-deficient matirx
>>> ns, nn = solve(M) # finds the 1d nullspace and a minimal integer basis vector
>>> M
Matrix([
[-7, -7, -7, -12],
[ 0,  6,  0,   6],
[ 4,  1,  4,  -3],
[-4, -7, -4,  -9]])
>>> ns
Matrix([
[-1],
[ 0],
[ 1],
[ 0]])
>>> M*ns
Matrix([
[0],
[0],
[0],
[0]])
>>> M = create_rd(40) # we can go to higher dimensions
>>> ns, nn = solve(M) # but solutions quickly become unwieldy
>>> ns
Matrix([
[  8620150337],
[-48574455644],
[ -6216916999],
[-14703127270],
 < - snip - >
Источник

Вот решение с использованием ortools:

      s = Solver('')
A = [[1, -1, 0], [0, 2, -1], [2, 0, -1]]
vars = [s.IntVar(1, 2**32) for _ in A[0]]
for row in A:
    s.Add(s.ScalProd(vars, row) == 0)
s.NewSearch(s.Phase(vars, s.CHOOSE_FIRST_UNBOUND, s.ASSIGN_MIN_VALUE))
s.NextSolution()
print([var.Value() for var in vars])
Источник

Этот вопрос довольно неформальный, как видно из комментариев. Не зная вашего определения наименьшего (примеры: l1-норма, l2-норма), трудно ответить на вашу конкретную проблему.

Общая проблема связана с решением системы диофантовых уравнений, но их трудно решить (в общем случае), и программного обеспечения не так много.

Естественным подходом является использование Integer-программирования, которое является NP-сложным, но коммерческое и некоторые решатели с открытым исходным кодом очень и очень мощные.

В numpy/scipy нет встроенного метода, который бы решал вашу проблему без огромных модификаций (например: реализация некоторого собственного алгоритма, основанного на numpy/scipy). К сожалению, в numpy/scipy нет IP-решателя.

Давайте предположим:

  • наименьшая - это сумма переменных (в большинстве случаев это неправильный подход; l2-норма более популярна, но ее сложнее сформулировать и она требует более мощных решателей)
  • Вы хотите запретить нулевой вектор
  • x is nonnegative

Вот несколько простых реализаций на основе IP с использованием целлюлозы и numpy. Я не очень люблю целлюлозу, но ее легко установить (pip install pulp) на всех популярных системах.

Код

from pulp import *
import numpy as np

EPS = 1e-3

""" Input """
A = np.array([[1,-1,0],[0,2,-1],[2,0,-1]])
b = np.array([0,0,0])

""" MIP """
# Variables
x = np.empty(b.shape[0], dtype=object)
for i in range(b.shape[0]):
    x[i] = LpVariable("x" + str(i), lowBound=0, upBound=None, cat='Integer')

# Problem
prob = LpProblem("prob", LpMinimize)

# Objective
prob += np.sum(x)

# Constraints
for row in range(A.shape[0]):
    prob += np.dot(A[row], x) == b[row]

prob += np.sum(x) >= EPS  # forbid zero-vector

# Solve
status = prob.solve()
print(LpStatus[status])
print([value(x_) for x_ in x])

Выход

Optimal
[1.0, 1.0, 2.0]
Источник

Возможно, это дикая идея, но звучит так, будто вы пытаетесь создать решение для ограничения.

Minizinc - это универсальный решатель ограничений. Возможно, можно выразить ваше ограничение таким образом, чтобы миницинк мог его решить?

Затем кажется, что есть библиотека Python для взаимодействия с ней: https://pypi.python.org/pypi/pymzn/

Источник
Другие вопросы по тегам