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??

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
pinv(A)
pinv(B)
(A'*A)\A'
# 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
end

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

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Â´
end

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

Careful, the â€śidentityâ€ť (BC)^+=C^+B^+ does not always hold (Wikipedia), so we canâ€™t conclude A^+=T^{-1}\Lambda^{-1}S^{-1}=:G. It is true that we have GAG=A and AGA=G, so G is a (not necessarily unique) reflexive generalized inverse of A. However, we may possibly fail to satisfy the Hermitianness axioms (GA)^*=GA, (AG)^*=AG, and thus G might not be the pseudoinverse.

Also, regarding the previously mentioned n, I wanted to add that the previously mentioned \det(A^\top A) (a.k.a. product of singular values) is indeed a valid n when A^\top A has full column rank. If A^\top A does not have full column rank, however, then \det(A^\top A)=0, which you might agree is not what weâ€™re looking for â€” we want something strictly positive! However, in this case, we can still take the product of the nonzero singular values. Denoting this quantity (the product of the nonzero singular values) by \text{volvol}(A) (â€śsquared volumeâ€ť), we can show that we again have n\mid\text{volvol}(A).

(This post and the above post are a slightly modified version of a discussion at Math.StackExchange.)