Matrix product verification (Mathematica)

From LiteratePrograms

A frequently used introductory example of randomized algorithms is matrix product verification, where we are given two matrices A and B and a matrix C and asked whether AB = C. With currently known matrix multiplication algorithms, the obvious approach of multiplying out AB takes Ω(n^2.37) time. Simple randomized algorithms that instead compute A(Bx) and compare it to Cx for certain randomly chosen vectors x can show equality with high probability in only O(n²) time.

These implementations are based primarily on R. Frievald's 1979 paper "Fast probabilistic algorithms" as surveyed in Tracy Kimbel and Rakesh Sinha's A Probabilistic Algorithm for Verifying Matrix Products Using O(n²) Time and log₂n + O(1) Random Bits.

In Mathematica, which already has built in support for matrix multiplication, loops, and choosing random integers, it's straightforward to implement Frievald's original algorithm. We first choose a vector with as many rows as B has columns, selecting each element randomly as either -1 and 1 with 50% chance each:

<<compute random vector>>=
x = Table[Random[Integer,{0,1}]*2 - 1, {i, 1, Length[B[[1]]]}]

Note that our use of pseudorandom values rather than actual random values for this vector could potentially make the test fail in certain very carefully chosen cases; for the test cases we consider it seems to work fine.

Next, we compute A(Bx) and compare it to Cx. If these are unequal, then AB ≠ C and we return False. Otherwise, AB = C with probability one-half, and we return True:

<<single vector product verification>>=
ProbMatrixProductTest[A_,B_,C_] :=
    Block[{x}, compute random vector;
               Return[A.(B.x) == C.x]]

To be accurate with probability close to 1, we must repeat the test a number of times. We do this using Mathematica's ForAll operator, which will yield true only if the given expression evaluates to true for every argument in the list in its first argument:

<<multiple vector product verification>>=
ProbMatrixProductVerify[A_,B_,C_,k_] :=
    ForAll[Range[k], ProbMatrixProductTest[A, B, C]]

We demonstrate with some randomly chosen matrices. We use small integer matrix entries to avoid false negatives due to numerical rounding errors, and replace C with M to avoid a conflict with the built-in symbol:

<<create test matrices>>=
A = Table[Random[Integer, {1, 100}], {i, 1, 10}, {j, 1, 10}];
B = Table[Random[Integer, {1, 100}], {i, 1, 10}, {j, 1, 10}];
M = A.B;
Map[MatrixForm, {A, B, M}]

Both of the above functions will always return True when run on A, B, and M. If we slightly perturb one entry of M, we'll find that the test will consistently return False for sufficiently large values of k:

<<test positive and negative verification>>=
Print[ProbMatrixProductVerify[A, B, M, 20]];
i = Random[Integer,{1,10}]; j = Random[Integer,{1,10}];
M[[i, j]] = M[[i, j]] + 1;
Print[ProbMatrixProductVerify[A, B, M, 20]];

This will print True followed by False (with very high probability). We place this implementation in a notebook file for reuse:

<<MatrixProductVerification.nb>>=
single vector product verification
multiple vector product verification
create test matrices
test positive and negative verification

Download code