[Numpy-discussion] Trick for fast

josef.pktd@gmai... josef.pktd@gmai...
Fri Feb 3 12:51:27 CST 2012


On Fri, Feb 3, 2012 at 1:29 PM, santhu kumar <mesanthu@gmail.com> wrote:
> Hello all,
>
> Thanks for lovely solutions. I have sat on it for some time and wrote it
> myself :
>
> n =x.shape[0]
> ea = np.array([1,0,0,0,1,0,0,0,1])
> inert = ((np.tile(ea,(n,1))*((x*x).sum(axis=1)[:,np.newaxis]) -
> np.hstack([x*x[:,0][:,np.newaxis],x*x[:,1][:,np.newaxis],x*x[:,2][:,np.newaxis]]))*mass[:,np.newaxis]).sum(axis=0)
> inert.shape = 3,3
>
> Does the trick and reduces the time from over 45 secs to 3 secs.
> I do want to try einsum but my numpy is little old and it does not have it.
>
> Thanks Sebastian (it was tricky to understand your code for me) and Josef
> (clean).

Isn't the entire substraction of the first term just to set the
diagonal of the result to zero.

It looks to me now just like the weighted dot product and setting the
diagonal to zero. That shouldn't take 3 secs unless you actual
dimensions are huge.

Josef




>
>
> On Fri, Feb 3, 2012 at 12:00 PM, <numpy-discussion-request@scipy.org> wrote:
>>
>> Send NumPy-Discussion mailing list submissions to
>>        numpy-discussion@scipy.org
>>
>> To subscribe or unsubscribe via the World Wide Web, visit
>>        http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> or, via email, send a message with subject or body 'help' to
>>        numpy-discussion-request@scipy.org
>>
>> You can reach the person managing the list at
>>        numpy-discussion-owner@scipy.org
>>
>> When replying, please edit your Subject line so it is more specific
>> than "Re: Contents of NumPy-Discussion digest..."
>>
>>
>> Today's Topics:
>>
>>   1. Re: Trick for fast (josef.pktd@gmail.com)
>>   2. Re: Trick for fast (S?ren Gammelmark)
>>   3. Re: Trick for fast (Sebastian Berg)
>>
>>
>> ----------------------------------------------------------------------
>>
>> Message: 1
>> Date: Fri, 3 Feb 2012 09:10:28 -0500
>> From: josef.pktd@gmail.com
>> Subject: Re: [Numpy-discussion] Trick for fast
>> To: Discussion of Numerical Python <numpy-discussion@scipy.org>
>> Message-ID:
>>
>>  <CAMMTP+B2yQ9Th5yhDEuogS6yWjuQpAnZZQ1GcU1GqRdvDk7itw@mail.gmail.com>
>> Content-Type: text/plain; charset=ISO-8859-1
>>
>>
>> On Fri, Feb 3, 2012 at 8:44 AM, Alan G Isaac <alan.isaac@gmail.com> wrote:
>> > On 2/3/2012 5:16 AM, santhu kumar wrote:
>> >> x = nX3 vector.
>> >> mass = nX1 vector
>> >> inert = zeros((3,3))
>> >> for i in range(n):
>> >> ? ? ? ?ri = x[i,:].reshape(1,3)
>> >> ? ? ? ?inert = inert + mass[i,]*(sum(ri*ri)*eye(3) - dot(ri.T,ri))
>>
>> >>
>> >
>> >
>> > This should buy you a bit.
>> >
>> > xdot = (x*x).sum(axis=1)
>> > for (massi,xi,xdoti) in zip(mass.flat,x,xdot):
>> > ? ? ? temp = -np.outer(xi,xi)
>> > ? ? ? temp.flat[slice(0,None,4)] += xdoti
>> > ? ? ? inert += massi*temp
>>
>> >
>> > Alan Isaac
>>
>> maybe something like this, (self contained example and name spaces to
>> make running it easier)
>>
>> import numpy as np
>> n = 15
>> x = np.arange(n*3.).reshape(-1,3) #nX3 vector.
>> mass = np.linspace(1,2,n)[:,None] #nX1 vector
>> inert = np.zeros((3,3))
>> for i in range(n):
>>      ri = x[i,:].reshape(1,3)
>>      inert = inert + mass[i,]*(sum(ri*ri)*np.eye(3) - np.dot(ri.T,ri))
>> print inert
>>
>> print  np.diag((mass * x**2).sum(0))  - np.dot(x.T, mass*x)
>>
>> [[     0.  -16755.  -17287.5]
>>  [-16755.       0.  -17865. ]
>>  [-17287.5 -17865.       0. ]]
>> [[     0.  -16755.  -17287.5]
>>  [-16755.       0.  -17865. ]
>>  [-17287.5 -17865.       0. ]]
>>
>> Josef
>>
>>
>> >
>> > _______________________________________________
>> > NumPy-Discussion mailing list
>> > NumPy-Discussion@scipy.org
>> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
>> ------------------------------
>>
>> Message: 2
>> Date: Fri, 3 Feb 2012 15:43:25 +0100
>> From: S?ren Gammelmark <gammelmark@gmail.com>
>> Subject: Re: [Numpy-discussion] Trick for fast
>> To: Discussion of Numerical Python <numpy-discussion@scipy.org>
>> Message-ID:
>>
>>  <CAJO1x6qV0ME_2Xgr4oU8q=waekT2qXQ=iComzbbQNGzB_MPKkQ@mail.gmail.com>
>> Content-Type: text/plain; charset="iso-8859-1"
>>
>>
>> What about this?
>>
>> A = einsum("i,ij->", mass, x ** 2)
>> B = einsum("i,ij,ik->jk", mass, x, x)
>> I = A * eye(3) - B
>>
>> /S?ren
>>
>>
>> On 3 February 2012 15:10, <josef.pktd@gmail.com> wrote:
>>
>> > On Fri, Feb 3, 2012 at 8:44 AM, Alan G Isaac <alan.isaac@gmail.com>
>> > wrote:
>> > > On 2/3/2012 5:16 AM, santhu kumar wrote:
>> > >> x = nX3 vector.
>> > >> mass = nX1 vector
>> > >> inert = zeros((3,3))
>> > >> for i in range(n):
>> > >>        ri = x[i,:].reshape(1,3)
>> > >>        inert = inert + mass[i,]*(sum(ri*ri)*eye(3) - dot(ri.T,ri))
>> > >>
>> > >
>> > >
>> > > This should buy you a bit.
>> > >
>> > > xdot = (x*x).sum(axis=1)
>> > > for (massi,xi,xdoti) in zip(mass.flat,x,xdot):
>> > >       temp = -np.outer(xi,xi)
>> > >       temp.flat[slice(0,None,4)] += xdoti
>> > >       inert += massi*temp
>> > >
>> > > Alan Isaac
>> >
>> > maybe something like this, (self contained example and name spaces to
>> > make running it easier)
>> >
>> > import numpy as np
>> > n = 15
>> > x = np.arange(n*3.).reshape(-1,3) #nX3 vector.
>> > mass = np.linspace(1,2,n)[:,None] #nX1 vector
>> > inert = np.zeros((3,3))
>> > for i in range(n):
>> >      ri = x[i,:].reshape(1,3)
>> >       inert = inert + mass[i,]*(sum(ri*ri)*np.eye(3) - np.dot(ri.T,ri))
>> > print inert
>> >
>> > print  np.diag((mass * x**2).sum(0))  - np.dot(x.T, mass*x)
>> >
>> > [[     0.  -16755.  -17287.5]
>> >  [-16755.       0.  -17865. ]
>> >  [-17287.5 -17865.       0. ]]
>> > [[     0.  -16755.  -17287.5]
>> >  [-16755.       0.  -17865. ]
>> >  [-17287.5 -17865.       0. ]]
>> >
>> > Josef
>> >
>> >
>> > >
>> > > _______________________________________________
>> > > NumPy-Discussion mailing list
>> > > NumPy-Discussion@scipy.org
>> > > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> > _______________________________________________
>> > NumPy-Discussion mailing list
>> > NumPy-Discussion@scipy.org
>> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> >
>> -------------- next part --------------
>> An HTML attachment was scrubbed...
>> URL:
>> http://mail.scipy.org/pipermail/numpy-discussion/attachments/20120203/d1faa546/attachment-0001.html
>>
>> ------------------------------
>>
>> Message: 3
>> Date: Fri, 03 Feb 2012 16:14:04 +0100
>> From: Sebastian Berg <sebastian@sipsolutions.net>
>> Subject: Re: [Numpy-discussion] Trick for fast
>> To: Discussion of Numerical Python <numpy-discussion@scipy.org>
>> Message-ID: <1328282044.2830.9.camel@sebastian-laptop>
>> Content-Type: text/plain; charset="UTF-8"
>>
>>
>> I guess Einsum is much cleaner, but I already had started with this and
>> maybe someone likes it, this is fully vectorized and uses a bit of funny
>> stuff too:
>>
>> # The dot product(s), written using broadcasting rules:
>> a = -(x.reshape(-1,1,3) * x[...,None])
>>
>> # Magic, to avoid the eye thing, takes all diagonal elements as view,
>> maybe there is a cooler way for it:
>> diagonals = np.lib.stride_tricks.as_strided(a, (a.shape[0], 3),
>> (a.dtype.itemsize*9, a.dtype.itemsize*4))
>>
>> # Add the x**2 (s is a view on the diagonals), the sum is broadcasted.
>> diagonals += (sum(x**2, 1))[:,None]
>>
>> # And multiply by mass using broadcasting:
>> a *= mass[...,None]
>>
>> # And sum up all the intermediat results:
>> inert = a.sum(0)
>>
>> print inert
>>
>> Regards,
>>
>> Sebastian
>>
>> On Fri, 2012-02-03 at 15:43 +0100, S?ren Gammelmark wrote:
>> > What about this?
>> >
>> >
>> > A = einsum("i,ij->", mass, x ** 2)
>> > B = einsum("i,ij,ik->jk", mass, x, x)
>> > I = A * eye(3) - B
>> >
>> >
>> > /S?ren
>>
>> >
>> > On 3 February 2012 15:10, <josef.pktd@gmail.com> wrote:
>> >         On Fri, Feb 3, 2012 at 8:44 AM, Alan G Isaac
>> >         <alan.isaac@gmail.com> wrote:
>> >         > On 2/3/2012 5:16 AM, santhu kumar wrote:
>> >         >> x = nX3 vector.
>> >         >> mass = nX1 vector
>> >         >> inert = zeros((3,3))
>> >         >> for i in range(n):
>> >         >>        ri = x[i,:].reshape(1,3)
>> >         >>        inert = inert + mass[i,]*(sum(ri*ri)*eye(3) -
>> >         dot(ri.T,ri))
>> >         >>
>> >         >
>> >         >
>> >         > This should buy you a bit.
>> >         >
>> >         > xdot = (x*x).sum(axis=1)
>> >         > for (massi,xi,xdoti) in zip(mass.flat,x,xdot):
>> >         >       temp = -np.outer(xi,xi)
>> >         >       temp.flat[slice(0,None,4)] += xdoti
>> >         >       inert += massi*temp
>> >         >
>> >         > Alan Isaac
>> >
>> >
>> >         maybe something like this, (self contained example and name
>> >         spaces to
>> >         make running it easier)
>> >
>> >         import numpy as np
>> >         n = 15
>> >         x = np.arange(n*3.).reshape(-1,3) #nX3 vector.
>> >         mass = np.linspace(1,2,n)[:,None] #nX1 vector
>> >         inert = np.zeros((3,3))
>> >         for i in range(n):
>> >              ri = x[i,:].reshape(1,3)
>> >
>> >              inert = inert + mass[i,]*(sum(ri*ri)*np.eye(3) -
>> >         np.dot(ri.T,ri))
>> >         print inert
>> >
>> >         print  np.diag((mass * x**2).sum(0))  - np.dot(x.T, mass*x)
>> >
>> >         [[     0.  -16755.  -17287.5]
>> >          [-16755.       0.  -17865. ]
>> >          [-17287.5 -17865.       0. ]]
>> >         [[     0.  -16755.  -17287.5]
>> >          [-16755.       0.  -17865. ]
>> >          [-17287.5 -17865.       0. ]]
>> >
>> >         Josef
>> >
>> >
>> >         >
>> >         > _______________________________________________
>> >         > NumPy-Discussion mailing list
>> >         > NumPy-Discussion@scipy.org
>> >         > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> >         _______________________________________________
>> >         NumPy-Discussion mailing list
>> >         NumPy-Discussion@scipy.org
>> >         http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> >
>> >
>> >
>> > _______________________________________________
>> > NumPy-Discussion mailing list
>> > NumPy-Discussion@scipy.org
>> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
>>
>>
>> ------------------------------
>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
>> End of NumPy-Discussion Digest, Vol 65, Issue 11
>> ************************************************
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>


More information about the NumPy-Discussion mailing list