[SciPy-dev] MCMC, Kalman Filtering, AI for SciPy?
Fernando Perez
Fernando.Perez at colorado.edu
Thu Sep 30 15:00:29 CDT 2004
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
More information about the Scipy-dev
mailing list