如何向量化这个 linalg.lstq() 操作？

Question

我正在尝试使用Python3和NumPy实现多频相位展开算法。我有 7 个形状为

(1080, 1920)

的单通道（灰度）图像。将它们沿第三轴堆叠后，我得到

(1080, 1920, 7)

。

我有形状为

A

的移位矩阵

(7, 7)

（对于所有像素都是固定的）。对于所有像素，我都有不同的强度值（形状为

r

的数组

(1, 7)

）。

要通过最小化 L2 范数为 ||r - Au|| 来求解每个像素，我可以这样做：

# Just for illustration
A = np.random.randn(7, 7)
r = np.random.randn(7, 1)

# Solved for each pixel location
u = np.linalg.lstsq(a = A, b = r, rcond = None)

这可以使用 2 个 for 循环来实现：

for y in range(0, image_height - 1):
    for x in range(0, image_width - 1):
        # code here

这是低效的。如何将其编写为高效的 NumPy 代码？

Answer 1

IIUC 你有一个

A

矩阵和

1080*1920

RHS。 NumPy 的

lstsq

函数针对这种情况进行了矢量化。

因此我们将您的

img

重塑为形状

(1080,1920, 7)

至

import numpy as np
rng = np.random.default_rng(234892359843)
img = rng.random((1080, 1920, 7))
A = rng.random((7, 7))
b = img.reshape((-1, 7)).T  # reshape to (7, 1080*1920)
x, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
x = (x.T).reshape(img.shape)  # return to original shape

然后对于任何有效索引

i

，

j

：

np.testing.assert_allclose(x[i, j], np.linalg.lstsq(A, img[i, j], rcond=None)[0])

我假设你的矩阵

A

是奇异的；否则您可以使用与

solve

相同的代码。

那么这应该成立：

np.testing.assert_allclose(A @ x[i, j], img[i, j])