Broadcasting in Python with permutations
<p>I understand that <code>transpose</code> on an <code>ndarray</code> is intended to be the equivalent of matlab's <code>permute</code> function however I have a specific usecase that doesn't work simply. In matlab I have the following:</p>

<pre><code>C = @bsxfun(@times, permute(A,[4,2,5,1,3]), permute(B, [1,6,2,7,3,4,5])
</code></pre>

<p>where A is a 3D tensor of shape NxNxM and B is a 5D tensor of shape NxNxMxPxP. The above function is meant to vectorize looped kronecker products. I'm assuming that Matlab is adding 2 singleton dimensions for both A and B which is why it's able to rearrange them. I'm looking to port this code over to Python <strike>but I don't think it has the capability of adding these extra dimensions.</strike>. I found <a href="http://stackoverflow.com/questions/17835121/is-this-the-best-way-to-add-an-extra-dimension-to-a-numpy-array-in-one-line-of-c">this</a> which successfully adds the extra dimensions however the broadcasting is not working the same matlab's <code>bsxfun</code>. I have attempted the obvious translation (yes I am using numpy for these <code>ndarray</code>'s and functions):</p>

<pre><code>A = A[...,None,None]
B = B[...,None,None]
C = transpose(A,[3,1,4,0,2])*transpose(B,[0,5,1,6,2,3,4])
</code></pre>

<p><strike>and I get the following error:</p>

<pre><code>return transpose(axes)
ValueError: axes don't match array
</code></pre>

<p>My first guess is to do a <code>reshape</code> on A and B to add in those singleton dimensions?</strike></p>

<p>I now get the following error:</p>

<pre><code>mults = transpose(rho_in,[3,1,4,0,2])*transpose(proj,[0,5,1,6,2,3,4])
ValueError: operands could not be broadcast together with shapes (1,9,1,9,8) (9,1,9,1,8,40,40)
</code></pre>

<p>EDIT: Amended my question to be less about adding singleton dimensions but more about correctly broadcasting this matlab multiplication in python.</p>

<p>The huge difference between MATLAB and numpy is that the former uses column-major format for its arrays, while the latter row-major. The corollary is that implicit singleton dimensions are handled differently.</p>

<p>Specifically, MATLAB explicitly ignores trailing singleton dimensions: <code>rand(3,3,1,1,1,1,1)</code> is actually a <code>3x3</code> matrix. Along these lines, you can use <code>bsxfun</code> to operate on two arrays if their <em>leading</em> dimensions match: <code>NxNxM</code> is implicitly <code>NxNxMx1x1</code>  which is compatible with <code>NxNxMxPxP</code>.</p>

<p>Numpy, on the other hand, <a href="http://docs.scipy.org/doc/numpy-1.10.1/user/basics.broadcasting.html#general-broadcasting-rules" rel="nofollow">allows implicit singletons <em>up front</em></a>. You need to <code>permute</code> your arrays in a way that their <em>trailing</em> dimensions match up, for instance shape <code>(40,40,9,1,9,1,8)</code> with shape <code>(1,9,1,9,8)</code>, and the result should be of shape <code>(40,40,9,9,9,9,8)</code>.</p>

<p>Dummy example:</p>

<pre><code>&gt;&gt;&gt; import numpy as np
&gt;&gt;&gt; (np.random.rand(40,40,9,1,9,1,8)+np.random.rand(1,9,1,9,8)).shape
(40, 40, 9, 9, 9, 9, 8)
</code></pre>

<p>Note that what you're trying to do can probably be done using <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html" rel="nofollow"><code>numpy.einsum</code></a>. I suggest taking a closer look at that. An example of what I mean: from your question I gathered that you want to perform this: take elements <code>A[1:N,1:N,1:M]</code> and <code>B[1:N,1:N,1:M,1:P,1:P]</code> and construct a new array <code>C[1:N,1:N,1:N,1:N,1:M,1:P,1:P]</code> such that</p>

<pre><code>C[i1,i2,i3,i4,i5,i6,i7] = A[i2,i4,i5]*B[i1,i3,i5,i6,i7]
</code></pre>

<p>(your specific index order might vary). If this is correct, you can indeed use <code>numpy.einsum()</code>:</p>

<pre><code>&gt;&gt;&gt; a = np.random.rand(3,3,2)
&gt;&gt;&gt; b = np.random.rand(3,3,2,4,4)
&gt;&gt;&gt; np.einsum('ijk,lmkno-&gt;limjkno',a,b).shape
(3, 3, 3, 3, 2, 4, 4)
</code></pre>

<p>Two things should be noted, though. Firstly, the above operation will be very memory-intense, which should be expected for cases of vectorization (where you usually win CPU time at the expense of memory need). Secondly, you should seriously consider rearranging the data model when you port your code. The reason that broadcasting works differently in the two languages is intricately connected to the column-major/row-major difference. This also implies that in MATLAB you should work with <em>leading</em> indices first, since <code>A(:,i2,i3)</code> corresponds to a contiguous block of memory, while <code>A(i1,i2,:)</code> does not. Conversely, in numpy <code>A[i1,i2,:]</code> is contiguous, while <code>A[:,i2,i3]</code> is not.</p>

<p>These considerations suggest that you should set up the logistics of your data such that vectorized operations preferably work with leading indices in MATLAB and trailing indices in numpy. You could still use <code>numpy.einsum</code> to perform the operation itself, however your dimensions should be in a different (possibly reverse) order compared to MATLAB, at least if we assume that both versions of the code use an optimal setup.</p>

python matlab vectorization broadcasting permute