[SciPy-dev] MCMC, Kalman Filtering, AI for SciPy?

[Fernando Perez]

eric jones schrieb:
> Of course, we should get some testing feedback on the main platforms 
> with the current CVS before doing this.  Besides that, does anyone have 
> any specific (small) things to address before we do so?  If so, what is 
> the timeframe?

I have 2 things:

1. I'd like to have the implementation of Frobenius norms fixed.  The current 
code, in scipy/linalg/basic.py does this:


        elif ord in ['fro','f']:
             X = scipy_base.mat(x)
             return sqrt(sum(diag(X.H * X)))

This is a straightforward implementation of the formula:

frob_norm(A) = sqrt(Trace(A*hermitian_conjugate(A))),     [1]

where A is a _matrix_, so that * is a true matrix product (not element-wise).

Note, however, that this is equivalent to:

frob_norm(A) = sqrt(sum_i(sum_j(|A_ij|**2)))             [2]

But [1] is an O(N**3) calculation, while [2] is O[N**2] for an NxN matrix, 
since [1] requires a matrix multiplication.

The following will produce the right answer for complex matrices:

sqrt(sum((a*conjugate(a)).flat)).real

Here's a quick test of the difference, for a 1024x1024 complex matrix:

In [27]: timing (1,scipy.linalg.norm,a)
Out[27]: 7.7238250000000015

In [28]: def fnorm(a):
    ....:     return sqrt(sum((a*conjugate(a)).flat)).real
    ....:

In [29]: timing (1,fnorm,a)
Out[29]: 0.12298099999999934

In [30]: scipy.linalg.norm(a)-fnorm(a)
Out[30]: -2.2737367544323206e-13


Untested code valid for real/complex should probably be:

if a.typecode()=='D':
   return sqrt(sum((a*conjugate(a)).flat)).real
else:
   return sqrt(sum((a**2).flat))


2. The matrixmultiply != bug I recently reported is _very_ serious for people 
working with large matrices.  I also reported it in the numpy list, but nobody 
replied.  Since scipy can trivially work around it with the one-liner I 
showed, I think this really should be done.

Code for #1 is here (but test it further), and I posted the one-liner for #2 
yesterday, so these two (IMHO important) fixes should take no time for a 
committer.

Best,

f