From: davis at space.mit.edu (John E. Davis)
Date: 19 May 1999 10:47:27 GMT
Subject: my loop is too slow
References: <Pine.SGI.3.91.990510155017.9052A-100000@dft> <joe-110519991325429671@chinacat.salk.edu> <Oc1_2.1163$G3.34761@news14.ispnews.com>
Message-ID: <slrn7k4gsf.hgp.davis@mygir.davis.net>
X-UID: 1714                                               

On Tue, 11 May 1999 14:24:39 -0600, Tim Hochberg <hochberg at psn.net>
wrote:
>FWIW: Taking into account Michael Haggerty's observation on that this is
>really a matrixmultiply of (c.f)(c.transpose(f)), I reimplemented the
>original loop as:
>
>from Numeric import *
>a = a + dot(c, f) * dot(c, transpose(f))

Are you sure about this?  The original code was

  a[k,i]=a[k,i]+c[k,h]*c[k,j]*f[j,i]*f[i,h]

summed over all indices.  I would write this as:

  a_ki = a_ki + c_kh c_kj f_ji f_ih
       = a_ki + c_kh (c#f)_ki f_ih
       = a_ki + (c#f)_ki (c#f')_ki

where `#' denotes matrix multiplication and ' represents the
transpose.  So, I would code this as:

  a = a + (c#f)*(c#f')

where `*' represents an element by element multiply and `+' represents
an element by element addition.

--John





