решить прямоугольную матрицу в Python, чтобы получить решение с произвольными параметрами - PullRequest
Новая
       27

решить прямоугольную матрицу в Python, чтобы получить решение с произвольными параметрами

3 голосов
/ 16 августа 2011

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

print matrix_a

print vector_b

print len(matrix_a),len(matrix_a[0])

дает: 

 [[0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0],
 [0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],
 [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0],
 [0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
 [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0],
 [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0],
 [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]]

[2, 0, 2, 1, 2, 1, 1, 1, 1, 1, 1]

11 26

мой полный код на http://codepad.org/tSgstpYe

Как видите, у меня есть система Ax = b. Я знаю, что каждое значение x x1, x2 .. должно быть либо 1, либо 0, и я ожидаю, что с этим ограничением должно быть только одно решение для системы.

На самом деле я ожидаю, что ответ x = [0,1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0 , 0,0,1,0,0,0,0]

Я посмотрел на http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html#numpy.linalg.solve но, похоже, он принимает только квадратные матрицы.

Любая помощь в решении этой системы была бы великолепна!

Ответы [ 3 ]

3 голосов
/ 17 августа 2011

Вот простая реализация (с жестко заданными пороговыми значениями), но она дает решение, которое вы ищете с тестовыми данными.

Он основан на Итеративно переоцененных наименьших квадратах . 

from numpy import abs, diag, dot, zeros
from numpy.linalg import lstsq, inv, norm

def irls(A, b, p= 1):
    """Solve least squares problem min ||x||_p s.t. Ax= b."""
    x, p, e= zeros((A.shape[1], 1)), p/ 2.- 1, 1.
    xp= lstsq(A, b)[0]
    for k in xrange(1000):
        if e< 1e-6: break
        Q= dot(diag(1./ (xp** 2+ e)** p), A.T)
        x= dot(dot(Q, inv(dot(A, Q))), b)
        x[abs(x)< 1e-1]= 0
        if norm(x- xp)< 1e-2* e** .5: e*= 1e-1
        xp= x
    return k, x.round()
2 голосов
/ 16 августа 2011

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

def testSolution(a, b, x):
  result = 0
  for i in range(len(b)):
    n = 0
    for j in range(len(a[i])):
      n += a[i][j] * x[j]
    if n < b[i]:
      result = -1
    elif n > b[i]:
      return 1
  return result

def solve(a, b):
  def solveStep(a, b, result, step):
    if step >= len(result):
      return False

    result[step] = 1
    areWeThere = testSolution(a, b, result)
    if areWeThere == 0:
      return True
    elif areWeThere < 0 and solveStep(a, b, result, step + 1):
      return True
    result[step] = 0
    return solveStep(a, b, result, step + 1)

  result = map(lambda x: 0, range(len(a[0])))
  if solveStep(a, b, result, 0):
    return result
  else:
    return None

matrix_a = [[0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0],
 [0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],
 [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0],
 [0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
 [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0],
 [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0],
 [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]]

vector_b = [2, 0, 2, 1, 2, 1, 1, 1, 1, 1, 1]

print solve(matrix_a, vector_b)

Это должно было проверить только 1325 возможных векторов с вашим входом, что намного меньше, чем все возможные результаты (67 миллионов).В худшем случае, конечно, 67 миллионов тестов.

1 голос
/ 16 августа 2011

Пусть Ax = b будет системой, а A|b будет расширенной матрицей, тогда

Есть 3 варианта

  1. Нет решений, если rank(A) < rank(A|b)
  2. Толькоодно решение, если rank(A) = rank(A|b) = n
  3. Бесконечное число решений, если rank(A) = rank(A|b) < n

, где n - число неизвестных.

...