如何使用矢量化代码解决许多超定线性方程组? [英] how to solve many overdetermined systems of linear equations using vectorized codes?

问题描述

我需要求解线性方程组Lx = b，其中x始终是向量(3x1数组)，L是Nx3数组，b是Nx1向量. N通常在4到10之类的范围内.使用

I need to solve a system of linear equations Lx=b, where x is always a vector (3x1 array), L is an Nx3 array, and b is an Nx1 vector. N usually ranges from 4 to something like 10. I have no problems solving this using

scipy.linalg.lstsq(L，b)

scipy.linalg.lstsq(L,b)

但是，我需要做很多次(大约200x200 = 40000次)，因为x实际上是与图像中每个像素相关联的东西.因此，x实际上存储在PxQx3数组中，其中P和Q类似于200-300，最后一个数字'3'表示向量x.现在，我只遍历每一行每一行，并逐一求解方程.如何有效地求解这40000个不同的线性方程组，可能使用某些矢量化技术或其他特殊方法?

However, I need to do this many times (something like 200x200=40000 times) as x is actually something associated with each pixel in an image. So x is actually stored in an PxQx3 array where P and Q is something like 200-300, and the last number '3' refers to the vector x. Right now I just loop through each column and row and solve the equation one-by-one .How do I solve those 40000 different overdetermined systems of linear equations efficiently, probably using some vectorization techniques or other special methods?

谢谢

推荐答案

通过使用 numpy .linalg.svd 可以，因此您可以自己实现lstsq:

You can gain some speed by making use of the stack of matrices feature of numpy.linalg routines. This doesn't yet work for numpy.linalg.lstsq, but numpy.linalg.svd does, so you can implement lstsq yourself:

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倍， 对于那个问题的大小.麻烦是否值得取决于 问题.

Some timing suggests the stacked version is around 6x faster here, for that problem size. Whether it's worth the trouble depends on the problem.

这篇关于如何使用矢量化代码解决许多超定线性方程组?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！