## Julia Community ðŸŸ£ is a community of amazing users

A fresh approach to technical computing.

Vinod V

Posted on

# Integral matrices with integral inverses using Sherman-Morrison formula

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$ .

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)