# Matrix product verification (Mathematica)

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**(**B**x) and compare it to **C**x for certain randomly chosen vectors *x* can show equality with high probability in only *O*(*n*^{2}) 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*^{2}) Time and log_{2}*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**(**B**x) and compare it to **C**x. 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 |