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.

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]

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 *)

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

```

• Prev by Date: Re: Newbie Question
• Next by Date: Does a Mathematica IRC channel exist?
• Previous by thread: Re: Moving notebooks from Windows to Mac
• Next by thread: Does a Mathematica IRC channel exist?