[Numpy-discussion] Trick for fast
josef.pktd at gmail.com
josef.pktd at gmail.com
Fri Feb 3 13:51:27 EST 2012
On Fri, Feb 3, 2012 at 1:29 PM, santhu kumar <mesanthu at 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 at scipy.org> wrote:
>>
>> Send NumPy-Discussion mailing list submissions to
>> numpy-discussion at 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 at scipy.org
>>
>> You can reach the person managing the list at
>> numpy-discussion-owner at 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 at 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 at gmail.com
>> Subject: Re: [Numpy-discussion] Trick for fast
>> To: Discussion of Numerical Python <numpy-discussion at scipy.org>
>> Message-ID:
>>
>> <CAMMTP+B2yQ9Th5yhDEuogS6yWjuQpAnZZQ1GcU1GqRdvDk7itw at 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 at 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 at 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 at gmail.com>
>> Subject: Re: [Numpy-discussion] Trick for fast
>> To: Discussion of Numerical Python <numpy-discussion at scipy.org>
>> Message-ID:
>>
>> <CAJO1x6qV0ME_2Xgr4oU8q=waekT2qXQ=iComzbbQNGzB_MPKkQ at 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 at gmail.com> wrote:
>>
>> > On Fri, Feb 3, 2012 at 8:44 AM, Alan G Isaac <alan.isaac at 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 at scipy.org
>> > > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> > _______________________________________________
>> > NumPy-Discussion mailing list
>> > NumPy-Discussion at 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 at sipsolutions.net>
>> Subject: Re: [Numpy-discussion] Trick for fast
>> To: Discussion of Numerical Python <numpy-discussion at scipy.org>
>> Message-ID: <1328282044.2830.9.camel at 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 at gmail.com> wrote:
>> > On Fri, Feb 3, 2012 at 8:44 AM, Alan G Isaac
>> > <alan.isaac at 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 at scipy.org
>> > > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> > _______________________________________________
>> > NumPy-Discussion mailing list
>> > NumPy-Discussion at scipy.org
>> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> >
>> >
>> >
>> > _______________________________________________
>> > NumPy-Discussion mailing list
>> > NumPy-Discussion at scipy.org
>> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
>>
>>
>> ------------------------------
>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at 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 at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
More information about the NumPy-Discussion
mailing list