[SciPy-user] compare scipy and matlab

[Pauli Virtanen]

Mon, 08 Jun 2009 16:19:12 +0800, oyster kirjoitti:
[clip]
> what I want to know is the truth and availabilty to use scipy in serious
> works. Because I have met scipy bug but it has been fixed soon; and now
> I met 2(maybe 1) bugs with sympy, so I regard opensoftware(or any
> software?) with some suspicion.
[clip]
> ####### speed #########
[clip]

The algorithms in Matlab perms.m and the one you wrote are not 
comparable. The Matlab one is vectorized, whereas the one you wrote is a 
naive implementation, and this should account for much of the performance
difference.

> ####### correctness #######
[clip]

>>> import numpy, scipy
>>> Hilb = lambda N: scipy.mat( [ [ 1.0/(i + j + 1) for j in range(N) ] for i in range(N) ] )
>>> A = Hilb(256)
>>> P, L, U = scipy.linalg.lu(A)
>>> scipy.linalg.norm( scipy.mat(P)*A - scipy.mat(L)*scipy.mat(U) ) / scipy.linalg.norm(A)
0.56581114936540999
>>> Q, R = scipy.linalg.qr(A)
>>> scipy.linalg.norm(scipy.mat(Q))
16.000000000000011

The default matrix norm is the Frobenius norm, not the 2-norm.
If you write the above using the correct norm, you get the same
as from Matlab. Also, note that scipy.linalg.lu uses a different
definition of the permutation matrix than Matlab's LU:

>>> scipy.linalg.norm(A - scipy.mat(P)*scipy.mat(L)*scipy.mat(U), 2) / scipy.linalg.norm(A, 2)
1.8016647238506838e-17
>>> scipy.linalg.norm(scipy.mat(Q), 2)
1.0000000000000007

Based on the above, I do not see any problems with either reliability or speed.

-- 
Pauli Virtanen