Julia Community 🟣

Vinod V
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+uvT)βˆ’1=Bβˆ’1βˆ’Bβˆ’1uvTBβˆ’11+vTBβˆ’1u...(1)(B + uv^T)^{-1} = B^{-1} - \frac{B^{-1} uv^T B^{-1}}{1 + v^T B^{-1} u} ... (1) ,
where BB is an invertible matrix of order nn and uu and vv are column vectors of size nn .

Ref : Sherman–Morrison formula - Wikipedia

Let B=IB = I and uTv=0u^T v = 0 in (1)(1) .The vectors uu and vv are orthogonal.

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

Let u=[1βˆ’414]u = \begin{bmatrix} 1 \\ -4 \\ 1 \\ 4 \end{bmatrix}

and v=[βˆ’4βˆ’4βˆ’4βˆ’2]v = \begin{bmatrix} -4 \\ -4 \\ -4 \\ -2 \end{bmatrix} . It can be verified that uTv=0u^Tv = 0 .

Then A=I+uvT=[βˆ’3βˆ’4βˆ’4βˆ’216768βˆ’4βˆ’4βˆ’3βˆ’2βˆ’16βˆ’16βˆ’16βˆ’7]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Γ—44 \times 4 linear system Ax=bAx = b can be written as (with b=[βˆ’1βˆ’68βˆ’5]b = \begin{bmatrix} -1 \\ -6 \\ 8 \\ -5 \end{bmatrix} )

βˆ’3x1βˆ’4x2βˆ’4x3βˆ’2x4=βˆ’1-3x_1 - 4x_2 - 4x_3 - 2x_4 = -1
16x1+7x2+6x3+8x4=616x_1 + 7x_2 + 6x_3 + 8x_4 = 6
βˆ’4x1βˆ’4x2βˆ’3x3βˆ’2x4=8-4x_1 - 4x_2 - 3x_3 - 2x_4 = 8
βˆ’16x1βˆ’16x2βˆ’16x3βˆ’7x4=βˆ’5-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 x1=βˆ’161,x2=314,x3=8,x4=βˆ’165x_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.

Flowchart to make unimodular matrices

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)
Enter fullscreen mode Exit fullscreen mode

Top comments (0)