Matrix inversion modulo p
Vector 18:2, page 106
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 106
I read Michel Dumontier’s articles Matrix Inversion Modulo p and Inversion of Boolean Matrices in Vol. 18 No. 1 with great interest. However, something did not seem quite right, and after a bit of thought I realized what it was. M. Dumontier is calculating matrix inverses modulo 256, yet he has resorted to multiprecision arithmetic, for which he has had to switch from APL to J. But, surely, when performing calculations modulo 256 you never need to exceed 65536 (i.e. \(256^2\)) in any one value, for which APL should be quite sufficient.
On further investigation, I noticed that the problem is twofold:
- The inverse matrix is calculated by dividing the adjoint (transpose matrix of the cofactors) by the determinant, which is a very inefficient method.
- Numbers should be reduced modulo \(p\) at every opportunity, and not just at the end of the calculation.
The fact that we are working modulo \(p\) should not worry us. When we are working in a field, such as \(\mathbb{R}\), the usual methods involve division somewhere along the way, and this division can be thought of as multiplication by the multiplicative inverse of the divisor. When working modulo \(p\) we also multiply by the multiplicative inverse, but this time modulo \(p\). In all cases if the multiplicative inverse does not exist then the matrix will not have an inverse.
A Better Function¶
The best way of finding an inverse matrix is by using elementary transformations, which are:
- The interchange of two rows (not used in this algorithm)
- The multiplication of a row by a scalar
- The addition of a scalar multiple of one row to another.
The work is carried out on two matrices, one of which is the matrix to be inverted, and the other is the identity matrix. The idea is to reduce our matrix to the identity matrix by applying these elementary transformations, and to end up with the inverse matrix by applying the same transformations to the identity matrix. This works because each of these elementary transformations on the rows of the matrix is equivalent to premultiplying by a suitable matrix (which we do not need to know), giving \(T_1T_2\ldots T_nA = I\) from which we see that \(T_1T_2\ldots T_n = A^{-1}\).
We therefore have the following function:
z←p matinvmod m;i;j;r;x
[1] ⍝ Inverse of matrix m modulo p by elementary transformations
[2] r←1↑⍴m ⍝ Rank
[3] z←(⍴m)⍴(r+1)↑1 ⍝ Start with identity matrix
[4] ⍝ Get the matrix into upper triangular form
[5] :For j :In ⍳r
[6] i←j
[7] ⍝ Find the multiplicative inverse of the jth diagonal element
[8] ⍝ If the inverse is 0, add another row to this one and try again
[9] :While 0=x←m[j;j]Inv_mod p
[10] ⎕SIGNAL(r=i←i+1)/11 ⍝ There isn't one so no inverse mod p
[11] x←(m[j;j]+m[i;j])Inv_mod p
[12] :If 0≠x
[13] m[j;]←m[j;]+m[i;] ⋄ z[j;]←z[j;]+z[i;]
[14] :EndIf
[15] :EndWhile
[16] ⍝ Scalar multiply the jth row of both matrices by x
[17] m[j;]←p|m[j;]×x ⋄ z[j;]←p|z[j;]×x
[18] ⍝ Position [j;j] should now be 1. Subtract row j from the other rows
[19] :For i :In (j+1)↓⍳r
[21] m[i;]←p|m[i;]-m[j;]×x←m[i;j] ⋄ z[i;]←p|z[i;]-z[j;]×x
[22] :EndFor
[23] :EndFor
[24] ⍝ Now work backwards to clear the upper triangle
[25] :For j :In ⌽1↓⍳r
[26] :For i :In ⍳j
[27] z[i;]←p|z[i;]-z[j;]×m[i;j]
[28] ⍝ m[i;]←m[i;]-m[j;]×m[i;j] ⍝ Uncomment for testing
[29] :EndFor
[30] :EndFor
To calculate the multiplicative inverse of a number modulo \(p\) I have used function Inv_mod. This function originally appeared in Vector Vol. 12 No. 4, p.75, in one of a series of articles I wrote on multiprecision arithmetic in APL. Since we are not dealing with multiprecision here, I have removed the multiprecision path from the program, and I have also rewritten it to use control structures, rather than branches and labels. Note that this function returns the impossible value 0 if the inverse does not exist, and this value can be caught and acted upon in the calling function.
z←x Inv_mod n;a;b;N;rn;r0;r1
[1] ⍝ Calculate the multiplicative inverse of x mod n
[2] ⍝ If the inverse does not exist, return 0
[3] a←⍬ ⋄ N←n
[4] z←0 ⋄ x←N|x
[5] :If x=1
[6] z←1
[7] :Else
[8] :While 0≠x
[9] a←a,⊃b←(0,x)⊤n
[10] n←x ⋄ x←1⊃b
[11] :EndWhile
[12] ⍝ At this point n holds the GCD of the original n and x.
[13] ⍝ If this is not 1 then there is no mult. inverse.
[14] :If 1=n
[15] r0←1 0 ⋄ r1←0 1
[16] ⍝ Recurrence relation is r(n)=r(n-2)-q.r(n-1)
[17] ⍝ i.e. rn ← r0 - (⊃a).r1
[18] :While 1<⍴a
[19] rn←r0-r1×⊃a
[20] r0←r1 ⋄ r1←rn ⋄ a←1↓a
[21] :EndWhile
[22] z←N|1⊃rn
[23] :EndIf
[24] :EndIf
The function matinvmod will work very efficiently on the rank-10 matrix shown on p.140 of Vol. 18 No. 1 without the necessity for multiprecision arithmetic. I also applied it to matrix a on p.139, and I’m afraid to say that I discovered that the matrix b shown there is not the inverse of a. (I think this b must have slipped in by accident, because the other calculations imply that b should be the inverse of a as shown below.)
a
33 192 116 135 55 11
173 173 238 97 132 211
8 13 135 171 1 97
17 106 175 150 237 215
134 23 166 106 178 232
194 66 12 187 83 161
256 matinvmod a
79 194 6 146 174 129
30 190 240 162 13 150
246 128 199 74 9 153
73 90 137 61 87 211
11 179 96 240 171 86
186 153 237 33 96 228
256|a+.×256 matinvmod a
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
b
175 50 54 130 174 65
46 166 248 202 13 54
198 72 47 82 9 185
121 146 33 53 87 179
123 139 24 136 171 182
186 153 237 33 96 228
256|a+.×b
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
80 136 40 200 1 32
0 0 0 0 0 1
The Boolean Case¶
This function will also work just as well on boolean matrices, by setting \(p\) to 2. However, in the boolean case there is a fair bit of baggage that can be removed, as follows:
z←boolmatinv m;i;j;r;x
[1] ⍝ Inverse of boolean matrix m by elementary transformations
[2] r←1↑⍴m ⍝ Rank
[3] z←(⍴m)⍴(r+1)↑1 ⍝ Start with identity matrix
[4] ⍝ Get the matrix into upper triangular form
[5] :For j :In ⍳r
[6] i←j
[7] ⍝ Find the multiplicative inverse of the jth diagonal element
[8] ⍝ This will be 1 if the element is 1, make it 0 otherwise
[9] ⍝ If the inverse is 0, add another row to this one and try again
[10] :While 0=x←m[j;j]
[11] ⎕SIGNAL(r=i←i+1)/11 ⍝ There isn't one
[12] x←m[j;j]≠m[i;j]
[13] :If 0≠x
[14] m[j;]←m[j;]≠m[i;] ⋄ z[j;]←z[j;]≠z[i;]
[15] :EndIf
[16] :EndWhile
[17] ⍝ Scalar multiply the jth row of both matrices by x
[18] m[j;]←m[j;]∧x ⋄ z[j;]←z[j;]∧x
[19] ⍝ Position [j;j] should now be 1. Subtract row j from the other rows
[20] :For i :In (j+1)↓⍳r
[21] m[i;]←m[i;]≠m[j;]∧x←m[i;j] ⋄ z[i;]←z[i;]≠z[j;]∧x
[22] :EndFor
[23] :EndFor
[24] ⍝ Now work backwards to clear the upper triangle
[25] :For j :In ⌽1↓⍳r
[26] :For i :In ⍳j
[27] z[i;]←z[i;]≠z[j;]∧m[i;j]
[28] :EndFor
[29] :EndFor