Matlab‘外积’的Python实现

Question

我正在尝试将以下关于矩阵外积的 Matlab 代码 sn-p 重写为 python 代码，

function Y = matlab_outer_product(X,x)
A = reshape(X, [size(X) ones(1,ndims(x))]);
B = reshape(x, [ones(1,ndims(X)) size(x)]);
Y = squeeze(bsxfun(@times,A,B));
end

我对python代码的一对一翻译如下（考虑numpy数组和matlab矩阵的形状如何排列），

def python_outer_product(X, x):
    X_shape = list(X.shape)
    x_shape = list(x.shape)
    A = X.reshape(*list(np.ones(np.ndim(x),dtype=int)),*X_shape)
    B = x.reshape(*x_shape,*list(np.ones(np.ndim(X),dtype=int)))
    Y = A*B
    return Y.squeeze()

然后尝试输入，例如，

matlab_outer_product([1,2],[[3,4];[5,6]])
python_out_product(np.array([[1,2]], np.array([[3,4],[5,6]])))

输出不完全匹配。在matlab中，它输出

output(:,:,1) = [[3,5];[6,10]]
output(:,:,2) = [[4,6];[8,12]]

在python中，它输出

output = array([
       [[ 3,  6],
        [ 4,  8]],

       [[ 5, 10],
        [ 6, 12]]
])

它们几乎相同，但并不完全相同。我想知道代码有什么问题以及如何更改python代码以匹配matlab输出？

Answer 1

血淋淋的细节（因为我的 MATLAB 内存太旧了）：

八度

>> X = [1,2];
>> x = [[3,4];[5,6]];
>> A = reshape(X, [size(X) ones(1,ndims(x))]);
>> B = reshape(x, [ones(1,ndims(X)) size(x)]);
>> A
A =

   1   2

>> B
B =

ans(:,:,1,1) =  3
ans(:,:,2,1) =  5
ans(:,:,1,2) =  4
ans(:,:,2,2) =  6

>> bsxfun(@times,A,B)
ans =

ans(:,:,1,1) =

   3   6

ans(:,:,2,1) =

    5   10

ans(:,:,1,2) =

   4   8

ans(:,:,2,2) =

    6   12

>> squeeze(bsxfun(@times,A,B))
ans =

ans(:,:,1) =

    3    5
    6   10

ans(:,:,2) =

    4    6
    8   12

您从 (1,2) 和 (2,2) 开始，将第二个扩展为 (1,1,2,2)。 bsxfun 产生一个 (1,2,2,2)，它被压缩到 (2,2,2)。

A 被X 重整为[1 2 1 1]，但是两个外部尺寸 1 尺寸被挤出，导致没有变化。

这个 MATLAB 外部有点复杂，使用 bsxfun 执行 (1,2,1,1) 与 (1,1,1,2) 的元素乘法。至少在 Octave 中是一样的

A.*B

在 numpy 中

In [77]: X
Out[77]: array([[1, 2]])    # (1,2)
In [78]: x
Out[78]: 
array([[3, 4],              # (2,2)
       [5, 6]])

请注意，MATLAB/Octave x 在展平时具有元素 (3,5,4,6)，而 numpy ravel 为 [3,4,5,6]。

在 numpy 中我可以简单地做：

In [79]: X[:,:,None,None]*x
Out[79]: 
array([[[[ 3,  4],          (1,2,2,2)
         [ 5,  6]],

        [[ 6,  8],
         [10, 12]]]])

或者没有X的额外尺寸1维度：

In [84]: (X[0,:,None,None]*x)
Out[84]: 
array([[[ 3,  4],
        [ 5,  6]],

       [[ 6,  8],
        [10, 12]]])

In [85]: (X[0,:,None,None]*x).ravel()
Out[85]: array([ 3,  4,  5,  6,  6,  8, 10, 12])

将其与 Octave ravel 进行比较

>> squeeze(bsxfun(@times,A,B))(:)'
ans =

    3    6    5   10    4    8    6   12

我们可以向 numpy 添加转置

In [96]: (X[0,:,None,None]*x).transpose(2,1,0).ravel()
Out[96]: array([ 3,  6,  5, 10,  4,  8,  6, 12])
In [97]: (X[0,:,None,None]*x).transpose(2,1,0)
Out[97]: 
array([[[ 3,  6],
        [ 5, 10]],

       [[ 4,  8],
        [ 6, 12]]])

至少在 numpy 中，我们可以通过多种方式调整维度顺序，因此我不会尝试提出最佳建议。我仍然认为编写对 numpy 来说“自然”的代码比盲目地匹配 MATLAB 顺序要好。

再试一次

我在上面意识到 MATLAB 只是在使用 A*.B (1,2,1,1) 数组 (1,1,1,2)，其中额外的 1 被添加到“广播”中。

使用转置到最外层的相同维度（以numpy开头）

In [5]: X = X.T; x = x.T
In [6]: X.shape
Out[6]: (2, 1)
In [7]: x.shape
Out[7]: (2, 2)
In [8]: x
Out[8]: 
array([[3, 5],
       [4, 6]])
In [9]: x.ravel()
Out[9]: array([3, 5, 4, 6])   # compare with MATLAB (:)'

具有相同维度扩展的元素乘法：

In [10]: X[None,None,:,:]*x[:,:,None,None]
Out[10]: 
array([[[[ 3],
         [ 6]],

        [[ 5],
         [10]]],


       [[[ 4],
         [ 8]],

        [[ 6],
         [12]]]])
In [11]: _.shape
Out[11]: (2, 2, 2, 1)         # compare with octave (1,2,2,2)
In [12]: __.squeeze()
Out[12]: 
array([[[ 3,  6],
        [ 5, 10]],

       [[ 4,  8],
        [ 6, 12]]])

ravel 和 Octave 一样：

In [13]: ___.ravel()
Out[13]: array([ 3,  6,  5, 10,  4,  8,  6, 12])

expand_dims 可以用来代替索引。在内部它使用reshape:

In [15]: np.expand_dims(X,(0,1)).shape
Out[15]: (1, 1, 2, 1)
In [16]: np.expand_dims(x,(2,3)).shape
Out[16]: (2, 2, 1, 1)