如何将 2D numpy 数组与 3D 数组矩阵相乘以得到 3D 数组？

Question

我正在解决光度立体问题，其中我有 "n" 个光源，每个光源有 3 个通道（红、绿、蓝）。因此光阵列的形状为 nx3：lights.shape = nx3 我有对应于每个照明条件的图像。图片尺寸为 hxw（高 x 宽），images.shape = nxhxw

我想将图像中的每个像素矩阵乘以一个形状为 3 x n 的矩阵，并得到另一个形状为 3xhxw 的数组，这些将是图像上每个像素的法向量。

形状：

图片 : (n_ims, h, w)
灯光：(n_ims, 3)

S = lights
S_pinv =  np.linalg.inv(S.T@S)@S.T  # pinv is pseudo inverse, S_pinv.shape : (n_ims,3)
b = S_pinv @ images  # I want (3xn @ nxhxw = 3xhxw)

但是我收到这个错误：

ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 100 is different from 3)

Answer 1

问题是 numpy 将多维数组视为矩阵堆栈，并且总是假定最后两个维度是线性 space 维度。这意味着点积将无法通过折叠 3d 数组的 first 维度来工作。

相反，您可以做的最简单的事情是将 3d 数组重新整形为 2d 数组，进行矩阵乘法，然后重新整形为 3d 数组。这也将利用优化的 BLAS 代码，这是 numpy 的一大优势。

import numpy as np 

S_pinv = np.random.rand(3, 4)
images = np.random.rand(4, 5, 6)

# error: 
# (S_pinv @ images).shape 
res_shape = S_pinv.shape[:1] + images.shape[1:]  # (3, 5, 6) 
res = (S_pinv @ images.reshape(images.shape[0], -1)).reshape(res_shape)
print(res.shape)  # (3, 5, 6)

所以我们做 (3,n) x (n, h*w) -> (3, h*w) 而不是 (3,n) x (n,h,w)，我们将其重塑回 (3, h, w)。重塑是免费的，因为这并不意味着对内存中的数据进行任何实际操作（只是对数组下的单个内存块的重新解释），正如我所说，适当的矩阵乘积在 numpy 中得到了高度优化。

既然您要求 更直观的 方法，这里有一个使用 numpy.einsum 的替代方法。它可能会更慢，但如果你稍微习惯它的符号，它会非常透明：

res_einsum = np.einsum('tn,nhw -> thw', S_pinv, images)
print(np.array_equal(res, res_einsum))  # True

此符号命名输入数组的每个维度：对于 S_pinv，第一和第二维度分别命名为 t 和 n，类似地 n、h 和 w 对应 images。输出设置为具有维度 thw，这意味着输出形状中不存在的任何剩余维度将在与输入数组相乘后相加。这正是您所需要的。

正如您在评论中指出的那样，您还可以转置数组，以便 np.dot 在正确的位置找到正确的维度。但这也会很慢，因为这可能会导致内存中的副本，或者至少是对数组的次优循环。

我使用以下定义进行了快速时序比较：

def reshaped(S_pinv, images): 
    res_shape = S_pinv.shape[:1] + images.shape[1:] 
    return (S_pinv @ images.reshape(images.shape[0], -1)).reshape(res_shape)

def einsummed(S_pinv, images): 
    return np.einsum('tn,nhw -> thw', S_pinv, images) 

def transposed(S_pinv, images): 
    return (S_pinv @ images.transpose(2, 0, 1)).transpose(1, 2, 0)

下面是使用 IPython 的内置 %timeit 魔法和一些更实际的数组大小的计时测试：

>>> S_pinv = np.random.rand(3, 100) 
... images = np.random.rand(100, 200, 300) 
... args = S_pinv, images 
... %timeit reshaped(*args) 
... %timeit einsummed(*args) 
... %timeit transposed(*args)                                          
5.92 ms ± 460 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
15.9 ms ± 190 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
44.5 ms ± 329 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Answer 2

答案是 np.swapaxes

import numpy as np

q= np.random.random([2, 5,5])
q.shape

w = np.random.random([3,2])
w.shape

w@q

我们有 ValueError 但是

import numpy as np

q= np.random.random([5, 2,5])
q.shape

w = np.random.random([3,2])
w.shape

res = (w@q).swapaxes(0,1)
res.shape # =[3, 5, 5]

Answer 3

一个简单的方法是 np.inner； inner 沿最后一个轴缩小并保留所有其他轴；因此它取决于一个完美匹配的转置：

n,h,w = 10,384,512
images = np.random.randint(1,10,(n,h,w))
S_pinv = np.random.randint(1,10,(n,3))

res_inr = np.inner(images.T,S_pinv.T).T
res_inr.shape
# (3, 384, 512)

同样，使用转置 matmul 实际上是正确的：

res_mml = (images.T@S_pinv).T
assert (res_mml==res_inr).all()

这两个似乎与@AndrasDeak 的einsum 方法大致同样快。

特别是，它们不如 reshaped matmul 快（不足为奇，因为单个直接 matmul 必须是目前最优化的操作之一）。他们以速度换取便利。

Answer 4

这基本上就是 np.einsum 的目的。

而不是：

b = S_pinv @ images

使用

b = np.einsum('ij, ikl -> jkl', S_pinv, images)

在这种情况下 i = n_ims、j = 3、k = h 和 l = w

因为是单次收缩，你也可以用np.tensordot()

b = np.tensordot(S_pinv.T, images, axes = 1)

或者，

b = np.tensordot(S_pinv, images, axes = ([0], [0]))

如何将 2D numpy 数组与 3D 数组矩阵相乘以得到 3D 数组？

How to matrix-multiply a 2D numpy array with a 3D array to give a 3D array?

python

numpy

numpy-ndarray

array-broadcasting