Items tagged with sparse


I thought I would share some code for computing sparse matrix products in Maple.  For floating point matrices this is done quickly, but for algebraic datatypes there is a performance problem with the builtin routines, as noted in

Download spm.txt

The code is fairly straightforward in that it uses op(1,A) to extract the dimensions and op(2,A) to extract the non-zero elements of a Matrix or Vector, and then loops over those elements.  I included a sparse map function for cases where you want to map a function (like expand) over non-zero elements only.

# sparse matrix vector product
spmv := proc(A::Matrix,V::Vector)
local m,n,Ae,Ve,Vi,R,e;
n, m := op(1,A);
if op(1,V) <> m then error "incompatible dimensions"; end if;
Ae := op(2,A);
Ve := op(2,V);
Vi := map2(op,1,Ve);
R := Vector(n, storage=sparse);
for e in Ae do
n, m := op(1,e);
if member(m, Vi) then R[n] := R[n] + A[n,m]*V[m]; end if;
end do;
return R;
end proc:
# sparse matrix product
spmm := proc(A::Matrix, B::Matrix)
local m,n,Ae,Be,Bi,R,l,e,i;
n, m := op(1,A);
i, l := op(1,B);
if i <> m then error "incompatible dimensions"; end if;
Ae := op(2,A);
Be := op(2,B);
R := Matrix(n,l,storage=sparse);
for i from 1 to l do
Bi, Be := selectremove(type, Be, (anything,i)=anything);
Bi := map2(op,[1,1],Bi);
for e in Ae do
n, m := op(1,e);
if member(m, Bi) then R[n,i] := R[n,i] + A[n,m]*B[m,i]; end if;
end do;
end do;
return R;
end proc:
# sparse map
smap := proc(f, A::{Matrix,Vector})
local B, Ae, e;
if A::Vector then
B := Vector(op(1,A),storage=sparse):
B := Matrix(op(1,A),storage=sparse):
end if;
Ae := op(2,A);
for e in Ae do
B[op(1,e)] := f(op(2,e),args[3..nargs]);
end do;
return B;
end proc:

As for how it performs, here is a demo inspired by the original post.

n := 674;
k := 6;
A := Matrix(n,n,storage=sparse):
for i to n do
  for j to k do
    A[i,irem(rand(),n)+1] := randpoly(x):
  end do:
end do:
V := Vector(n):
for i to k do
  V[irem(rand(),n)+1] := randpoly(x):
end do:
C := CodeTools:-Usage( spmv(A,V) ):  # 7ms, 25x faster
CodeTools:-Usage( A.V ):  # 174 ms
B := Matrix(n,n,storage=sparse):
for i to n do
  for j to k do
    B[i,irem(rand(),n)+1] := randpoly(x):
  end do:
end do:
C := CodeTools:-Usage( spmm(A,B) ):  # 2.74 sec, 50x faster
CodeTools:-Usage( A.B ):  # 2.44 min
# expand and collect like terms
C := CodeTools:-Usage( smap(expand, C) ):


I have use a sparsematrixplot to identify the patterns in a matrix.

In order to facilitate the readability of my sparsematrixplot, i would like to change the labels of the rows/ columns axis.

Instead of the numeric graduation, i would like to add :

- for my rows, the labels : [eq1,eq2, ..., ]
- for my columns, the labels : [q1,q2,q3,...,]

Do you have ideas how I can do that ?

Thank you for your help


Consider the following code, which just generates two "identical" matrices, differing only in their requested storage type, and then does some simple manipulations.

# Define matrix using sparse storage
   testM:= Matrix( 40,40,
# Define identical(?) matrix with rectangular storage
   nm:= Matrix( 40,40,
# Define procedure to return some matrix properties
   matData:= proc( myMat::Matrix)
                            return op(3, myMat)[2], # check storage type
                                      myMat[5, 1..-1], # get 5-th row
                                      add(myMat[5, 1..-1]); # add elements in 5-th row
                    end proc:
# Get properies of the two matrices - should be identical
# but check result of adding elements in the 5-th row

The matData procedure ought to produce the same results for the two matrices, with the exception of the storrage type. But the 'add()' command does not. The 'myMat[5, 1..-1]' command produces the same vector, the 5-th row - but stick an add() wrapper around it and all hell breaks loose.

Is this a bug or am I missing something?

Suggestions such as avoiding sparse data storage are not really acceptable: the above is a much simplified version of my original problem where I was using graph theory to play with a "cost function" and (with G a graph) the command,


returned a sparse-storage matrix - and I didn't notice. There appears to be no option on the WeightMatrix() command to control the storage tyoe of the returned matrix. Result was that all subsequent code based on slicing/dicing/and particularly 'add()ing' sub-blocks of this weight matrix fell apart

Don't get me wrong: I can sort of accept that the weight matrix of minimal spanning tree would (hopefully) be mainly zeros so sparse-storage might be a good default option but I don't see why the results of a command such as

add(myMat[5, 1..-1])

should vary depending on the internal storage used for the matrix, particularly when I have no control over the storage type being adopted



       I have a huge sparse matrix (1.5e6 x 1.5e6) of both symbolic or non-symbolic form. When I export matrix using ExportMatrix command, it is not writing the matrix in order of rows and columns. To further explain the problem,

          Let A be that huge sparse matrix. I am using


We expect the output to be,

 1 1 0.1

 1 10 0.65




But it is writing the output as


100 1 0.25

100 25 0.65

Hence it is not writing in order. What is the way we can force it to write in order? I do understand using for loop will consume more time.


I have a small Finite Element program in Maple. I want to accelerate the Matrix assembly. For Matlab this works like described here ( under "Triplet Sparse Storage". There, two indexing Vectors are used to place the values into a sparse Matrix with the Matlab function sparse. If index pairs are duplicate, they are added.


Now I want to ask if this is also possible in Maple?...

I'm trying to compute the rank of large, non-quadratic matrices modulo p=2,3,5 and 7.

For matrix-sizes of (approx.) 20000 x 20000 everything works fine with Rank(p,M) of LinearAlgebra[Modular].

But for (approx.) 40000 x 40000 even the parsing of the matrix is not finished after 5 days (for 20k x 20k this takes 1 or 2 hours) . The system I use has enough ram, so no swaping is involved.


First question: I expected the parsing to be O(n^2), why is this not the case?

I am dissapointed. I recently upgraded to Maple 15 and I bought a license for the NAG libraries thinking that I could use NAG mark 8 chapter F12 to solve my Eigenvalue problem with large sparse matrix, in reality, there seems to be no implementation of these NAG routines in Maple. is there a way to get around this problem? I know I can also connect to MATLAB to take advantage of the integrated sparse handling but I want to avoid buying another piece of software(you know we scientist aren't millionaries).

I'm trying to compute the determinant of a 19x19 matrix, each with entries consisting of multivariate polynomials (polys in 6 vars, of degree at most 8).  Specifically I want to factor the determinant.  Running this on Maple 15 quickly gives an "object too large" error in the determinant function, after using only 700 MB of RAM in about 2 minutes.  My machine has 3GB of RAM, so I am confused as to why this object is "too large". 




How to load Maple-NAG connector to execut this example ?

Is it an additional program to install ?

n := 4:
nnz := 11:
la := 22:
lfill := 1:
dtol := 0:
pstrat := "Nag_SparseNsym_CompletePiv":
milu := "Nag_SparseNsym_UnModFact":
a := Vector[row](la,[1.0,1.0,-1.0,2.0,2.0,3.0,-2.0,1.0,-2.0,1.0,1.0],datatype=float[8]):
irow := Vector[row](22,[1,1,2,2,2,3,3,4,4,4,4],datatype=integer[kernelopts('wordsize')/8]):
icol := Vector[row...

I am using finite element analysis and have 47008 nodes, so require a 47008 by 47008 matrix.  I wish to create the matrix using a procedure and using if statements to fill in the matrix with the equations (mostly repeating and every element in each row is 0 excpet for 4 or 5).  Anyway, I have tested my solution method on a much smalller matrix and know it works.  My question is if there is a way to create such a large matrix without getting a memory error, and...

Our previous article described the design of fast algorithms for multiplying and dividing sparse polynomials. We have integrated these algorithms into the expand and divide commands of Maple 14. In this post I want to talk a bit about what you might see when you try Maple 14. Keep in mind that the product isn't released yet and I don't work for Maplesoft, so general disclaimers apply. Nevertheless, one of the first things you may notice is this.

task manager with maple 14

In our previous article we described a packed representation for sparse polynomials is designed for scalability and high performance. The expand and divide commands in Maple 14 use this representation internally to multiply and divide polynomials with integer coefficients, converting to and from Maple's generic data structure described here. In this post I want to show you how these algorithms work and why they are fast. It's a critical stepping stone for our next topic, which is parallelization.

sdmp multiplication

Our first article introduced Maple's polynomial data structure and explained how Maple spends a lot of time working with monomials. To multiply polynomials having n and m terms, Maple must construct, simplify, hash, and sort all nm pairwise products to determine what monomials are equal. This work is performed even if the result has far fewer than nm terms, making it a rather inefficient way to multiply large multivariate polynomials. This article describes a new data structure for multivariate polynomials that is being added to Maple for a future release.

sdmp packed arrays

9xyz  -  4yz  -  6xyz  -  8x  -  5

This is the first in a series of short informal articles about our efforts to speed up polynomial arithmetic in Maple. We begin with an example of how polynomials are represented in Maple right now.

Maple sum structure

9xyz  -  4yz  -  6xyz  -  8x  -  5

When you enter a polynomial in Maple, it creates a generic data structure like the one above. In Maple's representation this polynomial is a sum of terms that is 11 words of memory long where each word is either 32 or 64 bits. For each term it stores a pointer to a monomial followed by a coefficient.

I was reminded of this by another thread.

It is faster to add in-place a large size storage=sparse float[8] Matrix into a new empty storage=rectangular float[8] Matrix than it is to convert it that way using the Matrix() or rtable() constructors.

Here's an example. First I'll do it with in-place Matrix addition. And then after that with a call to Matrix(). I measure the time to execute as well as the increase in bytes-allocated and bytes-used.

> with(LinearAlgebra):

> N := 500:
> A := RandomMatrix(N,'density'=0.1,
>                   'outputoptions'=['storage'='sparse',
>                                    'datatype'=float[8]]):

> st,ba,bu := time(),kernelopts(bytesalloc),kernelopts(bytesused):

> B := Matrix(N,'datatype'=float[8]):
> MatrixAdd(B,A,'inplace'=true):

> time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
                            0.022, 2489912, 357907
1 2 Page 1 of 2