Python 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 ，它不要求矩阵是非奇异的，并且与逆矩阵类似，理论上可以恢复 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

Cutoff for small singular values. Singular values smaller (in modulus) than rcond * largest_singular_value (again, in modulus) are set to zero. Broadcasts against the stack of matrices

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 向上计数，也有大量的组合解决方案。

如果有足够的时间，下面的脚本会找到所有这些内容。在实践中，我没有考虑超出 32x32 矩阵的情况，在下面所示的运行中产生了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