MathGroup Archive 2002

[Date Index] [Thread Index] [Author Index]

Search the Archive

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




  • 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?