[SciPy-Dev] scipy.sparse versus pysparse
nicky van foreest
vanforeest at gmail.com
Wed Jul 23 16:08:16 EDT 2014
Hi Moritz,
Sure. Please see below. I included an extra time stamp to analyse the
results in slightly more detail. It turns out that the matrix-vector
multiplications are roughly the same in scipy.stats and pysparse, but that
building the matrices in pysparse is way faster.
I compute the stationary distribution vector of a Markov chain. The details
of the algo are not really important. I guess that the logic of the code is
easy to understand, if not, please let me know. There are three nearly
identical methods to compute the distribution vector, the first uses
scipy.stats and the * operator to compute a vector times a matrix, the
second uses scipy.sparse dot(), and the third uses pysparse. You might
want to skip the code and jump right away to results which I include below
the code
====
Code
from numpy import ones, zeros, empty
import scipy.sparse as sp
import pysparse
from pylab import matshow, savefig
from scipy.linalg import norm
import time
labda, mu1, mu2 = 1., 1.1, 1.01
N1, N2 = 400, 400
size = N1*N2
eps = 1e-3
maxIterations = 1e5
print "size = ", size
def state(i,j):
return j*N1 + i
def fillOffDiagonal(Q):
# labda
for i in range(0,N1-1):
for j in range(0,N2):
Q[(state(i,j),state(i+1,j))]= labda
# mu2
for i in range(0,N1):
for j in range(1,N2):
Q[(state(i,j),state(i,j-1))]= mu2
# mu1
for i in range(1,N1):
for j in range(0,N2-1):
Q[(state(i,j),state(i-1,j+1))]= mu1
#print "ready filling"
def computePiMethod1():
"""
based on scipy.sparse, naive matrix-vector multiplication
"""
e0 = time.time()
Q = sp.dok_matrix((size,size))
fillOffDiagonal(Q)
# Set the diagonal of Q such that the row sums are zero
Q.setdiag( -Q*ones(size) )
# Compute a suitable stochastic matrix by means of uniformization
l = min(Q.values())*1.001 # avoid periodicity, see trivedi's book
P = sp.eye(size, size) - Q/l
P = P.tocsr()
pi = zeros(size); pi1 = zeros(size)
pi[0] = 1;
n = norm(pi - pi1,1); i = 0;
e1 = time.time()
while n > eps and i < maxIterations:
pi1 = pi*P
pi = pi1*P # avoid copying pi1 to pi
n = norm(pi - pi1,1); i += 1
print "Method 1: ", e1-e0, time.time() - e1, i
return pi
def computePiMethod2():
"""
based on scipy.sparse, dot multiplication
"""
e0 = time.time()
Q = sp.dok_matrix((size,size))
fillOffDiagonal(Q)
# Set the diagonal of Q such that the row sums are zero
Q.setdiag( -Q*ones(size) )
# Compute a suitable stochastic matrix by means of uniformization
l = min(Q.values())*1.001 # avoid periodicity, see trivedi's book
P = sp.eye(size, size) - Q/l
P = P.transpose()
P = P.tocsr()
pi = zeros(size); pi1 = zeros(size)
pi[0] = 1;
n = norm(pi - pi1,1); i = 0;
e1 = time.time()
while n > eps and i < maxIterations:
pi1 = P.dot(pi)
pi = P.dot(pi1)
n = norm(pi - pi1,1); i += 1
print "Method 2: ", e1-e0, time.time() - e1, i
return pi
def computePiMethod3():
"""
based on pysparse
"""
e0 = time.time()
Q = pysparse.spmatrix.ll_mat(size,size)
fillOffDiagonal(Q)
# fill diagonal
x = empty(size)
Q.matvec(ones(size),x)
Q.put(-x)
# uniformize
l = min(Q.values())*1.001
P = pysparse.spmatrix.ll_mat(size,size)
P.put(ones(size))
P.shift(-1./l, Q)
# Compute pi
P = P.to_csr()
pi = zeros(size); pi1 = zeros(size)
pi[0] = 1;
n = norm(pi - pi1,1); i = 0;
e1 = time.time()
while n > eps and i < maxIterations:
P.matvec_transp(pi,pi1)
P.matvec_transp(pi1,pi)
n = norm(pi - pi1,1); i += 1
print "Method 3: ", e1-e0, time.time() - e1, i
return pi
def plotPi(pi):
pi = pi.reshape(N2,N1)
matshow(pi)
savefig("pi.png")
if __name__ == "__main__":
pi1 = computePiMethod1()
pi2 = computePiMethod2()
pi3 = computePiMethod3()
d1 = norm(pi1-pi2,1)
d2 = norm(pi1-pi3,1)
print d1, d2from numpy import ones, zeros, empty
import scipy.sparse as sp
import pysparse
from pylab import matshow, savefig
from scipy.linalg import norm
import time
labda, mu1, mu2 = 1., 1.1, 1.01
N1, N2 = 400, 400
size = N1*N2
eps = 1e-3
maxIterations = 1e5
print "size = ", size
def state(i,j):
return j*N1 + i
def fillOffDiagonal(Q):
# labda
for i in range(0,N1-1):
for j in range(0,N2):
Q[(state(i,j),state(i+1,j))]= labda
# mu2
for i in range(0,N1):
for j in range(1,N2):
Q[(state(i,j),state(i,j-1))]= mu2
# mu1
for i in range(1,N1):
for j in range(0,N2-1):
Q[(state(i,j),state(i-1,j+1))]= mu1
#print "ready filling"
def computePiMethod1():
"""
based on scipy.sparse, naive matrix-vector multiplication
"""
e0 = time.time()
Q = sp.dok_matrix((size,size))
fillOffDiagonal(Q)
# Set the diagonal of Q such that the row sums are zero
Q.setdiag( -Q*ones(size) )
# Compute a suitable stochastic matrix by means of uniformization
l = min(Q.values())*1.001 # avoid periodicity, see trivedi's book
P = sp.eye(size, size) - Q/l
P = P.tocsr()
pi = zeros(size); pi1 = zeros(size)
pi[0] = 1;
n = norm(pi - pi1,1); i = 0;
e1 = time.time()
while n > eps and i < maxIterations:
pi1 = pi*P
pi = pi1*P # avoid copying pi1 to pi
n = norm(pi - pi1,1); i += 1
print "Method 1: ", e1-e0, time.time() - e1, i
return pi
def computePiMethod2():
"""
based on scipy.sparse, dot multiplication
"""
e0 = time.time()
Q = sp.dok_matrix((size,size))
fillOffDiagonal(Q)
# Set the diagonal of Q such that the row sums are zero
Q.setdiag( -Q*ones(size) )
# Compute a suitable stochastic matrix by means of uniformization
l = min(Q.values())*1.001 # avoid periodicity, see trivedi's book
P = sp.eye(size, size) - Q/l
P = P.transpose()
P = P.tocsr()
pi = zeros(size); pi1 = zeros(size)
pi[0] = 1;
n = norm(pi - pi1,1); i = 0;
e1 = time.time()
while n > eps and i < maxIterations:
pi1 = P.dot(pi)
pi = P.dot(pi1)
n = norm(pi - pi1,1); i += 1
print "Method 2: ", e1-e0, time.time() - e1, i
return pi
def computePiMethod3():
"""
based on pysparse
"""
e0 = time.time()
Q = pysparse.spmatrix.ll_mat(size,size)
fillOffDiagonal(Q)
# fill diagonal
x = empty(size)
Q.matvec(ones(size),x)
Q.put(-x)
# uniformize
l = min(Q.values())*1.001
P = pysparse.spmatrix.ll_mat(size,size)
P.put(ones(size))
P.shift(-1./l, Q)
# Compute pi
P = P.to_csr()
pi = zeros(size); pi1 = zeros(size)
pi[0] = 1;
n = norm(pi - pi1,1); i = 0;
e1 = time.time()
while n > eps and i < maxIterations:
P.matvec_transp(pi,pi1)
P.matvec_transp(pi1,pi)
n = norm(pi - pi1,1); i += 1
print "Method 3: ", e1-e0, time.time() - e1, i
return pi
def plotPi(pi):
pi = pi.reshape(N2,N1)
matshow(pi)
savefig("pi.png")
if __name__ == "__main__":
pi1 = computePiMethod1()
pi2 = computePiMethod2()
pi3 = computePiMethod3()
d1 = norm(pi1-pi2,1)
d2 = norm(pi1-pi3,1)
print d1, d2
==============================
Results
nicky at chuck:~/myprogs/python/queueing/tandemQueueMDP$ python tandemqueue.py
size = 40000
Method 1: 4.31593680382 0.387089014053 285
Method 2: 4.27599096298 0.273495912552 285
Method 3: 0.0856800079346 0.267058134079 285
0.0 5.05082123061e-15
The first number after "method1:" represents the time it takes to fill the
matrix, the second number is the time involved to carry out the
multiplications, and the third is (twice) the number of mutliplications
involved. The second number is (nearly) the same for all three methods,
hence the multiplication time is about the same. There is, however, a huge
difference in the first number, ie, the time required to build the matrix.
It takes about 4 sec for scipy.stats, and 8e-2 sec for pysparse.
The last row prints the difference between the results (the stationary
distributions vectors) as obtained by all three methods. Luckily the
results are the same, up to rounding.
Here is a try with a bigger matrix:
nicky at chuck:~/myprogs/python/queueing/tandemQueueMDP$ python tandemqueue.py
size = 160000
Method 1: 17.4650111198 1.80849194527 285
Method 2: 17.5270321369 1.54912996292 285
Method 3: 0.382800102234 1.63900899887 285
0.0 5.87925665394e-15
Again the same result.
Do you perhaps have any explanation for this?
Thanks
Nicky
On 23 July 2014 16:51, Moritz Beber <moritz.beber at gmail.com> wrote:
> Hey,
>
>
> On Wed, Jul 23, 2014 at 4:40 PM, nicky van foreest <vanforeest at gmail.com>
> wrote:
>
>>
>>
>> I am doing some testing between scipy.sparse and pysparse on my ubuntu
>> machine. Some testing reveals that pysparse is about 9 times faster in
>> matrix-vector multiplication that scipy.sparse. Might there be anything
>> specific I forgot to do during scipy's installation (I just ran apt-get
>> install python-scipy)? Is there another simple explanation for this
>> difference? I prefer to use scipy.sparse for its cleaner api, but a factor
>> 9 in speed is considerable.
>>
>>
> Could you post your benchmarking code somewhere (or show here), please?
>
> Cheers,
> Moritz
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20140723/681b5262/attachment.html>
More information about the SciPy-Dev
mailing list