Блог

Я пытаюсь переписать следующий фрагмент кода Matlab о внешнем продукте матриц в код 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?

Во всех кровавых подробностях (так как моя память 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). По крайней мере, в октаве это то же самое, что

 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/Октава 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]]]])
 

или без дополнительного размера 1 размер X :

 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])
 

сравните это с октавой равеля

 >> 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]]])
 

равель-это то же самое, что Октава:

 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)