OK… I’ve made a function called `rrqr`

– it is meant to stand for “rank revealing QR” decomposition. (I.e., I *think* my code produces what is known as rank revealing, although I don’t have time to prove it…). Currently, it only works with `Rational`

type matrices. The code is pretty rotten and probably highly inefficient, so if anyone is interested in (i) cleaning up the code, (ii) generalizing it to other data types, that would be super.

So, an example of use:

```
julia> A = [1 2 3 4 5; 3 1 5 7 8; 2 -1 2 3 3]
3×5 Array{Int64,2}:
1 2 3 4 5
3 1 5 7 8
2 -1 2 3 3
julia> A = Rational.(A)
julia> Q,R,p = rrqr(A)
```

Comment: `rrqr`

finds a triplet `Q,R,p`

where `A*P = QR`

where `p`

is the column permutation vector, related to `P`

by P = I_{:,p}. Observe that to stay within the bounds of rational numbers, `Q`

is *orthogonal* but **not*** orthonormal.

So:

```
julia> Q
3×3 Array{Rational{Int64},2}:
5//1 121//98 1//3
8//1 -11//49 -1//3
3//1 -143//98 1//3
julia> R
3×5 Array{Real,2}:
1//1 15//98 5//14 85//98 61//98
0//1 1//1 -7//11 -3//11 -1//11
0//1 0//1 0//1 0//1 0//1
```

Observe that the vector of squared norm of columns of `Q`

is:

```
julia> colnorm2(Q)
3-element Array{Rational{Int64},1}:
98//1
363//98
1//3
```

If I define a diagonal scaling matrix `S`

with `sqrt.(colnorm2(Q))`

on the diagonal so that `A*P=Q/S*S*R`

, then `Q/S`

is *orthonormal*, while `S*R`

will have values along the main diagonal with decreasing size until the diagonal becomes zero. This is what is meant by rank revealing QR, as far as I know. Good implementations of RRQR are faster than SVD, and gives more or less similar information.

From matrix `R`

, we see that the main diagonal of `R`

contains 2 nonzero elements. This means that \mathrm{rank} A = 2. A set of basis vectors for the *column space* of `A`

is then given by the 2 first columns of Q:

```
julia> rA = Int(sum(diag(R)))
julia> C_A = Q[:,1:rA]
3×2 Array{Rational{Int64},2}:
5//1 121//98
8//1 -11//49
3//1 -143//98
```

The null space of `R`

can be found from `x`

satisfying `R*x=0`

. Some simple manipulation leads to:

```
julia> N_R = [-R[1:rA,1:rA]\R[1:rA,rA+1:end];Matrix{Rational}(I,size(A,2)-rA,size(A,2)-rA)]
5×3 Array{Rational,2}:
-5//11 -10//11 -7//11
7//11 3//11 1//11
true//true 0//1 0//1
0//1 true//true 0//1
0//1 0//1 true//true
```

This is the null space of `R`

. To find the null space of `A`

, it is necessary with a technical swapping of rows [I wrote a simple utility function `permswap()`], leading to:

```
julia> N_A = N_R[permswap(p),:]
5×3 Array{Rational,2}:
true//true 0//1 0//1
7//11 3//11 1//11
0//1 0//1 true//true
0//1 true//true 0//1
-5//11 -10//11 -7//11
```

A quick test of the null space:

```
julia> Nf = LinearAlgebra.nullspace(Float64.(A))
5×3 Array{Float64,2}:
-0.562308 0.0665452 -0.623681
-0.164975 0.0838044 -0.469905
0.201196 -0.880102 -0.125491
0.640077 0.445378 -0.225899
-0.454328 0.124928 0.568712
julia> rank([N_A Nf])
3
```