ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Tue, 26 May 2020 18:00:39 +0200Parallelizing symbolic matrix*vector computations for sparse matrices and vectorshttps://ask.sagemath.org/question/51560/parallelizing-symbolic-matrixvector-computations-for-sparse-matrices-and-vectors/Hi everyone,
I am currently running computations which involve a lot of operations of the type matrix*vector, where both objects are sparse and the result can be expected to be sparse as well. The matrices' base field is the symbolic ring. If I drop certain normalizations I could also work over QQ, I think.
The above operations appear sequentially. Although there are some options for parallelizing these at least in part, the biggest acceleration would be achieved if the multiplication itself could be parallelized. Is there a way to do this in Sage?
Greetings and thank you very much.
Edit: Work in progress
Below I put some code that I produced so far which is very far from a solution but produces new slow-downs that I do not understand. The strategy is to save a sparse matrix by its nonzero rows which is facilitated by the following code snippet:
def default_to_row(M):
rows_of_M=M.rows()
relevant=[]
for r in range(len(rows_of_M)):
row=rows_of_M[r]
if dot_sparse(row,row)!=0:
relevant.append([r,row])
return relevant
The row*vector multiplications are facilitated by the function
# indexed_row of type [i,v] where v is a sparse vector and i is
# the row in the matrix it corresponds to.
def mult_sparse(indexed_row,v2):
res=0
for k in indexed_row[1].dict():
res=res+indexed_row[1][k]*v2[k]
return [indexed_row[0],res]
and regular vector multiplication is given by
def dot_sparse(v,w):
res=0
for k in v.dict():
res=res+v[k]*w[k]
return res
Now to the actual benchmarks:
import concurrent.futures
import time
import itertools
dim=5000
M=random_matrix(QQ,nrows=dim,ncols=dim,sparse=True,density=0.0001)
M_by_rows=default_to_row(M)
v=vector(QQ,dim,list(range(dim)),sparse=True)
print(len(M_by_rows))
Using the concurrent.futures library yields the following performance:
start = time.perf_counter()
with concurrent.futures.ProcessPoolExecutor() as executor:
start2 = time.perf_counter()
results = executor.map(mult_sparse,M_by_rows,itertools.repeat(v) )
finish2 = time.perf_counter()
coefficients={}
for result in results:
if result[1]!=0:
coefficients[result[0]]=result[1]
new_v=vector(QQ,v.degree(),coefficients,sparse=True)
finish = time.perf_counter()
print(f'Finished in {round(finish-start, 2)} seconds with Pooling.')
print(f'The executor alone took {round(finish2-start2, 2)} seconds.')
where the dimensions are adjusted so that the computation does not exceed my RAM capacities because every spawned subprocess takes copies (at least of v). A comparison with serial approaches, both the standard one and the serial version of my particular straegy, show that the parallel version is significantly slower (18 and 7 seconds in the above outputs) while the serial versions do not differ (0.02 seconds in both cases):
start = time.perf_counter()
new_v=M*v
finish = time.perf_counter()
print(f'Finished in {round(finish-start, 2)} second(s) for the regular computation')
start = time.perf_counter()
coefficients={}
for row in M_by_rows:
coefficients[row[0]]=mult_sparse(row,v)[1]
new_v=vector(QQ,v.degree(),coefficients,sparse=True)
finish = time.perf_counter()
print(f'Finished in {round(finish-start, 2)} second(s) in a serial computation')
print(f'Check if the vectors coincide:')
test=new_v-M*v
n=dot_sparse(test,test)
if n==0:
print(True)
else:
print(False)Robin_FTue, 26 May 2020 18:00:39 +0200https://ask.sagemath.org/question/51560/Kernel died when computing nullspace of a large matrixhttps://ask.sagemath.org/question/48650/kernel-died-when-computing-nullspace-of-a-large-matrix/I have a large sparse matrix with dimension roughly 100,000 by 70,000. Its only nonzero entries are 1 and -1. I'm trying to use the right_kernel command to calculate its kernel over the rational numbers. However, whenever I run the code, I get a message that says "the kernel appears to have died. It will restart automatically." Does this mean I'm running out of memory on my computer? Any ideas on how to address this? All this is done in a Jupyter notebook. Thank you!ClaudiaYunWed, 06 Nov 2019 18:15:57 +0100https://ask.sagemath.org/question/48650/Inverse of real sparse matrixhttps://ask.sagemath.org/question/47587/inverse-of-real-sparse-matrix/Dear sagemath community,
I'm a bit surprised in the result that I obtain when I compute the inverse of a real, sparse matrix. The problem occurs in Sagemath 8.7, (on windows).
In my case the matrix is
B=matrix(RR, [[1/6, -1/24, -1/30, 1/120,1/12, 0, 0, 0, 0],
[-1/24,1/60,1/60, 1/420, -1/24, 0, 0, 0, 0],
[-1/30,1/60, 2/105, 1/140, -1/20, 0, 0, 0, 0],
[1/120, 1/420, 1/140, 13/1260, -1/40, 0, 0, 0, 0],
[1/12, -1/24, -1/20, -1/40, 1/3, -1/24, -1/30, 1/120,1/12],
[0, 0, 0, 0, -1/24,1/60,1/60, 1/420, -1/24],
[0, 0, 0, 0, -1/30,1/60, 2/105, 1/140, -1/20],
[0, 0, 0, 0, 1/120, 1/420, 1/140, 13/1260, -1/40],
[0, 0, 0, 0,1/12, -1/24, -1/20, -1/40, 1/6]],sparse=true)
(B.inverse()*B).norm(1)
The result is 138.4999999999923, which is very far from the expected value: 1. If one replaces sparse=true by sparse=false, the result becomes reasonable.The problem also does not occur when I replace RR by QQ.
I suspect that this will not happen in all versions of Sagemath, since I stumbled upon this strange behaviour when I was executing a code that used to work fine years ago. Sadly I do not recall which version I used to have back then. I wanted to ask: what is currently the right way to compute the inverse of a sparse real matrix?
cheers,
IVIVSun, 25 Aug 2019 20:39:38 +0200https://ask.sagemath.org/question/47587/Concatenate two sparse matrices over a Polynomial Ringhttps://ask.sagemath.org/question/47205/concatenate-two-sparse-matrices-over-a-polynomial-ring/ Hi all,
I am wondering if there is a way to concatenate two sparse matrices over a polynomial ring. Here is a toy example that I am testing with:
S.<x,y> = PolynomialRing(QQ,2)
udict = {(0,0):1, (0,1):x, (1,1):y}
vdict = {(0,0):2*x,(1,0):y}
U = matrix(4,4,ud,sparse=True)
V = matrix(4,2,vd,sparse=True)
I would like to make a larger matrix that is $4 \times 6$, that has the block form $[U|V]$. I attempted to use the block_matrix command, but it said that the given number of rows and columns were incompatible with the number of submatrices.
In Python, for a sparse matrix (csr_matrix or coo_matrix), there is the command hstack from scipy.sparse. Is there a similar command where I can concatenate a list of these sparse matrices? This would be ideal for my general problem.
Thank you for your time!BarkThu, 18 Jul 2019 14:00:49 +0200https://ask.sagemath.org/question/47205/NEW VERSION: Sparse matrix stores something for each row?https://ask.sagemath.org/question/44709/new-version-sparse-matrix-stores-something-for-each-row/**NEW VERSION. This is what this question has boiled down to upon further investigation. Everything below this paragraph in boldface and the subsequent code is the old version and serves purely historical purposes =) I hope this was the right thing to do rather than creating a new question. The issue is that the amount of memory occupied by a sparse matrix seems to depend linearly on the number of rows in the matrix. In particular, my setup wont even let me create a zero $10^9\times 1$ sparse matrix but has no problem with a zero $1\times 10^9$ sparse matrix. This can all be illustrated by running the below code (and it's hardly the way things are supposed to be).**
print get_memory_usage()
B=matrix(QQ, 1, 1000000000, {})
print get_memory_usage()
B=matrix(QQ, 10000000, 1, {})
print get_memory_usage()
B=matrix(QQ, 1000000000, 1, {})
print get_memory_usage()
---------------------------------------------------------------------------------------------------------------------------------
I'm running the following code. Here I create a sparse $10^5\times 10^5$ identity matrix and then apply it repeatedly to a vector in $\mathbb R^{10^5}$ with 100 nonzero coordinates (which is stored as $10^5\times 1$ sparse matrix). Each time I add the result to a list until I run out of 2 GB of memory.
A=matrix(QQ, 100000, 100000, {(i,i):1. for i in range(100000)})
print get_memory_usage()
B=[matrix(QQ, 100000, 1, {(i,0):1. for i in range(100)})]
while (get_memory_usage()<2000): B.append(A*B[-1])
print len(B)
print get_memory_usage()
del B
del A
print get_memory_usage()
I'm receiving (on a freshly started kernel)
204.54296875
196
2003.51953125
1399.0234375
This raises two questions.
1. Why is there so much memory (1.4 GB) still in use after I successfully ran the code and deleted both objects I've created? That's a leak, right?
2. Why does deleting a list of 196 sparse matrices with 100 nonzero elements each free up 600 MB? Each such matrix should only take up a few KB, right?
I'm on Windows 8.1/SAGE 8.4.
**UPDATE**. Transposing the matrices, i.e. writing
...
B=[matrix(QQ, 1, 100000, {(0,i):1. for i in range(100)})]
for i in range(200): B.append(B[-1]*A)
...
seems to work well memory-wise, it returns
214.50390625
201
215.53125
186.62109375
However, it takes up much more time than the first version. This is probably due to the implementation of sparse matrix multiplication unfortunately containing a [loop over the columns of the right matrix](https://github.com/sagemath/sage/blob/6187d261eca3c980e575b53d1a31f580ba8cfdfd/src/sage/matrix/matrix_rational_sparse.pyx#L184). Is there any simple way around this high memory/low speed dilemma?
**UPDATE 2.** This might not be a memory leak and have more to do with the implementation of sparse matrices in general rather than their multiplication in particular. Apparently, a sparse matrix stores something for *each* of its rows as shown by
print get_memory_usage()
B=matrix(QQ, 1, 10000000, {})
print get_memory_usage()
B=matrix(QQ, 10000000, 1, {})
print get_memory_usage()
165.45703125
165.51953125
1082.68359375
This has got to be a known issue. I was not able to find a discussion, however. (This might be what is known as `csr_matrix` in `scipy` but why this would be chosen as the general standard here is beyond me.)imakhlinWed, 19 Dec 2018 01:20:27 +0100https://ask.sagemath.org/question/44709/Multiplying sparse matriceshttps://ask.sagemath.org/question/43979/multiplying-sparse-matrices/Can Sage multiply sparse matrices? If so, does it require some special syntax? I've been trying to run the following code:
M = Matrix(QQ, 1000001, 1000001, {(3,2): 27, (2,98): 71})
M2 = M*M
and it hasn't been working. If M is being converted to a dense matrix for multiplication then of course I wouldn't expect this to run. But it shouldn't be hard to do as sparse matrices, right?PWThu, 18 Oct 2018 20:49:55 +0200https://ask.sagemath.org/question/43979/span of a matrixspacehttps://ask.sagemath.org/question/42313/span-of-a-matrixspace/given a matrix subspace of sparse matrices over a ring Q
M=MatrixSpace(Q,1000,1000, sparse=True)
I need a way to define the subspace N of M given a list of matrices of M
I tried to turn matrices into lists and use the VectorSpace category, however my matrices are sparse matrices, and this takes very long
def lista(A): #converts a matrix into a list
return A.list()
def coor(P): # for the determined base B, it gives the coordinate vector of a matrix P[0]
A=base(P[1][1])
k=len(P[1][1])
V=VectorSpace(Q,2**(2*k), sparse=True) #I would like to work with MatrixSpace instead
B=[] #this will be my basis
for i in range(len(A)):
a=A[i][0]
B.append(V(lista(a)))
W=W=V.span_of_basis(B) #I dont know if this function exists in MatrixSpace
p=W.coordinate_vector(V(lista(P[0])))
return pgrburrulFri, 11 May 2018 07:47:52 +0200https://ask.sagemath.org/question/42313/sparse differential of a chain complexhttps://ask.sagemath.org/question/10415/sparse-differential-of-a-chain-complex/I'm trying to convert a simplicial complex into a chain complex and then extract the differential map. However the simplicial complex is rather large and sage runs out of memory when it tries to make the conversion. However, the differential should be a very sparse matrix and there should be plenty of memory to store a sparse representation of the differential. Is it possible to have sage compute the differential as a sparse matrix? Specifically the method which gives the out of memory error is `chain_complex()`.parsons.kyle.89Mon, 05 Aug 2013 18:16:29 +0200https://ask.sagemath.org/question/10415/