[Numpy-discussion] Help in speeding up accumulation in a matrix

[josef.pktd at gmail.com]

On Sat, Jan 29, 2011 at 5:30 PM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
>
>
> On Sat, Jan 29, 2011 at 2:01 PM, Nicolas SCHEFFER
> <scheffer.nicolas at gmail.com> wrote:
>>
>> Hi all,
>>
>> First email to the list for me, I just want to say how grateful I am
>> to have python+numpy+ipython etc... for my day to day needs. Great
>> combination of software.
>>
>> Anyway, I've been having this bottleneck in one my algorithms that has
>> been bugging me for quite a while.
>> The objective is to speed this part up. I've been doing tons of
>> optimization and parallel processing around that piece of code to get
>> a decent run time.
>>
>> The problem is easy. You want to accumulate in a matrix, a weighted
>> sum of other matrices. Let's call this function scale_and_add:
>> def scale_and_add_re(R,w,Ms):
>>    (nb_add,mdim,j)=np.shape(Ms)
>>    for i in range(nb_add):
>>        R+=w[i]*Ms[i]
>>    return R
>> This 'for' loop bugs me since I know this will slow things down.
>>
>
> I'd put the flattened matrices in a stack, weight in place, sum on the first
> index, and reshape. Something like
>
> In [1]: m = array([eye(3)]*4).reshape(4,-1)

using triu_indices instead of reshape would cut the memory consumption

ind1, ind2 = np.triu_indices(?,?)
m = m[:, ind1, ind2]

I have no idea how expensive the indexing is in this case

 Josef

>
> In [2]: m
> Out[2]:
> array([[ 1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.],
>        [ 1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.],
>        [ 1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.],
>        [ 1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.]])
>
> In [3]: w = array([1.,2.,3.,4.])
>
> In [4]: m *= w[:,None]
>
> In [5]: r = m.sum(0).reshape(3,3)
>
> In [6]: r
> Out[6]:
> array([[ 10.,   0.,   0.],
>        [  0.,  10.,   0.],
>        [  0.,   0.,  10.]])
>
> This should fit in memory I think, depending of course on how much memory
> you have.
>
> <snip>
>
> Chuck
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>