Commutant of set of matrices

I'm trying to compute the commutant of a set of matrices, so the set of matrices that commutes with all matrices in the original set.

I'm a complete newbie with mathematica, so I would be greatful for help.

In principle I thought it should be possible to define a matrix of variables M and then use M.B-B.M==0 for every B in my set of matrices to build a large set of equations which mathematica should be able to solve to obtain a basis for M. However, I am not quite sure how to execute it.

Thanks!

Take a random matrix:

n = 3;
M = RandomInteger[{-10, 10}, {n, n}]
(*    {{9, 1, 8},
{-7, 8, -3},
{7, 10, -7}}    *)


Build a matrix of unknowns:

B = Array[b, {n, n}]
vars = Flatten[B]

(*    {{b[1, 1], b[1, 2], b[1, 3]},
{b[2, 1], b[2, 2], b[2, 3]},
{b[3, 1], b[3, 2], b[3, 3]}}    *)

(*    {b[1, 1], b[1, 2], b[1, 3], b[2, 1], b[2, 2], b[2, 3], b[3, 1], b[3, 2], b[3, 3]}    *)


Build the equations for $$B$$ in matrix form, and look for the null space:

s = NullSpace[D[Flatten[M . B - B . M], {vars}]]
(*    {{1, 0, 0, 0, 1, 0, 0, 0, 1},
{183, 105, 77, -196, 72, -83, 0, 97, 0},
{-278, -953, 6, 1281, 735, 539, 679, 0, 0}}    *)


Arrange these solutions into $$n\times n$$ matrices:

q = Partition[#, n] & /@ s
(*    {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
{{183, 105, 77}, {-196,  72, -83}, {0, 97, 0}},
{{-278, -953, 6}, {1281, 735, 539}, {679, 0, 0}}}    *)


Try if they commute with $$M$$:

M . # - # . M & /@ q
(*    {{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}}    *)


Build an arbitrary matrix that commutes with $$M$$:

qq = {α, β, γ} . q
(*    {{α + 183 β - 278 γ, 105 β - 953 γ, 77 β + 6 γ},
{-196 β + 1281 γ, α + 72 β + 735 γ, -83 β + 539 γ},
{679 γ, 97 β, α}}                                      *)

M . qq - qq . M // Expand
(*    {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}    *)


Alternatively, go straight for the Solve way:

Solve[Thread[Flatten[M . B - B . M] == 0]]
(*    {{b[2, 1] -> -(147/109) b[1, 2] - 77/109 b[1, 3],
b[2, 2] -> b[1, 1] - 811/763 b[1, 2] + 6/763 b[1, 3],
b[2, 3] -> -(433/763) b[1, 2] - 232/763 b[1, 3],
b[3, 1] -> -(77/109) b[1, 2] + 105/109 b[1, 3],
b[3, 2] -> 6/763 b[1, 2] + 953/763 b[1, 3],
b[3, 3] -> b[1, 1] - 232/763 b[1, 2] - 1497/763 b[1, 3]}}    *)