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