<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
  <meta content="text/html;charset=ISO-8859-1" http-equiv="Content-Type">
  <title></title>
</head>
<body bgcolor="#ffffff" text="#000000">
It is slower.<br>
<br>
I did a little study on this issue since I got into the issue of
algorithms that can not be easily vectorized (like this one).<br>
On my PC an outer loop step took initially 17.3 seconds, and some
optimization brought it down to ~11 seconds.&nbsp; The dot product consumed
about 1/3 of the time. I estimate that objects creation/destruction
consumes most of the cpu time. It seems that this way comes nowhere
near cmex speed. I suspect that maybe blitz/boost may bridge the gap.<br>
<br>
My optimization (did not verify numeric results)<br>
<br>
from Numeric import *<br>
from RandomArray import *<br>
import time<br>
<br>
<br>
def run():<br>
<br>
&nbsp;&nbsp;&nbsp; x=random((100,1000))&nbsp;&nbsp;&nbsp; # 1000 input vectors<br>
<br>
&nbsp;&nbsp;&nbsp; numpats=x.shape[0]<br>
&nbsp;&nbsp;&nbsp; w=random(numpats)<br>
&nbsp;&nbsp;&nbsp; th=random((1,1))<br>
&nbsp;&nbsp;&nbsp; eta=array(0.001)<br>
&nbsp;&nbsp;&nbsp; tau = array(0.01)<br>
&nbsp;&nbsp;&nbsp; old_mx=0;<br>
&nbsp;&nbsp;&nbsp; for e in range(1):<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; rnd=randint(0,numpats,250000)<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; t1=time.time()<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if 1:&nbsp; # straight python<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; for i in xrange(len(rnd)):<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; pat=rnd[i]<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; xx = x[:,pat]<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; y = array(dot(xx,w))<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; w += eta*y*(xx-y*w)<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; th += tau * (y*y-th)<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; else: # pyrex<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dohebb(params,w,th,x,rnd)<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; print time.time()-t1<br>
<br>
if __name__ == '__main__':<br>
&nbsp;&nbsp;&nbsp; run()<br>
<br>
<br>
<br>
Bruce Southey wrote:
<blockquote
 cite="midbbcd77d00602220624m3ef00a75i3a93b73dc58fd0b6@mail.gmail.com"
 type="cite">
  <pre wrap="">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 <a class="moz-txt-link-rfc2396E" href="mailto:nadavh@visionsense.com">&lt;nadavh@visionsense.com&gt;</a> wrote:
  </pre>
  <blockquote type="cite">
    <pre wrap="">You may get a significant boost by replacing the line:
  w=w+ eta * (y*x - y**2*w)
with
  w *= 1.0 - eta*y*y
  w += 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:   <a class="moz-txt-link-abbreviated" href="mailto:numpy-discussion-admin@lists.sourceforge.net">numpy-discussion-admin@lists.sourceforge.net</a> on behalf of Bruce Southey
Sent:   Tue 21-Feb-06 17:15
To:     Brian Blais
Cc:     <a class="moz-txt-link-abbreviated" href="mailto:python-list@python.org">python-list@python.org</a>; <a class="moz-txt-link-abbreviated" href="mailto:numpy-discussion@lists.sourceforge.net">numpy-discussion@lists.sourceforge.net</a>; <a class="moz-txt-link-abbreviated" href="mailto:scipy-user@scipy.net">scipy-user@scipy.net</a>
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 <a class="moz-txt-link-rfc2396E" href="mailto:bblais@bryant.edu">&lt;bblais@bryant.edu&gt;</a> wrote:
    </pre>
    <blockquote type="cite">
      <pre wrap="">Hello,

I am trying to translate some Matlab/mex code to Python, for doing neural
simulations.  This application is definitely computing-time limited, and I need to
optimize at least one inner loop of the code, or perhaps even rethink the algorithm.
   The procedure is very simple, after initializing any variables:

1) select a random input vector, which I will call "x".  right now I have 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=dot(x,w)

3) modify the weight vector based on a matrix equation, like:

        w=w+ 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 amenable to
converting to a faster language.  This is my issue:

straight python, in the example posted below for 250000 inner-loop steps, takes 20
seconds for each outer-loop step.  I tried Pyrex, which should work very 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 there is
something I am doing wrong, because the C-code for both should run comparably.
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 want to be able
to use the Enthought version for Windows.  I develop on Linux, and haven't had a
chance to see if I can compile numpy using the Enthought Python for Windows.

If there is anything else anyone needs to know, I'll post it.  I put the main script,
and a dohebb.pyx code below.


                        thanks!

                                        Brian Blais

--
-----------------

             <a class="moz-txt-link-abbreviated" href="mailto:bblais@bryant.edu">bblais@bryant.edu</a>
             <a class="moz-txt-link-freetext" href="http://web.bryant.edu/~bblais">http://web.bryant.edu/~bblais</a>




# Main script:

from dohebb import *
import pylab as p
from Numeric import *
from RandomArray import *
import time

x=random((100,1000))    # 1000 input vectors

numpats=x.shape[0]
w=random((numpats,1));

th=random((1,1))

params={}
params['eta']=0.001;
params['tau']=100.0;
old_mx=0;
for e in range(100):

     rnd=randint(0,numpats,250000)
     t1=time.time()
     if 0:  # straight python
         for i in range(len(rnd)):
             pat=rnd[i]
             xx=reshape(x[:,pat],(1,-1))
             y=matrixmultiply(xx,w)
             w=w+params['eta']*(y*transpose(xx)-y**2*w);
             th=th+(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()


#=============================================

#         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=params['eta']  # learning rate
     tau=params['tau']  # used for variance estimate

     cdef double y
     num_iterations=rnd.dimensions[0]
     num_inputs=w.dimensions[0]

     # get the pointers
     wp=&lt;double *&gt;w.data
     xp=&lt;double *&gt;X.data
     rndp=&lt;int *&gt;rnd.data
     thp=&lt;double *&gt;th.data

     for it from 0 &lt;= it &lt; num_iterations:

         offset=rndp[it]*num_inputs

         # calculate the output
         y=0.0
         for i from 0 &lt;= i &lt; num_inputs:
             y=y+wp[i]*xp[i+offset]

         # change in the weights
         for i from 0 &lt;= i &lt; num_inputs:
             wp[i]=wp[i]+eta*(y*xp[i+offset] - y*y*wp[i])

         # estimate the variance
         thp[0]=thp[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!
<a class="moz-txt-link-freetext" href="http://sel.as-us.falkag.net/sel?cmd=lnk&kid=103432&bid=230486&dat=121642">http://sel.as-us.falkag.net/sel?cmd=lnk&amp;kid=103432&amp;bid=230486&amp;dat=121642</a>
_______________________________________________
Numpy-discussion mailing list
<a class="moz-txt-link-abbreviated" href="mailto:Numpy-discussion@lists.sourceforge.net">Numpy-discussion@lists.sourceforge.net</a>
<a class="moz-txt-link-freetext" href="https://lists.sourceforge.net/lists/listinfo/numpy-discussion">https://lists.sourceforge.net/lists/listinfo/numpy-discussion</a>

      </pre>
    </blockquote>
    <pre wrap="">
-------------------------------------------------------
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!
<a class="moz-txt-link-freetext" href="http://sel.as-us.falkag.net/sel?cmd=k&kid3432&bid#0486&dat1642">http://sel.as-us.falkag.net/sel?cmd=k&amp;kid3432&amp;bid#0486&amp;dat1642</a>
_______________________________________________
Numpy-discussion mailing list
<a class="moz-txt-link-abbreviated" href="mailto:Numpy-discussion@lists.sourceforge.net">Numpy-discussion@lists.sourceforge.net</a>
<a class="moz-txt-link-freetext" href="https://lists.sourceforge.net/lists/listinfo/numpy-discussion">https://lists.sourceforge.net/lists/listinfo/numpy-discussion</a>




    </pre>
  </blockquote>
  <pre wrap=""><!---->
  </pre>
</blockquote>
</body>
</html>