# [Numpy-discussion] Trick for fast

santhu kumar mesanthu@gmail....
Fri Feb 3 12:29:26 CST 2012

```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).

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
>
> 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"
>
>
> 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:
> >
> >
> > 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
> ************************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20120203/4b3ae03d/attachment.html
```