Nice function, and wonderful that it speeds some tasks up. <br><br>some feedback: the following notation is a little counter intuitive to me:<br><div><div>    &gt;&gt;&gt; np.einsum(&#39;i...-&gt;&#39;, a)</div><div>    array([50, 55, 60, 65, 70])</div>
<div>    &gt;&gt;&gt; np.sum(a, axis=0)</div>
<div>    array([50, 55, 60, 65, 70])<br>
</div>
Since there is nothing after the -&gt;, I expected a scalar not a vector. I might suggest &#39;i...-&gt;...&#39;<div><br>Just noticed also a typo in the doc:<br>
     order : &#39;C&#39;, &#39;F&#39;, &#39;A&#39;, or &#39;K&#39;<br>
<div>        Controls the memory layout of the output. &#39;C&#39; means it should</div>



        be Fortran contiguous. &#39;F&#39; means it should be Fortran contiguous,<br>should be changed to <br>
    order : &#39;C&#39;, &#39;F&#39;, &#39;A&#39;, or &#39;K&#39;
<div>        Controls the memory layout of the output. &#39;C&#39; means it should</div>

<div>        be C contiguous. &#39;F&#39; means it should be Fortran contiguous,<br>
</div>
<br></div><br>Hope this helps,<br>Jonathan<br></div><br><div class="gmail_quote">On Wed, Jan 26, 2011 at 2:27 PM, Mark Wiebe <span dir="ltr">&lt;<a href="mailto:mwwiebe@gmail.com" target="_blank">mwwiebe@gmail.com</a>&gt;</span> wrote:<br>



<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">I wrote a new function, einsum, which implements Einstein summation notation, and I&#39;d like comments/thoughts from people who might be interested in this kind of thing.<div>



<br></div><div>In testing it, it is also faster than many of NumPy&#39;s built-in functions, except for dot and inner.  At the bottom of this email you can find the documentation blurb I wrote for it, and here are some timings:</div>




<div><br></div><div><div>In [1]: import numpy as np</div><div>In [2]: a = np.arange(25).reshape(5,5)</div><div><br></div><div>In [3]: timeit np.einsum(&#39;ii&#39;, a)</div><div>100000 loops, best of 3: 3.45 us per loop</div>




<div>In [4]: timeit np.trace(a)</div><div>100000 loops, best of 3: 9.8 us per loop</div><div><br></div><div>In [5]: timeit np.einsum(&#39;ii-&gt;i&#39;, a)</div><div>1000000 loops, best of 3: 1.19 us per loop</div><div>In [6]: timeit np.diag(a)</div>




<div>100000 loops, best of 3: 7 us per loop</div><div><br></div><div>In [7]: b = np.arange(30).reshape(5,6)</div><div><br></div><div>In [8]: timeit np.einsum(&#39;ij,jk&#39;, a, b)</div><div>10000 loops, best of 3: 11.4 us per loop</div>




<div>In [9]: timeit np.dot(a, b)</div><div>100000 loops, best of 3: 2.8 us per loop</div><div><br></div><div>In [10]: a = np.arange(10000.)</div><div><div><br></div><div>In [11]: timeit np.einsum(&#39;i-&gt;&#39;, a)</div>




<div>10000 loops, best of 3: 22.1 us per loop</div><div>In [12]: timeit np.sum(a)</div><div>10000 loops, best of 3: 25.5 us per loop</div></div><div><br></div><div>-Mark</div><div><br></div><div>The documentation:</div><div>




<br></div><div><div>    einsum(subscripts, *operands, out=None, dtype=None, order=&#39;K&#39;, casting=&#39;safe&#39;)</div><div><br></div><div>    Evaluates the Einstein summation convention on the operands.</div><div><br>




</div><div>    Using the Einstein summation convention, many common multi-dimensional</div><div>    array operations can be represented in a simple fashion.  This function</div><div>    provides a way compute such summations.</div>




<div><br></div><div>    The best way to understand this function is to try the examples below,</div><div>    which show how many common NumPy functions can be implemented as</div><div>    calls to einsum.</div><div>    </div>




<div>    The subscripts string is a comma-separated list of subscript labels,</div><div>    where each label refers to a dimension of the corresponding operand.</div><div>    Repeated subscripts labels in one operand take the diagonal.  For example,</div>




<div>    ``np.einsum(&#39;ii&#39;, a)`` is equivalent to ``np.trace(a)``.</div><div>    </div><div>    Whenever a label is repeated, it is summed, so ``np.einsum(&#39;i,i&#39;, a, b)``</div><div>    is equivalent to ``np.inner(a,b)``.  If a label appears only once,</div>




<div>    it is not summed, so ``np.einsum(&#39;i&#39;, a)`` produces a view of ``a``</div><div>    with no changes.</div><div><br></div><div>    The order of labels in the output is by default alphabetical.  This</div><div>




    means that ``np.einsum(&#39;ij&#39;, a)`` doesn&#39;t affect a 2D array, while</div><div>    ``np.einsum(&#39;ji&#39;, a)`` takes its transpose.</div><div><br></div><div>    The output can be controlled by specifying output subscript labels</div>




<div>    as well.  This specifies the label order, and allows summing to be</div><div>    disallowed or forced when desired.  The call ``np.einsum(&#39;i-&gt;&#39;, a)``</div><div>    is equivalent to ``np.sum(a, axis=-1)``, and</div>




<div>    ``np.einsum(&#39;ii-&gt;i&#39;, a)`` is equivalent to ``np.diag(a)``.</div><div><br></div><div>    It is also possible to control how broadcasting occurs using</div><div>    an ellipsis.  To take the trace along the first and last axes,</div>




<div>    you can do ``np.einsum(&#39;i...i&#39;, a)``, or to do a matrix-matrix</div><div>    product with the left-most indices instead of rightmost, you can do</div><div>    ``np.einsum(&#39;ij...,jk...-&gt;ik...&#39;, a, b)``.</div>




<div><br></div><div>    When there is only one operand, no axes are summed, and no output</div><div>    parameter is provided, a view into the operand is returned instead</div><div>    of a new array.  Thus, taking the diagonal as ``np.einsum(&#39;ii-&gt;i&#39;, a)``</div>




<div>    produces a view.</div><div><br></div><div>    Parameters</div><div>    ----------</div><div>    subscripts : string</div><div>        Specifies the subscripts for summation.</div><div>    operands : list of array_like</div>




<div>        These are the arrays for the operation.</div><div>    out : None or array</div><div>        If provided, the calculation is done into this array.</div><div>    dtype : None or data type</div><div>        If provided, forces the calculation to use the data type specified.</div>




<div>        Note that you may have to also give a more liberal ``casting``</div><div>        parameter to allow the conversions.</div><div>    order : &#39;C&#39;, &#39;F&#39;, &#39;A&#39;, or &#39;K&#39;</div><div>        Controls the memory layout of the output. &#39;C&#39; means it should</div>




<div>        be Fortran contiguous. &#39;F&#39; means it should be Fortran contiguous,</div><div>        &#39;A&#39; means it should be &#39;F&#39; if the inputs are all &#39;F&#39;, &#39;C&#39; otherwise.</div><div>        &#39;K&#39; means it should be as close to the layout as the inputs as</div>




<div>        is possible, including arbitrarily permuted axes.</div><div>    casting : &#39;no&#39;, &#39;equiv&#39;, &#39;safe&#39;, &#39;same_kind&#39;, &#39;unsafe&#39;</div><div>        Controls what kind of data casting may occur.  Setting this to</div>




<div>        &#39;unsafe&#39; is not recommended, as it can adversely affect accumulations.</div><div>        &#39;no&#39; means the data types should not be cast at all. &#39;equiv&#39; means</div><div>        only byte-order changes are allowed. &#39;safe&#39; means only casts</div>




<div>        which can preserve values are allowed. &#39;same_kind&#39; means only</div><div>        safe casts or casts within a kind, like float64 to float32, are</div><div>        allowed.  &#39;unsafe&#39; means any data conversions may be done.</div>




<div><br></div><div>    Returns</div><div>    -------</div><div>    output : ndarray</div><div>        The calculation based on the Einstein summation convention.</div><div><br></div><div>    See Also</div><div>    --------</div>




<div>    dot, inner, outer, tensordot</div><div><br></div><div>    </div><div>    Examples</div><div>    --------</div><div><br></div><div>    &gt;&gt;&gt; a = np.arange(25).reshape(5,5)</div><div>    &gt;&gt;&gt; b = np.arange(5)</div>




<div>    &gt;&gt;&gt; c = np.arange(6).reshape(2,3)</div><div><br></div><div>    &gt;&gt;&gt; np.einsum(&#39;ii&#39;, a)</div><div>    60</div><div>    &gt;&gt;&gt; np.trace(a)</div><div>    60</div><div><br></div><div>    &gt;&gt;&gt; np.einsum(&#39;ii-&gt;i&#39;, a)</div>




<div>    array([ 0,  6, 12, 18, 24])</div><div>    &gt;&gt;&gt; np.diag(a)</div><div>    array([ 0,  6, 12, 18, 24])</div><div><br></div><div>    &gt;&gt;&gt; np.einsum(&#39;ij,j&#39;, a, b)</div><div>    array([ 30,  80, 130, 180, 230])</div>




<div>    &gt;&gt;&gt; np.dot(a, b)</div><div>    array([ 30,  80, 130, 180, 230])</div><div><br></div><div>    &gt;&gt;&gt; np.einsum(&#39;ji&#39;, c)</div><div>    array([[0, 3],</div><div>           [1, 4],</div><div>           [2, 5]])</div>




<div>    &gt;&gt;&gt; c.T</div><div>    array([[0, 3],</div><div>           [1, 4],</div><div>           [2, 5]])</div><div><br></div><div>    &gt;&gt;&gt; np.einsum(&#39;,&#39;, 3, c)</div><div>    array([[ 0,  3,  6],</div>




<div>           [ 9, 12, 15]])</div><div>    &gt;&gt;&gt; np.multiply(3, c)</div><div>    array([[ 0,  3,  6],</div><div>           [ 9, 12, 15]])</div><div><br></div><div>    &gt;&gt;&gt; np.einsum(&#39;i,i&#39;, b, b)</div>




<div>    30</div><div>    &gt;&gt;&gt; np.inner(b,b)</div><div>    30</div><div><br></div><div>    &gt;&gt;&gt; np.einsum(&#39;i,j&#39;, np.arange(2)+1, b)</div><div>    array([[0, 1, 2, 3, 4],</div><div>           [0, 2, 4, 6, 8]])</div>




<div>    &gt;&gt;&gt; np.outer(np.arange(2)+1, b)</div><div>    array([[0, 1, 2, 3, 4],</div><div>           [0, 2, 4, 6, 8]])</div><div><br></div><div>    &gt;&gt;&gt; np.einsum(&#39;i...-&gt;&#39;, a)</div><div>    array([50, 55, 60, 65, 70])</div>




<div>    &gt;&gt;&gt; np.sum(a, axis=0)</div><div>    array([50, 55, 60, 65, 70])</div><div><br></div><div>    &gt;&gt;&gt; a = np.arange(60.).reshape(3,4,5)</div><div>    &gt;&gt;&gt; b = np.arange(24.).reshape(4,3,2)</div>




<div>    &gt;&gt;&gt; np.einsum(&#39;ijk,jil-&gt;kl&#39;, a, b)</div><div>    array([[ 4400.,  4730.],</div><div>           [ 4532.,  4874.],</div><div>           [ 4664.,  5018.],</div><div>           [ 4796.,  5162.],</div>




<div>           [ 4928.,  5306.]])</div><div>    &gt;&gt;&gt; np.tensordot(a,b, axes=([1,0],[0,1]))</div><div>    array([[ 4400.,  4730.],</div><div>           [ 4532.,  4874.],</div><div>           [ 4664.,  5018.],</div>




<div>           [ 4796.,  5162.],</div><div>           [ 4928.,  5306.]])</div></div><div><br></div></div>
<br>_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org" target="_blank">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
<br></blockquote></div><br><br clear="all"><br>-- <br>Jonathan Rocher,<br><font color="#888888"><span>Enthought</span>, Inc.<br>
<a href="mailto:jrocher@enthought.com" target="_blank">jrocher@<span>enthought</span>.com</a><br>
1-512-536-1057<br>
<a href="http://www.enthought.com/" target="_blank">http://www.<span>enthought</span>.com</a></font><br><div style="display: inline;"></div><br>