0

I am having trouble understanding the documentation of np.einsum(). How are subscripts interpreted?

I am trying to write np.einsum('a...c,b...c', Y, conj(Y)) where Y is a matrix of shape C, F, T on the original python. Also, due to previous implementation differences my MATLAB Y is of size [F, T, C].

What does 'a...c,b...c' index in each component? I am confused.

How do I write the same instructions in MATLAB?

Dev-iL
  • 23,742
  • 7
  • 57
  • 99
havakok
  • 1,185
  • 2
  • 13
  • 45
  • 1
    Did you look for questions on this topic on SO? One example - [Element-wise matrix multiplication for multi-dimensional array](https://stackoverflow.com/q/55913093/3372061) – Dev-iL Feb 20 '20 at 10:18
  • I did but I missed your link. I am reading through it now. – havakok Feb 20 '20 at 10:21
  • It does not help me. For example, I do not understand why is 'mki' is `5,3,1`. what do these letters represent? are these indices of 'A'? If so, Arent they supposed to be `2,2,2000`? – havakok Feb 20 '20 at 10:29
  • Those letter represent the different dimension of the array. – obchardon Feb 20 '20 at 10:32
  • 1
    If you missed it, you must not have used the "search by tag" feature, [here's an example](https://stackoverflow.com/questions/tagged/matlab+numpy-einsum) of searching for [tag:matlab]+[tag:numpy-einsum] (which returns 5 question including this one and the one I linked). `5,3,1` is the order of dimensions. So if one has a matrix of size `[20,3,500,100,4,2]` it refers to the dimension of size 2, then the one of size 100 and finally the one of size 3 (considering python's 0-based indexing). – Dev-iL Feb 20 '20 at 10:32
  • So, does `a..c` refers to dimensions `0,1,2` and `b...c` refers to dimensions `1,2`? – havakok Feb 20 '20 at 10:38
  • 1
    Both `Y` and `Y.conj()` are 3d arrays of the same shape, so `a...c` is an implicit `aic` and `b...c` is an implicit `bic`, and the missing `->...` makes it equivalent to `aic, bic -> iab`. This in turn is the same as `(Y.transpose(1, 2, 0)[..., None] * Y.conj().transpose(1, 2, 0)[:, :, None, :]).sum(1)`: turn the first array's `[F, T, C]` shape into `[T, C, F, 1]` and the second into `[T, C, 1, F]` which broadcast to `[T, C, F, F]`, finally sum over the second dimension to get `[T, F, F]` which is hopefully exactly what you get. I'd be explicit in the `einsum` call to make it readable. – Andras Deak -- Слава Україні Feb 20 '20 at 14:28

1 Answers1

6

Quoting from the einsum documentation page:

To enable and control broadcasting, use an ellipsis. Default NumPy-style broadcasting is done by adding an ellipsis to the left of each term, like np.einsum('...ii->...i', a). To take the trace along the first and last axes, you can do np.einsum('i...i', a), or to do a matrix-matrix product with the left-most indices instead of rightmost, one can do np.einsum('ij...,jk...->ik...', a, b).

Later, an example is given:

>>> a = np.arange(25).reshape(5,5)
>>> np.einsum('...j->...', a)
array([ 10,  35,  60,  85, 110])

The equivalent MATLAB code to this example would be:

>> a = reshape(0:24, [5,5]).';
>> sum(a,2).'
ans =
    10    35    60    85   110

Several things to note:

  1. The ellipsis operator (...) should not be understood as "range", but as "whatever needs to be there".
  2. "Broadcasting" refers to automatic replication of an array along the relevant dimension so that the mathematical operation is defined. This is a feature that exists in MATLAB since R2016b (called "implicit expansion").
  3. You may notice several transpose operations (.') in the MATLAB equivalent. This is because numpy arrays are row-major while MATLAB array are column-major. Practically, while the underlying data has the same sequential order, a numpy array appears transposed compared to MATLAB. The transpositions were done so that arrays appear the same in intermediate stages.

Another example from these docs is:

>>> a = np.arange(6).reshape((3,2))
>>> b = np.arange(12).reshape((4,3))
>>> np.einsum('ki,jk->ij', a, b)
array([[10, 28, 46, 64],
       [13, 40, 67, 94]])
>>> np.einsum('ki,...k->i...', a, b)
array([[10, 28, 46, 64],
       [13, 40, 67, 94]])
>>> np.einsum('k...,jk', a, b)
array([[10, 28, 46, 64],
       [13, 40, 67, 94]])

Which can be written as follows in MATLAB:

A = reshape(0:5, [2 3]).';
B = reshape(0:11, [3 4]).';

A.' * B.'
permute(sum( permute(A, [3 1 2]) .* B,2), [3 1 2])
shiftdim(sum( shiftdim(A, -1) .* B, 2), 2)

Several things to note:

  1. When going from np.einsum('ki,jk->ij', a, b) to np.einsum('ki,...k->i...', a, b) you can see that the j-th dimension is replaced with .... The fact that both these examples have a -> in them, means it's in explicit mode.
  2. When going from np.einsum('ki,jk->ij', a, b) to np.einsum('k...,jk', a, b), you can see that now the i-th dimension is replaced with .... The omission of ->...j simply demonstrates implicit mode (where the output dimensions are ordered alphabetically).
Dev-iL
  • 23,742
  • 7
  • 57
  • 99