Let me explain how to make a square matrix with integer entries so that its inverse contains only integer entries. This kind of matrices might be useful for demonstrating Gaussian-Jordan elimination for solving a linear system without any fractions coming up in between.

The Sherman-Morrison formula is given by

$(B + uv^T)^{-1} = B^{-1} - \frac{B^{-1} uv^T B^{-1}}{1 + v^T B^{-1} u} ... (1)$
,

where
$B$
is an invertible matrix of order
$n$
and
$u$
and
$v$
are column vectors of size
$n$
.

Ref : Sherman–Morrison formula - Wikipedia

Let
$B = I$
and
$u^T v = 0$
in
$(1)$
.The vectors
$u$
and
$v$
are **orthogonal**.

If $A = I + uv^T$ , then $A^{-1} = I - uv^T$ . We can notice that if entries of both $u$ and $v$ are integers then entries of both $A$ and $A^{-1}$ will be integers. This could be useful for designing a consistent linear system $Ax = b$ with integer coefficient matrix. If the constant vector $b$ contains only integer entries, then the solution $x$ also contains only integer entries. Integer matrices with integer inverses are called unimodular matrices. The determinant of $A$ will be 1.

Let $u = \begin{bmatrix} 1 \\ -4 \\ 1 \\ 4 \end{bmatrix}$

and $v = \begin{bmatrix} -4 \\ -4 \\ -4 \\ -2 \end{bmatrix}$ . It can be verified that $u^Tv = 0$ .

Then $A = I + uv^T = \begin{bmatrix} -3 & -4 & -4 & -2 \\ 16 & 7 & 6 & 8 \\ -4 & -4 & -3 & -2 \\ -16 & -16 & -16 & -7 \end{bmatrix}$

The resulting $4 \times 4$ linear system $Ax = b$ can be written as (with $b = \begin{bmatrix} -1 \\ -6 \\ 8 \\ -5 \end{bmatrix}$ )

$-3x_1 - 4x_2 - 4x_3 - 2x_4 = -1$

$16x_1 + 7x_2 + 6x_3 + 8x_4 = 6$

$-4x_1 - 4x_2 - 3x_3 - 2x_4 = 8$

$-16x_1 - 16x_2 - 16x_3 - 7x_4 = -5$

Now the Gauss-Jordan elimination can be applied to the above linear system to get a unique solution as $x_1 = -161,x_2 = 314,x_3 = 8,x_4 = -165$

I have coded to make unimodular matrices in Julia. The flowchart of the program is given below.

```
using LinearAlgebra
using Primes
function getVectors(n,lim)
#get the first n
global u,v,p
while true
u = rand(-lim:lim,n,1)
v = rand(-lim:lim,n,1)
println("u = ",u)
println("v = ",v)
p = u'*v
if 0 == p
continue
end
if abs(p[1]) > 0
break
end
end
return u,v,p
end
function getCoefficientMatrix(n,lim)
#Return two matrices: A and its inverse Ainv
#Coefficient matrix A such that A = I + uv^T
#Ainv = I - uv^T
#lim = 5
#get random column vectors u and v whose elements are in the range [-lim,lim]
#p is the inner product of u and v, p = u^T v.
u,v,p = getVectors(n,lim)
println("p = u^T v = ",p[1])
#Add one more element to u and v such that the resulting vectors are orthogonal.
#Factorize p[1] into two numbers factor1 and factor2
mult = 1
if p[1] > 0
mult = -1
end
factors = factor(Vector,p[1]) #from the Primes package
println("Factors of p are ")
println(factors)
switches = rand([true,false],length(factors))
if length(factors) > 1
switches[1:2] .= [true,false]
end
factor1 = 1
factor2 = 1
k = 1
for f in factors
if true == switches[k]
factor1 *= f
else
factor2 *= f
end
k = k + 1
end
factor1 *= mult
println("factor1 = ",factor1)
println("factor2 = ",factor2)
#Append last element to u and v
u = vcat(u,factor1)
v = vcat(v,factor2)
#Now u and v are orthogonal
println("u = ",u)
println("v = ",v)
println("u'*v = ",u'*v)
A = I + u*v'
Ainv = I - u*v'
return A,Ainv
end
function getLatexExprLinearSystem(A,b)
m,n = size(A)
stri = "\$"
for i in 1:m
stri *= join([string(A[i,j])*"x_"*string(j) for j in 1:n], " + ")
stri *= " = " * string(b[i])
if m == i
stri *= "\$\n"
else
stri *= "\$\n\$"
end
end
stri = replace(stri,"+ -"=>" - ")
stri = replace(stri,"+ 1" => "+ ")
stri = replace(stri,".0x" => "x")
stri
end
n = 4 #Size of the linear system
limit = 4 #Approximate range of integers in the vectors u and v
A,Ainv = getCoefficientMatrix(n-1,limit)
println("The coefficient matrix A is")
display(A)
println("The inverse of coefficient matrix A is")
display(Ainv)
b = rand(-10:10,n) #the constant vector b of Ax = b
s = getLatexExprLinearSystem(A,b)
println("Latex expression for the linear system is :")
print(s)
println("The solution for the linear system is :")
display(Ainv*b)
```

## Discussion (0)