How to remove loops over inner()
Dear all: After using numpy for several weeks, I am very happy about it and deeply impressed about the performance improvements it brings in my python code. Now I have stumbled upon a problem, where I cannot use numpy to eliminate all my loops in python. Currently the return value of inner(a, b) is defined as inner(a, b)[I, J] = sum_k a[I, k] * b[J, k], for some super indices I and J. Somewhat more general is the tensordot() function that allows to specify over which axes K is summed over. However, if I understand numpy correctly, the following more general version is currently missing: inner(a, b, keep_axis=0)[H, I, J] = sum_k a[H, I, k] * b[H, J, k]. Here H would be an additional super index (specified via the keep_axis keyword), on which no outer product is taken, i.e., the same index is used for a[] and b[]. This more general definition would allow elimination of an extra level of loops. For example, I wish to calculate the following a = rand(200, 5, 2) b = rand(200, 4, 2) r = empty(a.shape[:1] + b.shape[1:1]) for h in range(a.shape[0]): r[h] = inner(a[h], b[h]) How could I eliminate the loop? It would be great if there would be the mentioned generalized version of the inner() [or tensordot()] function, since it would eliminate this loop and make my code much faster. What are your opinions? Would such a feature be desirable (or is it already implemented)? Thank you, Best, Hansres
On Nov 26, 2007 2:30 PM, HansAndreas Engel <engel@physics.harvard.edu> wrote:
Dear all:
After using numpy for several weeks, I am very happy about it and deeply impressed about the performance improvements it brings in my python code. Now I have stumbled upon a problem, where I cannot use numpy to eliminate all my loops in python.
Currently the return value of inner(a, b) is defined as inner(a, b)[I, J] = sum_k a[I, k] * b[J, k], for some super indices I and J. Somewhat more general is the tensordot() function that allows to specify over which axes K is summed over.
However, if I understand numpy correctly, the following more general version is currently missing: inner(a, b, keep_axis=0)[H, I, J] = sum_k a[H, I, k] * b[H, J, k]. Here H would be an additional super index (specified via the keep_axis keyword), on which no outer product is taken, i.e., the same index is used for a[] and b[].
This more general definition would allow elimination of an extra level of loops. For example, I wish to calculate the following a = rand(200, 5, 2) b = rand(200, 4, 2) r = empty(a.shape[:1] + b.shape[1:1]) for h in range(a.shape[0]): r[h] = inner(a[h], b[h]) How could I eliminate the loop? It would be great if there would be the mentioned generalized version of the inner() [or tensordot()] function, since it would eliminate this loop and make my code much faster.
What are your opinions? Would such a feature be desirable (or is it already implemented)?
Essentially, you want to operate on a stack of two dimensional arrays, correct? I'd be mildly supportive of something like this for tensordot; I'd prefer more descriptive name for keep_axis, but I don't know what it would be off the top of my head. In any event it should be XXX_axes and optionally take a sequence of axes so that more than one can be ignored. You could trivially build more specific functions on top of tensordot, so I don't see that inner needs to be changed as it's basically a convenience function anyway.  . __ . \ . . tim.hochberg@ieee.org
Essentially, you want to operate on a stack of two dimensional arrays, correct?
Yes, this is correct  and I also think that one should be able to provide a list of axes to be ignored.
I'd be mildly supportive of something like this for tensordot; I'd prefer more descriptive name for keep_axis, but I don't know what it would be off the top of my head. In any event it should be XXX_axes and optionally take a sequence of axes so that more than one can be ignored. You could trivially build more specific functions on top of tensordot, so I don't see that inner needs to be changed as it's basically a convenience function anyway.
Good point; indeed a generalization of tensordot() would be sufficient for me! I don't know the best naming here, ignore_axes might be a bit more descriptive. Would this feature only require a few lines of code? Would this be easy to implement for an expert on the intricacies of the numpy internals? Best, Hansres
participants (2)

HansAndreas Engel

Timothy Hochberg