my loop is too slow

[John E. Davis]

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