first of all, I must make it clear that I don't want to stir up any battle between languages! Never! 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. the following is written by lac. I hope he will not charge me for my post ;) ####### speed ######### # python
def xperm(L): if len(L) <= 1: yield L else: for i in range(len(L)): for p in xperm(L[:i]+L[i+1:]): yield [L[i]] + p
def test(i): for p in xperm(map(lambda x:x+1, range(i))): pass;
import timeit timeit.Timer('test(5)','from __main__ import test').timeit()
1263.9717730891925 # MatLab - timeit.m function [ret] = timeit(times, repeat) ret = zeros(1,repeat); for i = 1:repeat tic(); for j = 1:times perms([1:5]); end ret(i) = toc(); end
timeit(1000000,3)
ans = 468.0005 449.7487 466.5325 ####### correctness ####### # MatLab
A = hilb(256);
[L,U,P] = lu(A);
norm(P*A - L*U)/norm(A)
ans = 3.2428e-017
[Q,R] = qr(A)
norm(Q)
ans = 1.0000 # SciPy
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
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
On Monday 08 June 2009 11:19:12 oyster wrote:
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.
the following is written by lac. I hope he will not charge me for my post ;)
If he does, refuse to pay, for this code is over 8 times slower than it should be! See below.
####### speed #########
# python
def xperm(L):
if len(L) <= 1: yield L else: for i in range(len(L)): for p in xperm(L[:i]+L[i+1:]): yield [L[i]] + p
def test(i):
for p in xperm(map(lambda x:x+1, range(i))): pass;
import timeit timeit.Timer('test(5)','from __main__ import test').timeit()
1263.9717730891925
from numpy import arange from numpy.random import permutation And then: In [17]: %timeit test(5) # Your version 1000 loops, best of 3: 444 µs per loop In [18]: %timeit test(5) 1000 loops, best of 3: 444 µs per loop First get rid of the map (low hanging fruit): In [19]: def test2(i): for p in xperm(range(1, i+1)): pass ....: In [21]: %timeit test2(5) 1000 loops, best of 3: 429 µs per loop In [22]: %timeit test2(5) 1000 loops, best of 3: 430 µs per loop In [23]: %timeit test2(5) 1000 loops, best of 3: 431 µs per loop Now use numpy and sets: In [38]: def xperm2(L): tried = set() while len(tried) < len(L): new_perm = tuple(permutation(L)) if new_perm in tried: continue tried.add(new_perm) yield new_perm In [45]: %timeit test2(5) 10000 loops, best of 3: 86.9 µs per loop In [46]: %timeit test2(5) 10000 loops, best of 3: 86.2 µs per loop And finally numpy.arange FTW: In [47]: def test2(i): for p in xperm2(arange(1, i+1)): pass In [51]: %timeit test2(5) 10000 loops, best of 3: 51.9 µs per loop In [52]: %timeit test2(5) 10000 loops, best of 3: 51.4 µs per loop Anyone have a faster, vectorized way?
participants (3)
-
oyster
-
Pauli Virtanen
-
Yosef Meller