Sparse Matrix-Vector product using evaluation

Posted on 2006-04-11 04:04 By roman_pearce ( Maple Rating 4

508)

Categories:

This tip comes care of Dr. Michael Monagan at Simon Fraser University. Represent your sparse matrix as a list of rows, and represent each row as a linear equation in an indexed name. For example:

A := [[1,0,3],[2,0,0],[0,4,5]];
S := [ 1*x[1] + 3*x[3], 2*x[1], 4*x[2]+5*x[3] ];

To compute the product of the matrix A with a Vector X, assign x[i] := V[i] and evaluate. This can be done inside of a procedure because x is a table.

V := [7,8,9]:
for i to 3 do x[i] := V[i] end do:
eval(S); # result
for i to 3 do x[i] := evaln(x[i]) end do: # unassign

Don't forget to unassign. Note that this method of assignment is slow when V is a big list. In short, we need to use "for i in V", which does one nice linear pass through V. This is what I do:

n := 0;
count := proc() global n; n := n+1 end proc:
for i in V do x[count()] := i end do:
eval(S);
for i to 3 do x[i] := evaln(x[i]) od:

In a real procedure, n would be lexically scoped and not a global variable. Instead of an explicit loop you can also use seq and assign.

n := 0;
count := proc() global n; n := n+1 end proc:
seq(assign(x[count()],i), i=V);
eval(S);
seq(unassign('x[i]'), i=1..3);

Of course, where this trick really shines is in doing back substitution. You can easily substitute tens of thousands of variables into any number of equations, all in linear time.

Printer-friendly version
Login or register to post comments

Comments

Sparse Matrix-Vector Product: Application and Timing Comparison

Comment on 2006-04-13 10:50 from Joe Riel ( Maple Rating 4

1260)

This is an interesting technique. To see a benefit, I presume, the multiplication by the matrix must be done repeatedly for different vectors, otherwise the start-up cost is prohibitive. An obvious application is Matrix-Matrix multiplication. I put together a small procedure to test this:

MatrixMultiply := proc(A,B)
local m,n,p,X,Ax,i,j,x,col;
uses LinearAlgebra:-LA_Main;
    m := op([1,1],A);
    (n,p) := op(1,B);
    X := rtable( 1..n
                 , i -> x[i]
                 , 'subtype' = 'Vector'['column']
               );
    Ax := convert(MatrixVectorMultiply(A, X, 'outputoptions' = [])
                  , list);
    for j to p do
        for i to m do
            x[i] := B[i,j];
        end do;
        col[j] := eval(Ax);
        x := 'x';
    end do;
    rtable( 1..m, 1..p
            , ListTools:-Transpose([seq](col[j],j=1..p))
            , 'subtype' = 'Matrix'
          );
end proc:

Timing

Here is a direct comparison with Matrix multiplication using the LA:-Main subpackage (to avoid the minor overhead of type checking, etc):

with(LinearAlgebra):
(m,n,p) := (100,100,1000);
B := RandomMatrix(n,p):
A := RandomMatrix(m,n,density=1):
time(MatrixMultiply(A,B));
time(LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(A,B,'inplace'=false,'outputoptions'=[]));
A := RandomMatrix(m,n,density=1/10):
time(MatrixMultiply(A,B));
time(LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(A,B,'inplace'=false,'outputoptions'=[]));
A := RandomMatrix(m,n,density=1/100):
time(MatrixMultiply(A,B));
time(LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(A,B,'inplace'=false,'outputoptions'=[]));

The times were

density     MatrixMultiply   LA_Main
------------------------------------
   1            5.340         4.269 
  1/10          0.910         3.949
  1/100         0.450         3.840

For dense Matrices there is a minor penalty to this technique. For sparse Matrices there can be a considerable time saving.

Follow up
I corrected a bug in the code. A for-loop over an rtable, in the previous version, for v in B[1..-1,j] is not guaranteed to loop sequentially through the elements. It gets all the elements, but not necessarily in order of increasing index. If the selected element were not equal to B[i,j], then an incorrect result would be returned. There was no clear change in the execution time.

Actually I was thinking...

Comment on 2006-04-28 02:32 from roman_pearce ( Maple Rating 4

508)

Actually I was thinking more along the lines of using this as a new matrix data structure...