[Numpy-discussion] identity

[Lars Bittrich]

Hi,

a colleague made me aware of a speed issue with numpy.identity. Since he was 
using numpy.diag(numpy.ones(N)) before, he expected identity to be at least as 
fast as diag. But that is not the case.

We found that there was a discussion on the list (July, 20th; "My identity" by 
Keith Goodman). The presented solution was much faster. Someone wondered if 
the change was already made in the svn. 
But I got something different:

In [1]:import numpy

In [2]:numpy.__version__
Out[2]:'1.4.0.dev7301'

In [3]:numpy.identity??
[...]
def identity(n, dtype=None):
    """
    [...]
    """
    a = array([1]+n*[0],dtype=dtype)
    b = empty((n,n),dtype=dtype)

    # Note that this assignment depends on the convention that since the a
    # array is shorter than the flattened b array, then the a array will
    # be repeated until it is the appropriate size. Given a's construction,
    # this nicely sets the diagonal to all ones.
    b.flat = a
    return b

instead of (mail by Keith Goodman):

def myidentity(n, dtype=None):
    a = zeros((n,n), dtype=dtype)
    a.flat[::n+1] = 1
    return a


Did I look at the wrong place or is there a reason to keep the slow version of 
identity?

Lars