Питон numpy: linalg.pinv () слишком неточно

Question

В последнее время я работаю с numpy матрицами в алгоритме и столкнулся с проблемой:

Я использую всего 3 матрицы.

m1 = [[  3   2   2 ...   3   2   3]
      [  3   3   3 ...   2   2   2]
      [500 501 502 ... 625 626 627]
      ...
      [623 624 625 ... 748 749 750]
      [624 625 626 ... 749 750 751]
      [625 626 627 ... 750 751 752]]

m1является (128,128) единичной квадратной матрицей.Первые две строки представляют собой, казалось бы, случайные последовательности 2 и 3.Следующие строки заполняются алгоритмически, начиная с 500, добавляя по одному для каждой строки и для каждого столбца, начиная с третьей строки, первого столбца.

m2 = [[  2   3 500 ... 623 624 625]
      [  2   2 500 ... 623 624 625]
      [  3   2 500 ... 623 624 625]
      ...
      [  2   3 500 ... 623 624 625]
      [  2   2 500 ... 623 624 625]
      [  3   2 500 ... 623 624 625]]

m2 также является (128,128) сингулярной квадратной матрицей.На этот раз случайные последовательности относятся к первым двум столбцам.Остальная часть каждой строки заполнена 500, 501, 502, 503 и т. Д.

m3 = [[     790      784   157500 ...   196245   196560   196875]
      [     804      811   161000 ...   200606   200928   201250]
      [  180501   180411 36064000 ... 44935744 45007872 45080000]
      ...
      [  219861   219771 43936000 ... 54744256 54832128 54920000]
      [  220181   220091 44000000 ... 54824000 54912000 55000000]
      [  220501   220411 44064000 ... 54903744 54991872 55080000]]

m3 = m1*m2

Итак, я хотел восстановить m2используя m1 и m3.Теоретически все, что мне нужно было сделать, это выполнить следующий код m2 = (m1**-1)*m3.К сожалению, из-за того, что m1 была единственной матрицей, было невозможно рассчитать ее обратную величину, и, даже если бы это было возможно, матрица была слишком большой, вызывая многочисленные числовые неточности.

Вместо этого я решилиспользуйте Moore-Penrose Inverse из m1, который не требует, чтобы матрица была неособой, и, как и Inverse, делает теоретически возможным восстановление m2, используя np.linalg.pinv(m1) * m3.

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

[[  2.46207616   2.48959603 500.         ... 623.         624.
  625.        ]
 [  2.38612549   2.61197086 500.         ... 623.         624.
  625.        ]
 [  2.38711085   2.6125801  500.         ... 623.         624.
  625.        ]
 ...
 [  2.61998539   2.37184747 500.         ... 623.         624.
  625.        ]
 [  2.54403472   2.4942223  500.         ... 623.         624.
  625.        ]
 [  2.62195611   2.37306595 500.         ... 623.         624.
  625.        ]]

Как видите, вся «заполнитель» части m1 рассчитан правильно, с этим проблем нет.Однако, похоже, что с первыми двумя столбцами возникли проблемы, и округление чисел до 2 и 3 дает мне неправильный m2.

Я ищу способ сделать метод np.linalg.pinv()точнее с его вычислениями с плавающей точкой, чтобы он мог получить правильные значения для последовательностей, поскольку они очень важны.

Проведя некоторое исследование, я узнал, что у np.linalg.pinv() есть аргумент, называемый rcond, описываемый какследующее:

rcond: (…) array_like of float

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

rcond по умолчанию установлено на 1e-15.Я думал, что сокращение этого числа еще больше может помочь с неточностью.1e-16 было недостаточно, и, начиная с 1e-17, я получаю очень странные значения, такие как:

[[ 3.000e+00  3.000e+00  5.000e+02 ...  6.230e+02  6.240e+02  6.250e+02]
 [ 1.100e+01  4.000e+00  1.840e+02 ...  1.722e+03  2.032e+03  1.831e+03]
 [-3.000e+00 -5.000e+00 -4.030e+02 ... -1.232e+03 -7.400e+02 -1.272e+03]
 ...
 [ 1.100e+01  1.200e+01  2.164e+03 ...  4.030e+03  4.872e+03  1.873e+03]
 [-1.200e+01 -9.000e+00 -1.618e+03 ... -3.240e+03 -2.519e+03 -4.167e+03]
 [ 2.000e+01  2.600e+01  4.535e+03 ...  5.165e+03  5.881e+03  5.189e+03]]

Так что, по сути, я застрял, я не знаю, какповысить точность.Хуже всего то, что у меня есть модуль, который может значительно увеличить точность с плавающей точкой, он называется mpmath, а также есть матрицы, которые, как мне кажется, работают даже лучше с моим алгоритмом как numpy.Но mpmath не имеет метода для вычисления псевдообращений, а numpy не адаптирует свою собственную точность с плавающей запятой к значению, установленному с помощью mpmath.

Не могли бы вы предложить какие-либо предложения, которые я мог бы попробоватьполучить правильный m2, используя псевдообратный метод?

Answer 1

Ваши проблемы не имеют ничего общего с точностью или нет pinv.

Если вы заметите, что ваши матрицы имеют недостаточный ранг, m1 имеет ранг 4 или менее, m2 ранг 3или менее.Следовательно, ваша система m1@x = m3 крайне недоопределена, и восстановить невозможно m2.

Даже если мы добавим все, что знаем, о структуре m2, то есть первых двух столбцах 3 и2, остальные 500, считая вверх, есть комбинаторно большое количество решений.

Сценарий ниже находит их все, если есть достаточно времени.На практике я не смотрел дальше матриц 32х32, которые в приведенном ниже прогоне дали 15093381006 различных действительных реконструкций m2', которые удовлетворяют m1@m2' = m3, и структурных ограничений, которые я только что упомянул.

import numpy as np
import itertools

def make(n, offset=500):
    offset -= 2
    res1 = np.add.outer(np.arange(n), np.arange(offset, offset+n))
    res1[:2] = np.random.randint(2, 4, (2, n))
    res2 = np.add.outer(np.zeros(n, int), np.arange(offset, offset+n))
    res2[:, :2] = np.random.randint(2, 4, (n, 2))
    return res1, res2

def subsets(L, n, mn, mx, prepend=[]):
    if n == 0:
        if mx >= mn:
            yield prepend
    elif n == 1:
        for l in L[L.searchsorted(mn):L.searchsorted(mx, 'right')]:
            yield prepend + [l]
    else:
        ps = L.cumsum()
        ps[n:] -= ps[:-n]
        ps = ps[n-1:]
        for i in range(L.searchsorted(mn - np.sum(L[len(L)-n+1:])),
                       ps.searchsorted(mx, 'right')):
            yield from subsets(L[i+1:], n-1, mn - L[i], mx - L[i],
                               prepend = prepend + [L[i]])

def solve(m1, m3, ci=0, offset=500):
    n, n = m1.shape
    col = m3.T[ci]
    n3s = col[3] - col[2] - 2 * n
    six = col[2] - offset * (col[3] - col[2]) - n * (n-1)
    idx = np.lexsort(m1[:2])
    m1s = m1[:2, idx]
    sm = m1s[1].searchsorted(2.5)
    sl = m1s[0, :sm].searchsorted(2.5)
    sr = sm + m1s[0, sm:].searchsorted(2.5)
    n30 = n - sl - sr + sm
    n31 = n - sm
    n330 = col[0] - 4*n - 2*n3s - 2*n30
    n331 = col[1] - 4*n - 2*n3s - 2*n31
    for n333 in range(max(0, n330 - sm + sl, n331 - sr + sm),
                      min(n - sr, n330, n331) + 1):
        n332 = n330 - n333
        n323 = n331 - n333
        n322 = n3s - n332 - n323 - n333
        mx333 = six - idx[sl:sl+n332].sum() - idx[sm:sm+n323].sum() \
                - idx[:n322].sum()
        mn333 = six - idx[sm-n332:sm].sum() - idx[sr-n323:sr].sum() \
                - idx[sl-n322:sl].sum()
        for L333 in subsets(idx[sr:], n333, mn333, mx333):
            mx332 = six - np.sum(L333) - idx[sm:sm+n323].sum() \
                - idx[:n322].sum()
            mn332 = six - np.sum(L333) - idx[sr-n323:sr].sum() \
                - idx[sl-n322:sl].sum()
            for L332 in subsets(idx[sl:sm], n332, mn332, mx332,
                                prepend=L333):
                mx323 = six - np.sum(L332) - idx[:n322].sum()
                mn323 = six - np.sum(L332) - idx[sl-n322:sl].sum()
                for L323 in subsets(idx[sm:sr], n323, mn323, mx323,
                                    prepend=L332):
                    ex322 = six - np.sum(L323)
                    yield from subsets(idx[:sl], n322, ex322, ex322,
                                       prepend=L323)

def recon(m1, m3, ci=0, offset=500):
    n, n = m1.shape
    nsol = nfp = 0
    REC = []
    for i3s in solve(m1, m3, ci, offset):
        rec = np.full(n, 2)
        rec[i3s] = 3
        if not np.all(m3.T[ci] == m1@rec):
            print('!', rec, m3.T[ci], m1@rec)
            nfp += 1
        else:
            nsol += 1
            REC.append(rec)
    print('col', ci, ':',  nsol, 'solutions,', nfp, 'false positives')
    return np.array(REC)

def full_recon(m1, m3, offset=500, subsample=10):
    n, n = m1.shape
    col0, col1 = (recon(m1, m3, i, offset) for i in (0, 1))
    yield col0.shape[0], col1.shape[0]
    if not subsample is None:
        col0, col1 = (col[np.random.choice(col.shape[0], subsample)]
                      if col.shape[0] > subsample else col
                      for col in (col0, col1))
    print('col 0', col0)
    print('col 1', col1)
    for c0, c1 in itertools.product(col0, col1):
        out = np.add.outer(np.zeros(n, int), np.arange(offset-2, offset+n-2))
        out[:, :2] = np.c_[c0, c1]
        yield out

def check(m1, m3, offset=500, subsample=10):
    for cnt, m2recon in enumerate(full_recon(m1, m3, offset, subsample)):
        if cnt == 0:
            tot0, tot1 = m2recon
            continue
        assert np.all(m3 == m1@m2recon)
    print(cnt, 'solutions verified out of', tot0, 'x', tot1, '=', tot0 * tot1)

Пример прогона:

>>> m1, m2 = make(32)
>>> check(m1, m1@m2)
col 0 : 133342 solutions, 0 false positives
col 1 : 113193 solutions, 0 false positives
col 0 [[2 2 3 2 2 2 3 2 3 3 3 3 3 2 3 3 2 2 2 2 2 3 3 3 2 2 2 2 2 2 2 2]
 [2 2 3 2 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 3 3 3 2 2 2 3 2 3 2 2 2 2]
 [2 3 3 2 3 3 2 2 2 3 2 3 3 2 2 2 2 3 3 3 2 2 2 2 3 2 2 2 2 2 2 3]
 [2 2 3 3 2 2 3 3 3 2 2 3 2 3 3 3 2 2 2 2 2 2 2 3 2 3 3 2 2 2 2 2]
 [2 3 3 3 2 3 3 2 2 2 2 3 2 2 3 3 2 2 3 2 2 3 2 2 2 2 2 2 3 3 2 2]
 [3 2 2 2 3 2 3 3 3 2 3 2 2 2 2 3 3 3 2 2 2 2 3 3 2 3 2 2 2 2 2 2]
 [2 2 2 3 3 3 2 2 3 3 3 3 3 2 2 2 2 3 2 2 2 3 2 2 2 3 2 2 3 2 2 2]
 [2 2 2 3 3 2 3 3 3 2 3 2 2 3 3 3 2 3 2 2 2 2 2 2 2 2 2 3 2 3 2 2]
 [3 2 3 3 2 2 3 2 2 3 2 3 3 2 2 2 3 2 3 2 2 2 2 3 2 3 2 2 3 2 2 2]
 [3 2 2 3 3 2 3 3 2 2 2 3 2 3 2 2 3 2 3 3 2 2 2 2 2 3 2 2 2 2 2 3]]
col 1 [[2 2 2 2 3 3 3 3 3 3 3 3 2 2 3 2 2 2 3 3 2 3 2 2 2 3 2 3 2 2 2 3]
 [2 3 3 2 3 3 2 2 2 3 3 3 2 2 3 3 2 3 2 3 2 2 3 2 2 2 2 3 3 2 2 3]
 [3 2 3 2 2 3 3 2 2 3 3 3 2 2 3 2 3 2 3 2 3 2 3 3 2 2 2 2 2 3 3 2]
 [2 2 3 2 3 3 3 2 3 3 2 3 2 3 2 3 2 3 2 2 3 2 3 2 2 3 2 3 2 2 2 3]
 [3 3 2 2 2 2 3 2 3 3 3 3 2 3 3 2 2 3 3 2 2 2 2 2 3 2 2 3 3 3 2 2]
 [3 3 2 3 2 2 2 3 3 3 2 2 2 3 2 2 2 3 3 3 3 2 3 3 2 2 2 3 3 2 2 2]
 [2 3 2 3 2 2 3 3 2 3 3 3 2 3 2 2 3 3 2 3 2 2 3 2 3 2 2 3 2 2 3 2]
 [2 3 2 3 2 3 3 3 3 3 3 2 2 2 2 2 3 2 3 2 2 2 3 2 2 3 3 2 3 2 2 3]
 [3 2 3 2 2 3 3 3 2 2 3 3 2 2 3 3 2 3 2 2 3 2 2 3 2 2 3 2 3 2 2 3]
 [3 2 3 3 2 3 2 2 2 3 3 3 2 3 3 2 2 2 3 2 3 2 2 3 2 2 3 2 2 2 3 3]]
100 solutions verified out of 133342 x 113193 = 15093381006