Find the order $m$ of a matrix ${\bf A}$ such that ${\bf A}^m= {\bf 1}$

8

I have a square matrix ${\bf A}$ defined over the field $\mathbb Z_2$ and I want to find its order such that ${\bf A}^m=1$. I tried using

Solve[MatrixPower[A,m] == IdentityMatrix[Length[A]], m, Modulus-> 2]

but this seems to be taking forever. How should I find the order of the matrix as defined? For example, A = {{1, 0}, {1, 1}} has order $2$.

cleanplay

Posted 2018-05-08T20:38:05.210

Reputation: 578

3Could you show a small typical matrix so that people might experiment, confirm what you are seeing and try possible solutions before finding that they are doing something very different from you? – Bill – 2018-05-08T20:46:39.830

1Maybe something like: First@Position[ Table[Mod[MatrixPower[A, n], 2] == IdentityMatrix[2], {n, 4}], True]? Can be improved by stopping computations when True is detected. – anderstood – 2018-05-08T21:06:17.337

2Or: B = A; n = 0; While[(B != IdentityMatrix[2]) && (n < 100), B = Modulus[A.B, 2]; n++]; n – anderstood – 2018-05-08T21:10:22.450

1m = 1; While[m<=100&& Mod[MatrixPower[A, m++], 2] =!= IdentityMatrix[Length[A]]]; m - 1? – kglr – 2018-05-08T21:13:46.377

2also mmax = 1000; Block[{m = 1, B = A}, While[m++ <= mmax && (B = Mod[B.A, 2]) =!= IdentityMatrix[Length[A]]]; m]` – kglr – 2018-05-09T01:01:16.247

Answers

7

Using the example in the OP,

A = {{1, 0}, {1, 1}};

we compare the different options mentioned in the comments (plus one that is mine):

Table[If[Mod[MatrixPower[A, n], 2] == IdentityMatrix[Length@A], Return[n, Table]], {n, 1000}] // RepeatedTiming
(* {0.00007, 2}, by anderstood *)

Block[{m = 1}, While[m <= 100 && Mod[MatrixPower[A, m++], 2] =!= IdentityMatrix[Length[A]]]; m - 1] // RepeatedTiming
(* {0.00003061, 2}, by kglr *)

Length[NestWhileList[Mod[#.A, 2] &, A, Mod[#, 2] =!= IdentityMatrix[Length[A]] &]] // RepeatedTiming
(* {9.8*10^-6, 2}, by me *)

Block[{mmax = 1000, m = 1, B = A}, While[m++ <= mmax && (B = Mod[B.A, 2]) =!= IdentityMatrix[Length[A]]]; m] // RepeatedTiming
(* {5.5*10^-6, 2}, by kglr *)

For less trivial examples, consider

SeedRandom[2];
A = RandomInteger[{0, 1}, {5, 5}];
(* {0.001, 12}, by anderstood *)
(* {0.00088, 12}, by kglr *)
(* {0.000064, 12}, by me *)
(* {0.000046, 12}, by kglr *)

SeedRandom[2];
A = RandomInteger[{0, 1}, {10, 10}];
(* {0.04, 60}, by anderstood *)
(* {0.037, 60}, by kglr *)
(* {0.00051, 60}, by me *)
(* {0.00037, 60}, by kglr *)

AccidentalFourierTransform

Posted 2018-05-08T20:38:05.210

Reputation: 9 042

Changing Modulus to Mod gives basically the same answer as kglr number 2, and the timing is very similar: Block[{}, B = A; n = 0; While[(B =!= IdentityMatrix[5]) && (n < 1000), B = Mod[A.B, 2]; n++]; n + 1]. And here is an updated version of my first comment with a break: Table[If[Mod[MatrixPower[A, n], 2] == IdentityMatrix[Length@A], Return[n, Table]], {n, 1000}] // RepeatedTiming. – anderstood – 2018-05-09T15:20:03.603

7

The Code

I solve the more generic question of finding the power of a matrix under any mod-z with the following command:

Options[solver] = {matrixForm -> False};
solver[input_, mod_, iteration_, opts : OptionsPattern[]] := 
Module[{a, exp, exp2},
a /: Power[a, n_] := a;
a /: Times[a, n_] /; (n > mod - 1) := Mod[n, mod] a;
exp = Simplify[a input];
exp2 = NestWhileList[{#[[1]] + 1, exp.#[[2]]} &, {1, exp}, 
 And[(! Simplify[#[[2]] === a IdentityMatrix[Length@input]]), #[[
     1]] < iteration] &];
If[OptionValue[matrixForm] == True,
Grid[{#[[1]], #[[2]] // MatrixForm} & /@ (exp2 /. a -> 1)],
exp2 /. a -> 1
]
];

The code basically multiplies the given matrix with a scooping variable which is constrained such that its coefficients are in modulo-chosen number (second argument). Then we calculate its matrix powers till we find the identity matrix or we reach the chosen iteration number (third argument). We also have the option to make the output more visual.

Check of OP's matrix

solver[{{1, 0}, {1, 1}}, 2, 100] // Timing

{0.000244, {{1, {{1, 0}, {1, 1}}}, {2, {{1, 0}, {0, 1}}}}}

where we get a list of matrices $\{A,A^2,A^3,\cdots\}$. We can make it more visual as follows:

solver[{{1, 0}, {1, 1}}, 2, 100, matrixForm -> True]

$$\begin{array}{cc} 1 & \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right) \\ 2 & \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right) \\ \end{array}$$

We can also check modulo other numbers:

solver[{{1, 0}, {1, 1}}, 3, 100, matrixForm -> True] // Timing

$$\left\{0.000352, \begin{array}{cc} 1 & \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ \end{array} \right) \\ 2 & \left( \begin{array}{cc} 1 & 0 \\ 2 & 1 \\ \end{array} \right) \\ 3 & \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right) \\ \end{array} \right\}$$

Another example

A = RandomInteger[{0, 1}, {5, 5}]

{{1, 0, 0, 0, 1}, {0, 1, 0, 1, 0}, {0, 1, 0, 0, 1}, {0, 1, 1, 1, 0}, {0, 0, 0, 1, 0}}

solver[A, 2, 100, matrixForm -> True] // Timing

$$\left\{0.001191, \begin{array}{cc} 1 & \left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ \end{array} \right) \\ 2 & \left( \begin{array}{ccccc} 1 & 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \\ \end{array} \right) \\ 3 & \left( \begin{array}{ccccc} 1 & 1 & 1 & 0 & 1 \\ 0 & 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 1 & 1 & 0 & 1 \\ \end{array} \right) \\ 4 & \left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) \\ \end{array} \right\}$$

For the same matrix, we can find the necessary power in module-3 as well:

solver[A, 3, 100] // Timing

{0.006589, {{1, {{1, 0, 0, 0, 1}, {0, 1, 0, 1, 0}, {0, 1, 0, 0, 1}, {0, 1, 1, 1, 0}, {0, 0, 0, 1, 0}}}, {2, {{1, 0, 0, 1, 1}, {0, 2, 1, 2, 0}, {0, 1, 0, 2, 0}, {0, 0, 1, 2, 1}, {0, 1, 1, 1, 0}}}, {3, {{1, 1, 1, 2, 1}, {0, 2, 2, 1, 1}, {0, 0, 2, 0, 0}, {0, 0, 2, 0, 1}, {0, 0, 1, 2, 1}}}, {4, {{1, 1, 2, 1, 2}, {0, 2, 1, 1, 2}, {0, 2, 0, 0, 2}, {0, 2, 0, 1, 2}, {0, 0, 2, 0, 1}}}, {5, {{1, 1, 1, 1, 0}, {0, 1, 1, 2, 1}, {0, 2, 0, 1, 0}, {0, 0, 1, 2, 0}, {0, 2, 0, 1, 2}}}, {6, {{1, 0, 1, 2, 2}, {0, 1, 2, 1, 1}, {0, 0, 1, 0, 0}, {0, 0, 2, 2, 1}, {0, 0, 1, 2, 0}}}, {7, {{1, 0, 2, 1, 2}, {0, 1, 1, 0, 2}, {0, 1, 0, 0, 1}, {0, 1, 2, 0, 2}, {0, 0, 2, 2, 1}}}, {8, {{1, 0, 1, 0, 0}, {0, 2, 0, 0, 1}, {0, 1, 0, 2, 0}, {0, 0, 0, 0, 2}, {0, 1, 2, 0, 2}}}, {9, {{1, 1, 0, 0, 2}, {0, 2, 0, 0, 0}, {0, 0, 2, 0, 0}, {0, 0, 0, 2, 0}, {0, 0, 0, 0, 2}}}, {10, {{1, 1, 0, 0, 1}, {0, 2, 0, 2, 0}, {0, 2, 0, 0, 2}, {0, 2, 2, 2, 0}, {0, 0, 0, 2, 0}}}, {11, {{1, 1, 0, 2, 1}, {0, 1, 2, 1, 0}, {0, 2, 0, 1, 0}, {0, 0, 2, 1, 2}, {0, 2, 2, 2, 0}}}, {12, {{1, 0, 2, 1, 1}, {0, 1, 1, 2, 2}, {0, 0, 1, 0, 0}, {0, 0, 1, 0, 2}, {0, 0, 2, 1, 2}}}, {13, {{1, 0, 1, 2, 0}, {0, 1, 2, 2, 1}, {0, 1, 0, 0, 1}, {0, 1, 0, 2, 1}, {0, 0, 1, 0, 2}}}, {14, {{1, 0, 2, 2, 2}, {0, 2, 2, 1, 2}, {0, 1, 0, 2, 0}, {0, 0, 2, 1, 0}, {0, 1, 0, 2, 1}}}, {15, {{1, 1, 2, 1, 0}, {0, 2, 1, 2, 2}, {0, 0, 2, 0, 0}, {0, 0, 1, 1, 2}, {0, 0, 2, 1, 0}}}, {16, {{1, 1, 1, 2, 0}, {0, 2, 2, 0, 1}, {0, 2, 0, 0, 2}, {0, 2, 1, 0, 1}, {0, 0, 1, 1, 2}}}, {17, {{1, 1, 2, 0, 2}, {0, 1, 0, 0, 2}, {0, 2, 0, 1, 0}, {0, 0, 0, 0, 1}, {0, 2, 1, 0, 1}}}, {18, {{1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 0, 1, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 1}}}}}

where we see that it is 18th power, and we get the result in around 6 ms, which is fast enough.

Another example

A = RandomInteger[{0, 4}, {4, 4}]

{{2, 1, 1, 1}, {1, 2, 0, 2}, {4, 3, 2, 2}, {1, 4, 3, 2}}

With[{result = solver[A, 5, 100] // Timing}, {result[[1]],Last[result[[2]]]}]

$$\left\{0.02803,\left\{62,\left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)\right\}\right\}$$

We can compare the result with direct computation:

Mod[#, 5] & /@ Nest[#.A &, A, 61]

{{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}

which proves that indeed 62nd power of the matrix A above is identity matrix for modulo-5 integers.

Soner

Posted 2018-05-08T20:38:05.210

Reputation: 1 903

1For such short times, I think it's better to use RepeatedTiming instead of Timing (not that it's important to get an accurate timing here, just saying...) – AccidentalFourierTransform – 2018-05-09T18:30:51.057

2@AccidentalFourierTransform Thanks! – Soner – 2018-05-10T02:26:10.947

5

Another implementation of Soner's solver:

cf = Compile[{{A, _Integer, 2}, {B, _Integer, 2}, {mod, _Integer, 0}},
             Mod[A.B, mod], Parallelization -> True]
solver[A_, mod_, maxIter_: 1000] := 
 Block[{id = IdentityMatrix[Length@A], B = A, n = 0}, 
         While[(B =!= id) && (n < 1000), B = cf[A, B, mod]; n++]; n + 1]

Examples:

A = {{2, 1, 1, 1}, {1, 2, 0, 2}, {4, 3, 2, 2}, {1, 4, 3, 2}};
solver[A, 5] // RepeatedTiming
(* {0.00022, 62} *)

A = {{1, 0, 0, 0, 1}, {0, 1, 0, 1, 0}, {0, 1, 0, 0, 1}, {0, 1, 1, 1, 0},
     {0, 0, 0, 1, 0}};
solver[A, 2] // RepeatedTiming
(* {0.000018, 4} *)

solver[A, 3] // RepeatedTiming
(* {0.000074, 18} *)

anderstood

Posted 2018-05-08T20:38:05.210

Reputation: 12 994

1@AccidentalFourierTransform You're right. And with RepeatedTiming the time is greatly reduced; there must be some overhead. – anderstood – 2018-05-09T18:45:26.850

3

MatrixPowerMod from A Course in Computational Number Theory by Bressoud and Wagon, page 34, is quite fast.

MatrixPowerMod[a_, m_, p_] :=
   Block[{b = a, d = IntegerDigits[m, 2]},
      Do[
         b = Mod[b.b, p];
         If[d[[i]] == 1, b = Mod[b.a, p]],
         {i, 2, Length[d]}];
      b]

Then use some information from an answer to this question by loup blanc.

MatrixOrderMod[a_, p_] :=
   Block[{d, s, k, ii = IdentityMatrix[Length[a]]},
         d = Mod[Det[a], p];
         If[d == 0, 0,
            s = MultiplicativeOrder[d, p];
            k = s;
            While[MatrixPowerMod[a, k, p] =!= ii, k += s];
            k]
   ]

MatrixOrderMod is a brute force search, as others have implemented, but is faster in general because the matrix exponent is incremented by s instead of by 1. The code also tests for zero determinant, which corresponds to no solution. This test eliminates the need for an explicit upper bound on the exponent.

KennyColnago

Posted 2018-05-08T20:38:05.210

Reputation: 14 269

1Nice. I thought about the determinant test too, but I wasn't sure: there is a solution if and only if $\det(A)\neq 0$? the only if direction is obvious, but the other one -- not quite so, at least to me. – AccidentalFourierTransform – 2018-05-09T23:31:02.710

This can be compacted a bit: MatrixOrderMod[a_, p_] := Block[{k = 0, ii = IdentityMatrix[Length[a]], d, s}, d = Det[a, Modulus -> p]; If[d == 0, 0, s = MultiplicativeOrder[d, p]; While[Algebra`MatrixPowerMod[a, k += s, p] =!= ii]; k]]. – J. M.'s ennui – 2018-10-01T20:32:11.410