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

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)

Julia Community 🟣

Integral matrices with integral inverses using Sherman-Morrison formula

Top comments (0)