Sparse matrix manipulations
- To: mathgroup at smc.vnet.net
- Subject: [mg32716] Sparse matrix manipulations
- From: Hugh Goyder <goyder at rmcs.cranfield.ac.uk>
- Date: Thu, 7 Feb 2002 05:09:57 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
I would like to speed up my calculations by using sparse matrix methods
instead of using full matrices with many zero elements. Mathematica
provides a sparse method for solving smat . x == vec
?Developer`SparseLinearSolve
"SparseLinearSolve[smat, vec] solves a sparse linear system; the matrix smat
is represented in the form {{i1, j1}->a1, {i2, j2}->a2, ... }, so that the
element at position ik, jk has value ak and all unspecified elements are
taken to be zero."
I would like to do addition, multiplication and inversion with sparse
matrices that are numeric, have 97% or more zero elements and have
dimensions of 1000 x 1000 or more.
I wrote the following modules which are not as fast as I hoped. Are there
standard modules for these calculations? Are the best modules always
written for specific applications?
I am certain that you experts can do a lot better.
1.1 Addition
Add to one element, if it is non zero, or make a new element if it is zero.
SparseAddToElement[smat_, {i_, j_} -> val_] := Module[{new},
If[Flatten[Position[smat, ({i, j} -> a_)]] == {},
new = Join[smat, {({i, j} -> val)}],
new = smat /. ({i, j} -> t_) -> ({i, j} -> t + val)]; new]
1.2 Addition
Add sparse matrices m1 and m2
SparseAdd[m1_, m2_] := Module[{m = m1},
Do[m = SparseAddToElement[m, m2[[i]]], {i, Length[m2]}]; m]
2. Multiplication
Two sparse matrices m1 and m2 are multiplied together. The issue is to
select and pair off the row and column elements, if they are non zero.
SparseMultiply[m1_, m2_] :=
Module[{numrows = Max[#[[1, 1]] & /@ m1], numcols = Max[#[[1, 2]] & /@ m2],
rr, cc, ii, val, pos, result = {}},
Do[rr = Select[m1, #[[1, 1]] == i &]; cc = Select[m2, #[[1, 2]] == j &];
ii = 0; val = 0;
While[ii = ii + 1; ii <= Length[rr],
If[pos = Position[cc, {rr[[ii, 1, 2]], a_} -> b_]; pos != {},
{{jj}} = pos; val = val + rr[[ii, 2]] cc[[jj, 2]]]];
If[val != 0, result = {result, {i, j} -> val}],
{i, numrows}, {j, numcols}];
Flatten[result]]
3. Inverse
Solve the sparse equation A X == B where A and B are sparse but it is
anticipated that the unknown matrix X will be a full matrix. I use the
Developer`SparseLinearSolve repeatedly presumably redoing much of the
decomposition of A for each column of X.
SparseInverse[A_, B_] :=
Module[{numrows = Max[#[[1, 1]] & /@ A], numcols = Max[#[[1, 2]] & /@ B],
col, i, j, result = {}},
col = Table[{i, j}, {i, numrows}];
Do[result = {result,
Developer`SparseLinearSolve[
A, ((col /. j -> jj) /. B) //. {a_, b_} -> 0]}, {jj, numcols}];
Partition[Flatten[result], numcols]]
Test using tridiagonal matrices following the logic of Daniel Lichtblau on
a previous post.
In[14]:=
tridiagonal[n_]:={Table[Random[],{n-1}],Table[Random[],{n}],
Table[Random[],{n-1}]};
toSparseMat[td_]:=
With[{n=Length[td[[2]]]},
Join[Table[{i,i+1}\[Rule]td[[1,i]],{i,n-1}],
Table[{i,i}\[Rule]td[[2,i]],{i,n}],
Table[{i+1,i}\[Rule]td[[3,i]],{i,n-1}]]]
In[52]:=
n=100; (* I would like n to be 1000 or more *)
SeedRandom[1111];
td=tridiagonal[n]; A=toSparseMat[td];
td=tridiagonal[n]; B=toSparseMat[td];
In[56]:=
Timing[SparseMultiply[A,B];]
Out[56]=
{250.52 Second,Null}
In[57]:=
Timing[X=SparseInverse[A,B];]
Out[57]=
{21.37 Second,Null}
In[58]:=
(* Make regular matrix out of A *)
AA=(Table[{i,j},{i,n},{j,n}]/.A)/.{i_,j_}\[Rule]0;
(* Make sparse matrix out of -AA . Transpose[X] *)
B1=Chop[AA.Transpose[X]];B1neg=
Select[Flatten[
Table[{i,j}\[Rule]-B1[[i,j]],{i,n},{j,n}]],#[[2]] \[NotEqual] 0&];
In[62]:=
(* B + AA.X, should be zero *)
Timing[res =SparseAdd[B,B1neg];]
Out[62]=
{2.75 Second,Null}
In[63]:=
Select[res//Chop,#[[2]] \[NotEqual] 0&]
Out[63]=
{}
All guidance on the above gratefully received.
Hugh Goyder