Rational pseudo-inverse matrix

Given a generic non-square matrix A, one can compute the pseudo-inverse (aka Moore-Penrose inverse) via the LinearAlgebra.pinv function, or also from the LinearAlgebra.qr decomposition (pinv(A) = inv(qr(A).R) * qr(A).Q'). The result has floating-point eltype (it is computed with LAPACK I think). However, if A::AbstractMatrix{<:Integer}, I think one can prove that pinv(A) is a matrix of rational numbers, so that n * pinv(A) is a matrix of integers for some integer n.

Question: how can one find an n for a generic non-square A::AbstractMatrix{<:Integer}?

If A is a square invertible matrix, it is easy: n = det(A) (can be read off the explicit formula for the inverse in this case). I’m struggling to find the equivalent for a general A. Any ideas??

EDIT: typo corrected

The pseudo inverse is (a'a) \a', can it be decomposed into a subproblem using the known solution to the inverse of a'a?

Yes! Thanks @baggepinnen! n = det(a'a) is a valid n. It is not necessarily the minimum n, but it is always a good solution.

(This community is truly fantastic)


Its unclear to me if you need this proof to get around the floating-point eltype of pinv or for some bigger reason. If its to get around the eltype, you can convert the Ints to Rationals to prevent promotion to floats

using LinearAlgebra

A = rand(1:10, 5,5)
B = Rational.(A)

# promotes to float

# Rational
pinv_B = (B'*B)\B'
n = maximum(denominator.(pinv_B))

I think that in order to find the smallest n such that n*pinv(A) is a matrix of integers, one would have to find the least common multiple of the denominators of pinv(A). Luckily, the base function lcm can find the least common multiple of a series of numbers.

So something like:

function npinv(A)
    rA = Rational.(A)
    B = (transpose(rA)*rA)\transpose(rA)
    n = lcm([denominator(e) for e in B])

    return n, B

Would work I think

EDIT: typos in code

1 Like

Nice @DrPapa! That’s another way, avoiding pinv altogether and using (B'*B)\B'. However, I suspect that is going to be much slower than det(a'a), right? Also, it is not enough to take the maximum denominator, one would need the least common multiple of all denominators, I’d guess

EDIT: sniped by @Bruno_Amorim!

You are most likely correct. Either way, converting to Rational should work.

Doh! You are correct.

I think I’ll mark Bruno’s post as the solution, since it gives not only a valid n but the minimum n.

In my usecase I use StaticArrays as input, so it becomes important to avoid allocations. I leave here an uglier variation of @baggepinnen, @Bruno_Amorim and @DrPapa’s idea that manages to work without allocations and Rationals when fed a static matrix. It might perhaps be cleaned even further, but its already in the 100ns runtime range for a 3x2 matrix, so good enough for me…

function npinv(A)
    n = Int(det(A'A))
    qrA = qr(A)
    iA = inv(qrA.R) * qrA.Q'
    iAn = round.(Int, n .* iA)
    d = gcd(iAn) # common divisor to remove from iAn
    n´ = div(n, d, RoundNearest)
    iAn´ = div.(iAn, d, RoundNearest)
    return n´, iAn´
1 Like

Bit late to the party, but wanted to note that this could be answered neatly and stably (the lcm approach seems susceptible to floating point errors?) via the Smith normal decomposition (there are Julia implementations in Nemo.jl I believe and also at wildart/SmithNormalForm.jl). Specifically, the n above is the maximum of the invariant factors (or elementary divisors) of A, i.e. the maximum of the Smith normal form diagonal.

This is because the Smith decomposition A = S\Lambda T can be used to define the pseudoinverse of A, namely A^+ = T^{-1}\Lambda^+S^{-1}. Since S and T are integer unimodular matrices, as are their inverses - and the only non-integer bits are therefore those in the diagonal matrix \Lambda^+ (and they are just the inverses of the nonzero diagonal of \Lambda, i.e. the invariant factors of A - which are also just integers).