多维数组-如何解决许多线性方程组的超定系统使用矢量编码

Question

我需要解一个线性方程组Lx=b，其中x总是一个向量(3x1数组)，L是一个Nx3数组，b是一个Nx1向量。N通常在4到10之间。我使用

解决这个问题没有问题。

scipy.linalg.lstsq (L, b)

然而，我需要这样做很多次(比如200x200=40000次)，因为x实际上是与图像中的每个像素相关的东西。所以x实际上存储在PxQx3数组中，其中P和Q大概是200-300，最后一个数字"3"指的是向量x，现在我只是循环遍历每一列和每一行，逐个求解方程。我如何有效地求解这40000个不同的超定线性方程组，可能使用一些矢量化技术或其他特殊方法?

谢谢

Answer 1

您可以通过利用numpy.linalg例程的矩阵堆栈特性来获得一些速度。这还不能在numpy.linalg.lstsq中工作，但是numpy. linear .svd可以，所以您可以自己实现lstsq:

import numpy as np

def stacked_lstsq(L, b, rcond=1e-10):
    """
    Solve L x = b, via SVD least squares cutting of small singular values
    L is an array of shape (..., M, N) and b of shape (..., M).
    Returns x of shape (..., N)
    """
    u, s, v = np.linalg.svd(L, full_matrices=False)
    s_max = s.max(axis=-1, keepdims=True)
    s_min = rcond*s_max
    inv_s = np.zeros_like(s)
    inv_s[s >= s_min] = 1/s[s>=s_min]
    x = np.einsum('...ji,...j->...i', v,
                  inv_s * np.einsum('...ji,...j->...i', u, b.conj()))
    return np.conj(x, x)

def slow_lstsq(L, b):
    return np.array([np.linalg.lstsq(L[k], b[k])[0]
                     for k in range(L.shape[0])])    

def test_it():
    b = np.random.rand(1234, 3)
    L = np.random.rand(1234, 3, 6)
    x = stacked_lstsq(L, b)
    x2 = slow_lstsq(L, b)
    # Check
    print(x.shape, x2.shape)
    diff = abs(x - x2).max()
    print("difference: ", diff)
    assert diff < 1e-13

test_it()

一些时间表明堆叠版本在这里快了6倍左右，对于这个问题的大小。是否值得这么麻烦取决于这个问题。