[Numpy-discussion] Enhancing dot()

[Tim Hochberg]

Bill Baxter wrote:
> On 7/7/06, *Tim Hochberg* <tim.hochberg at cox.net 
> <mailto:tim.hochberg at cox.net>> wrote:
>
>     I'd caution here though that the H is another thing that's going to
>     encourage people to write code that's less accurate and slower than it
>     needs to be. Consider the simple equation:
>
>     [Side note: Here's something I periodically think about -- it
>     would be
>     nice if dot could take multiple arguments and chain them into a series
>     of matrix multiplies so that dot(dot(a, b), a.H) could be written
>     dot(a,
>     b, a.H). I'll use that here for clarity].
>
>     x = dot(a, b, a.H)
>
>     versus:
>
>     x  = dot(a.real, b, a.real.transpose()) + dot(a.imag, b,
>     a.imag.transpose())
>
>     The latter should be about twice as fast and won't have any pesky
>     imaginary residuals showing up from the finite precision of the
>     various
>     operations.
>
>
> The funny thing is that having a dot(a,b,c,...) would lead to the 
> exact same kind of hidden performance problems you're arguing against.
Not exactly arguing -- this isn't why I don't like H and friends -- just 
noting that this is one of the traps that people are likely to fall into 
when transferring equations to code.

> The cost of a series of matrix multiplications can vary dramatically 
> with the order in which you perform them.  E.g. A*B*C*v where A,B,C 
> are 2-dim and v is a column vector.  Then you should do (A*(B*(C*v))), 
> and definitely not ((A*B)*C)*v)
That's a good point.
> But that said, dot could maybe be made smart enough to figure out the 
> best order.   That would be a nice feature. 
> dot(A, dot(B, dot(C, v))) is pretty darn ugly.  I seem to remember 
> finding the best order being an example of a dynamic programming 
> problem, and probably O(N^2) in the number of terms.  But presumably N 
> is never going to be very big, and the N==2 (no choice) and N=3 (only 
> one choice) cases could be optimized.
You could probably just brute force it with out too much trouble. As you 
say O(n^2) is not that daunting when n is usually 2 or 3 and rarely over 
5. Perhaps cache the last few configurations tested since the same 
matrix size combinations will likely occur over and over.

Another interesting thought is one just wrapped an array with a wrapper 
that says "use the conjugate of what I'm wrapping", then dot(a, b H(a)) 
could do the write thing and only compute the real part of the 
equation,  assuming H(a) created the wrapper. That's starting to get 
complex though.


-tim