[Numpy-discussion] How to improve performance of slow tri*_indices calculations?

[josef.pktd at gmail.com]

On Mon, Jan 24, 2011 at 6:49 PM,  <josef.pktd at gmail.com> wrote:
> On Mon, Jan 24, 2011 at 4:29 PM, eat <e.antero.tammi at gmail.com> wrote:
>> Hi,
>>
>> Running on:
>> In []: np.__version__
>> Out[]: '1.5.1'
>> In []: sys.version
>> Out[]: '2.7.1 (r271:86832, Nov 27 2010, 18:30:46) [MSC v.1500 32 bit (Intel)]'
>>
>> For the reference:
>> In []: X= randn(10, 125)
>> In []: timeit dot(X.T, X)
>> 10000 loops, best of 3: 170 us per loop
>> In []: X= randn(10, 250)
>> In []: timeit dot(X.T, X)
>> 1000 loops, best of 3: 671 us per loop
>> In []: X= randn(10, 500)
>> In []: timeit dot(X.T, X)
>> 100 loops, best of 3: 5.15 ms per loop
>> In []: X= randn(10, 1000)
>> In []: timeit dot(X.T, X)
>> 100 loops, best of 3: 20 ms per loop
>> In []: X= randn(10, 2000)
>> In []: timeit dot(X.T, X)
>> 10 loops, best of 3: 80.7 ms per loop
>>
>> Performance of triu_indices:
>> In []: timeit triu_indices(125)
>> 1000 loops, best of 3: 662 us per loop
>> In []: timeit triu_indices(250)
>> 100 loops, best of 3: 2.55 ms per loop
>> In []: timeit triu_indices(500)
>> 100 loops, best of 3: 15 ms per loop
>> In []: timeit triu_indices(1000)
>> 10 loops, best of 3: 59.8 ms per loop
>> In []: timeit triu_indices(2000)
>> 1 loops, best of 3: 239 ms per loop
>>
>> So the tri*_indices calculations seems to be unreasonable slow compared to for
>> example calculations of inner products.
>>
>> Now, just to compare for a very naive implementation of triu indices.
>> In []: def iut(n):
>>   ..:     r= np.empty(n* (n+ 1)/ 2, dtype= int)
>>   ..:     c= r.copy()
>>   ..:     a= np.arange(n)
>>   ..:     m= 0
>>   ..:     for i in xrange(n):
>>   ..:         ni= n- i
>>   ..:         mni= m+ ni
>>   ..:         r[m: mni]= i
>>   ..:         c[m: mni]= a[i: n]
>>   ..:         m+= ni
>>   ..:     return (r, c)
>>   ..:
>>
>> Are we really calculating the same thing?
>> In []: triu_indices(5)
>> Out[]:
>> (array([0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]),
>>  array([0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]))
>> In []: iut(5)
>> Out[]:
>> (array([0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]),
>>  array([0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]))
>>
>> Seems so, and then its performance:
>> In []: timeit iut(125)
>> 1000 loops, best of 3: 992 us per loop
>> In []: timeit iut(250)
>> 100 loops, best of 3: 2.03 ms per loop
>> In []: timeit iut(500)
>> 100 loops, best of 3: 5.3 ms per loop
>> In []: timeit iut(1000)
>> 100 loops, best of 3: 13.9 ms per loop
>> In []: timeit iut(2000)
>> 10 loops, best of 3: 39.8 ms per loop
>>
>> Even the naive implementation is very slow, but allready outperforms
>> triu_indices, when n is > 250!
>>
>> So finally my question is how one could substantially improve the performance
>> of indices calculations?
>
> What's the timing of this version (taken from nitime) ? it builds a
> full intermediate array

I should have checked the numpy source first, that's exactly the
implementation of triu_indices in numpy 1.5.1

Josef

>
>    m = np.ones((n,n),int)
>    a = np.triu(m,k)
>    np.where(a != 0)
>
>>>> n=5
>>>> m = np.ones((n,n),int)
>>>> np.where(np.triu(m,0) != 0)
> (array([0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]),
> array([0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]))
>
> or maybe variations on it that all build a full intermediate matrix
>
>>>> np.where(1-np.tri(n,n,-1))
> (array([0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]),
> array([0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]))
>
>>>> np.where(np.subtract.outer(np.arange(n), np.arange(n))<=0)
> (array([0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]),
> array([0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]))
>
> Josef
>
>>
>>
>> Regards,
>> eat
>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>