Re: [Numpy-discussion] algorithm, optimization, or other problem?

[Bruce S. <bso...@gm...>]

Hi,
Actually it makes it slightly worse - given the responses on another
thread it is probably due to not pushing enough into C code.

Obviously use of blas etc will be faster but it doesn't change the
fact that removing the inner loop would be faster still.

Bruce

On 2/22/06, Nadav Horesh <na...@vi...> wrote:
> You may get a significant boost by replacing the line:
>   w=3Dw+ eta * (y*x - y**2*w)
> with
>   w *=3D 1.0 - eta*y*y
>   w +=3D eta*y*x
>
> I ran a test on a similar expression and got 5 fold speed increase.
> The dot() function runs faster if you compile with dotblas.
>
>   Nadav.
>
>
> -----Original Message-----
> From:   num...@li... on behalf of Bruce S=
outhey
> Sent:   Tue 21-Feb-06 17:15
> To:     Brian Blais
> Cc:     pyt...@py...; num...@li...; s=
cip...@sc...
> Subject:        Re: [Numpy-discussion] algorithm, optimization, or other =
problem?
> Hi,
> In the current version, note that Y is scalar so replace the squaring
> (Y**2) with Y*Y as you do in the dohebb function.  On my system
> without blas etc removing the squaring removes a few seconds (16.28 to
> 12.4). It did not seem to help factorizing Y.
>
> Also, eta and tau are constants so define them only once as scalars
> outside the loops and do the division outside the loop. It only saves
> about 0.2 seconds but these add up.
>
> The inner loop probably can be vectorized because it is just vector
> operations on a matrix. You are just computing over the ith dimension
> of X.  I think that you could be able to find the matrix version on
> the net.
>
> Regards
> Bruce
>
>
>
> On 2/21/06, Brian Blais <bb...@br...> wrote:
> > Hello,
> >
> > I am trying to translate some Matlab/mex code to Python, for doing neur=
al
> > simulations.  This application is definitely computing-time limited, an=
d I need to
> > optimize at least one inner loop of the code, or perhaps even rethink t=
he algorithm.
> >    The procedure is very simple, after initializing any variables:
> >
> > 1) select a random input vector, which I will call "x".  right now I ha=
ve it as an
> > array, and I choose columns from that array randomly.  in other cases, =
I may need to
> > take an image, select a patch, and then make that a column vector.
> >
> > 2) calculate an output value, which is the dot product of the "x" and a=
 weight
> > vector, "w", so
> >
> >         y=3Ddot(x,w)
> >
> > 3) modify the weight vector based on a matrix equation, like:
> >
> >         w=3Dw+ eta * (y*x - y**2*w)
> >                ^
> >                |
> >                +---- learning rate constant
> >
> > 4) repeat steps 1-3 many times
> >
> > I've organized it like:
> >
> > for e in 100:   # outer loop
> >      for i in 1000:  # inner loop
> >          (steps 1-3)
> >
> >      display things.
> >
> > so that the bulk of the computation is in the inner loop, and is amenab=
le to
> > converting to a faster language.  This is my issue:
> >
> > straight python, in the example posted below for 250000 inner-loop step=
s, takes 20
> > seconds for each outer-loop step.  I tried Pyrex, which should work ver=
y fast on such
> > a problem, takes about 8.5 seconds per outer-loop step.  The same code =
as a C-mex
> > file in matlab takes 1.5 seconds per outer-loop step.
> >
> > Given the huge difference between the Pyrex and the Mex, I feel that th=
ere is
> > something I am doing wrong, because the C-code for both should run comp=
arably.
> > Perhaps the approach is wrong?  I'm willing to take any suggestions!  I=
 don't mind
> > coding some in C, but the Python API seemed a bit challenging to me.
> >
> > One note: I am using the Numeric package, not numpy, only because I wan=
t to be able
> > to use the Enthought version for Windows.  I develop on Linux, and have=
n't had a
> > chance to see if I can compile numpy using the Enthought Python for Win=
dows.
> >
> > If there is anything else anyone needs to know, I'll post it.  I put th=
e main script,
> > and a dohebb.pyx code below.
> >
> >
> >                         thanks!
> >
> >                                         Brian Blais
> >
> > --
> > -----------------
> >
> >              bb...@br...
> >              http://web.bryant.edu/~bblais
> >
> >
> >
> >
> > # Main script:
> >
> > from dohebb import *
> > import pylab as p
> > from Numeric import *
> > from RandomArray import *
> > import time
> >
> > x=3Drandom((100,1000))    # 1000 input vectors
> >
> > numpats=3Dx.shape[0]
> > w=3Drandom((numpats,1));
> >
> > th=3Drandom((1,1))
> >
> > params=3D{}
> > params['eta']=3D0.001;
> > params['tau']=3D100.0;
> > old_mx=3D0;
> > for e in range(100):
> >
> >      rnd=3Drandint(0,numpats,250000)
> >      t1=3Dtime.time()
> >      if 0:  # straight python
> >          for i in range(len(rnd)):
> >              pat=3Drnd[i]
> >              xx=3Dreshape(x[:,pat],(1,-1))
> >              y=3Dmatrixmultiply(xx,w)
> >              w=3Dw+params['eta']*(y*transpose(xx)-y**2*w);
> >              th=3Dth+(1.0/params['tau'])*(y**2-th);
> >      else: # pyrex
> >          dohebb(params,w,th,x,rnd)
> >      print time.time()-t1
> >
> >
> > p.plot(w,'o-')
> > p.xlabel('weights')
> > p.show()
> >
> >
> > #=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
> >
> > #         dohebb.pyx
> >
> > cdef extern from "Numeric/arrayobject.h":
> >
> >    struct PyArray_Descr:
> >      int type_num, elsize
> >      char type
> >
> >    ctypedef class Numeric.ArrayType [object PyArrayObject]:
> >      cdef char *data
> >      cdef int nd
> >      cdef int *dimensions, *strides
> >      cdef object base
> >      cdef PyArray_Descr *descr
> >      cdef int flags
> >
> >
> > def dohebb(params,ArrayType w,ArrayType th,ArrayType X,ArrayType rnd):
> >
> >
> >      cdef int num_iterations
> >      cdef int num_inputs
> >      cdef int offset
> >      cdef double *wp,*xp,*thp
> >      cdef int *rndp
> >      cdef double eta,tau
> >
> >      eta=3Dparams['eta']  # learning rate
> >      tau=3Dparams['tau']  # used for variance estimate
> >
> >      cdef double y
> >      num_iterations=3Drnd.dimensions[0]
> >      num_inputs=3Dw.dimensions[0]
> >
> >      # get the pointers
> >      wp=3D<double *>w.data
> >      xp=3D<double *>X.data
> >      rndp=3D<int *>rnd.data
> >      thp=3D<double *>th.data
> >
> >      for it from 0 <=3D it < num_iterations:
> >
> >          offset=3Drndp[it]*num_inputs
> >
> >          # calculate the output
> >          y=3D0.0
> >          for i from 0 <=3D i < num_inputs:
> >              y=3Dy+wp[i]*xp[i+offset]
> >
> >          # change in the weights
> >          for i from 0 <=3D i < num_inputs:
> >              wp[i]=3Dwp[i]+eta*(y*xp[i+offset] - y*y*wp[i])
> >
> >          # estimate the variance
> >          thp[0]=3Dthp[0]+(1.0/tau)*(y**2-thp[0])
> >
> >
> >
> >
> >
> >
> >
> > -------------------------------------------------------
> > This SF.net email is sponsored by: Splunk Inc. Do you grep through log =
files
> > for problems?  Stop!  Download the new AJAX search engine that makes
> > searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
> > http://sel.as-us.falkag.net/sel?cmd=3Dlnk&kid=3D103432&bid=3D230486&dat=
=3D121642
> > _______________________________________________
> > Numpy-discussion mailing list
> > Num...@li...
> > https://lists.sourceforge.net/lists/listinfo/numpy-discussion
> >
>
>
> -------------------------------------------------------
> This SF.net email is sponsored by: Splunk Inc. Do you grep through log fi=
les
> for problems?  Stop!  Download the new AJAX search engine that makes
> searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
> http://sel.as-us.falkag.net/sel?cmd=3Dk&kid=103432&bid#0486&dat=121642
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>
>
>