A(BC)=(AB)C
It's a property I use all the time, but typically never consciously think about. For instance, in my sparse matrix array classes, the sparse matrices are stored in reverse order because they are multiplied in reverse order. This is for efficiency reasons. If you have a set of matrices, {Ai}, and you want to multiply them through and then multiply the result with a vector, as follows:
A1A2A3 ... An-1An v (1)
then this:
A1(A2(A3( ... (An-1(An v)) ... ))) (2)
is more efficient than multiplying them through in the order in which they are written:
(( ... ((A1A2)A3) ... An-1)An v) (3)
because each intermediate result is a vector, rather than a matrix. Similarly, if v is a full matrix instead of a vector, (2) is still more efficient than (3) because multiplying a sparse matrix with a full one is more efficient than multiplying two sparse matrices.
Is associativity related to linearity? No, because we can define a new, linear operator, lets call it, #:
A#B = ABT = Σk aikbik
that loses the associative property. Now lets rewrite (1) in terms of #:
A1#(A2#(A3( ... (An-1#(An vT)T)T ... )T)T)T
= A1#vT#An#An-1# ... #A3#A2
= (vT#An#An-1# ... #A1#A2#A1)T (4)
since:
(A#B)T = B#A
So why is this useful? From a strictly mathematical standpoint, the new operator is quite useless: it loses the useful and important property of associativity, without producing any new advantages. Numerically, however, it is extremely useful: multiplying a sparse matrix with the transpose of a sparse matrix is much more straightforward than straight multiplication since you move through the non-zero elements of both matrices in order. This is assuming that both are sorted by subscript in row-major order.
Suppose, once again, that v is a matrix. Even if v is sparse, it's usually more efficient to convert it to a full matrix and multiply through using (2) since the result will become full after only a few iterations. Not always: if all the matrices are tri-diagonal, for instance, then the result will also be tri-diagonal.
So lets assume that not only is v a sparse matrix, but that the result will also be sparse. Clearly, (4) is now more efficient than (2) since we now need only 2 matrix transpositions instead of n, where n is the number of A's. Since matrix multiplication is O(n2) (where n is the number of non-zero elements in each matrix, assuming it is about the same for each) while matrix transposition is O(n log n) (for a generalized sparse matrix: swap the subscripts, then re-sort them in row-major order) the gains will be quite modest.
When constructing parallel algorithms, it might be better to stick with (2) since we can multiply through in just about any order we want: e.g., multiply, in order, every second pair of consecutive matrices, then take the results and again multiply by consecutive pairs and so on.
I mention multiplying an array of sparse matrices with a full matrix. This is also an interesting case, although mostly from a numerical perspective. The intuitive approach would be to have two full matrices: one for the result, and one for the multiplicand. We multiply each of the A's in turn, swapping the result with the multiplicand in between. A second approach would be to decompose the multiplicand into columns and multiply each of the columns in turn, all the way through.
It turns out that for large problems, the latter approach is much, much faster, even though it violates at least one optimization principle. If we have two loops, both of which have the same limits, convert them to one loop. For example:
for (int i=0; i < n; i++) a[i]=b[i]*b[i];
for (int i=0; i < n; i++) r+=a[i];
should be converted to:
for (int i=0; i < n; i++) {
a[i]=b[i]*b[i];
r+=a[i];
}
The column-by-column method is better because it uses less memory and because it operates on a smaller chunk of memory for longer. Whereas the first method accesses a whole full matrix at each multiplication step, the second method accesses only one column of the matrix. Therefore there is less paging, either between virtual memory and RAM or between core memory and cache. Like multiplication between two sparse matrices, this method is also more efficient when multiplying with the transpose, assuming row-major storage order once again for both the sparse matrices and the full matrices.
It does well to remind us that while high-level languages like C and C++ do a good job isolating the programmer from the hardware, you don't always end up with the best code by considering the computer as this perfect, abstract machine, instead of a real one with real limitations and real, hacked-together solutions and work-arounds.
It is worth comparing the code for both methods, so I will post it in full. Here is method #1:
template
void sparse_array
real **result) {
real **med;
real **swap;
med=allocate_matrix
sparse_a[0]->mat_mult(cand, result, m);
for (long i=1; i
sparse_a[i]->mat_mult(result, med, m);
swap=med;
med=result;
result=swap;
}
if (nsparse % 2 == 0) {
copy_matrix(result, med, m, m);
}
delete_matrix(med);
}
Here is method #2, applied, however, to the transpose of the multiplicand:
template
void sparse_array::mat_mult_tt(real **cand,
real **result) {
printf("%6.1f%%", 0.);
for (index_t i=0; i
printf("\b\b\b\b\b\b\b%6.1f%%", 100.*i/m);
fflush(stdout);
vect_mult(cand[i], result[i]);
}
printf("\b\b\b\b\b\b\b%6.1f%%\n", 100.);
}
There are a few of things worth noting. In the first method, if the number of sparse matrices is odd, the result and the multiplicand are swapped an odd number of times. For efficiency, only the memory locations are swapped, not the actual values so the memory locations for the multiplicand and the result end up swapped. Therefore we need to copy the values back to the correct location in memory before the function can return. In the second method, we print out a crude progress meter in the form of the percentage of rows that have been multiplied through.
The second version is not only much simpler, it is also easier to parallelize since each row multiplication is completely independent.
No comments:
Post a Comment