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