Rational pseudo-inverse matrix

General Usage
1

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

2

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

3

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)

2 Likes
4

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))
2 Likes
5

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

Would work I think

EDIT: typos in code

1 Like
6

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!

7

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

Doh! You are correct.

8

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
1 Like
9

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

3 Likes
10

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.

For example, if

A=\begin{pmatrix}1&1\\0&0\end{pmatrix},\enspace S=\begin{pmatrix}1&0\\0&1\end{pmatrix}=I,\enspace \Lambda=\begin{pmatrix}1&0\\0&0\end{pmatrix},\enspace T=\begin{pmatrix}1&1\\0&1\end{pmatrix},

then S,T are invertible and A=S\Lambda T, yet

GA=T^{-1}\Lambda^+S^{-1}S\Lambda T=T^{-1}\Lambda^+\Lambda T=\begin{pmatrix}1&1\\0&0\end{pmatrix}

is not Hermitian. (So G cannot be the pseudoinverse.)

However, matrix decompositions are a good idea here — everything works if we just use the rank decomposition instead of the Smith decomposition.

1 Like
11

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

1 Like