The block Wiedemann algorithm for computing kernel vectors of a matrix over a finite field is a generalisation of an algorithm due to Don Coppersmith.
Let
M
be an
n
×
n
square matrix over some finite field F, let
x
b
a
s
e
be a random vector of length
n
, and let
x
=
M
x
b
a
s
e
. Consider the sequence of vectors
S
=
[
x
,
M
x
,
M
2
x
,
…
]
obtained by repeatedly multiplying the vector by the matrix
M
; let
y
be any other vector of length
n
, and consider the sequence of finite-field elements
S
y
=
[
y
⋅
x
,
y
⋅
M
x
,
y
⋅
M
2
x
…
]
We know that the matrix
M
has a minimal polynomial; by the Cayley–Hamilton theorem we know that this polynomial is of degree (which we will call
n
0
) no more than
n
. Say
∑
r
=
0
n
0
p
r
M
r
=
0
. Then
∑
r
=
0
n
0
y
⋅
(
p
r
(
M
r
x
)
)
=
0
; so the minimal polynomial of the matrix annihilates the sequence
S
and hence
S
y
.
But the Berlekamp–Massey algorithm allows us to calculate relatively efficiently some sequence
q
0
…
q
L
with
∑
i
=
0
L
q
i
S
y
[
i
+
r
]
=
0
∀
r
. Our hope is that this sequence, which by construction annihilates
y
⋅
S
, actually annihilates
S
; so we have
∑
i
=
0
L
q
i
M
i
x
=
0
. We then take advantage of the initial definition of
x
to say
M
∑
i
=
0
L
q
i
M
i
x
b
a
s
e
=
0
and so
∑
i
=
0
L
q
i
M
i
x
b
a
s
e
is a hopefully non-zero kernel vector of
M
.
The natural implementation of sparse matrix arithmetic on a computer makes it easy to compute the sequence S in parallel for a number of vectors equal to the width of a machine word – indeed, it will normally take no longer to compute for that many vectors than for one. If you have several processors, you can compute the sequence S for a different set of random vectors in parallel on all the computers.
It turns out, by a generalization of the Berlekamp–Massey algorithm to provide a sequence of small matrices, that you can take the sequence produced for a large number of vectors and generate a kernel vector of the original large matrix. You need to compute
y
i
⋅
M
t
x
j
for some
i
=
0
…
i
max
,
j
=
0
…
j
max
,
t
=
0
…
t
max
where
i
max
,
j
max
,
t
max
need to satisfy
t
max
>
d
i
max
+
d
j
max
+
O
(
1
)
and
y
i
are a series of vectors of length n; but in practice you can take
y
i
as a sequence of unit vectors and simply write out the first
i
max
entries in your vectors at each time t.
