# Rational pseudo-inverse matrix

**URL:** <https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773>\
**Category:** General Usage\
**Created:** [July 27, 2020, 4:54pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773 "2020-07-27T16:54:42Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![pablosanjose](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablosanjose/32/7006_2.png) [@pablosanjose](https://discourse.julialang.org/u/pablosanjose)\
**Post date:** [July 27, 2020, 4:54pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/1 "2020-07-27T16:54:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [July 27, 2020, 5:12pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/2 "2020-07-27T17:12:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![pablosanjose](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablosanjose/32/7006_2.png) [@pablosanjose](https://discourse.julialang.org/u/pablosanjose)\
**Post date:** [July 27, 2020, 5:18pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/3 "2020-07-27T17:18:01Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [July 27, 2020, 6:05pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/4 "2020-07-27T18:05:01Z")

</div>

> [@pablosanjose](#):
>
> 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` .

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

```julia
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))

```

---

<div class="post-metadata">

**Author:** ![Bruno\_Amorim](https://avatars.discourse-cdn.com/v4/letter/b/f6c823/32.png) [@Bruno\_Amorim](https://discourse.julialang.org/u/Bruno_Amorim)\
**Post date:** [July 27, 2020, 6:17pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/5 "2020-07-27T18:17:20Z")

</div>

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:

```julia
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

---

<div class="post-metadata">

**Author:** ![pablosanjose](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablosanjose/32/7006_2.png) [@pablosanjose](https://discourse.julialang.org/u/pablosanjose)\
**Post date:** [July 27, 2020, 6:17pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/6 "2020-07-27T18:17:29Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [July 27, 2020, 6:37pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/7 "2020-07-27T18:37:14Z")

</div>

> [@pablosanjose](#):
>
> However, I suspect that is going to be much slower than `det(a'a)` , right?

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

> [@pablosanjose](#):
>
> Also, it is not enough to take the maximum denominator, one would need the least common multiple of all denominators

Doh! You are correct.

---

<div class="post-metadata">

**Author:** ![pablosanjose](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablosanjose/32/7006_2.png) [@pablosanjose](https://discourse.julialang.org/u/pablosanjose)\
**Post date:** [July 27, 2020, 6:53pm UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/8 "2020-07-27T18:53:09Z")

</div>

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…

```julia
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

```

---

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [April 25, 2021, 3:20am UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/9 "2021-04-25T03:20:43Z")

</div>

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](https://en.m.wikipedia.org/wiki/Smith_normal_form) (there are Julia implementations in Nemo.jl I believe and also at [wildart/SmithNormalForm.jl](https://github.com/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).

---

<div class="post-metadata">

**Author:** ![thatsme314](https://avatars.discourse-cdn.com/v4/letter/t/dfb087/32.png) [@thatsme314](https://discourse.julialang.org/u/thatsme314)\
**Post date:** [September 16, 2022, 5:58am UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/10 "2022-09-16T05:58:26Z")

</div>

Careful, the “identity” (BC)^+=C^+B^+ does not always hold ([Wikipedia](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse#Products)), 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](https://en.wikipedia.org/wiki/Rank_factorization) instead of the Smith decomposition.

---

<div class="post-metadata">

**Author:** ![thatsme314](https://avatars.discourse-cdn.com/v4/letter/t/dfb087/32.png) [@thatsme314](https://discourse.julialang.org/u/thatsme314)\
**Post date:** [September 16, 2022, 6:53am UTC](https://discourse.julialang.org/t/rational-pseudo-inverse-matrix/43773/11 "2022-09-16T06:53:56Z")

</div>

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](https://math.stackexchange.com/questions/854578/some-questions-about-the-pseudoinverse-of-a-matrix/854773).)