# QR-like factorization preserving type?

Is there a factorization in the LinearAlgebra system which works similar to QR factorization in factoring A into A= T*R where T is a square, full rank matrix while R is an upper triangular matrix?

In particular, R should have types rational if the types of A are integers or rational? (If floats, one can just as well use QR factorization.) [I’d guess this is possible using the Gram Schmidt procedure without normalizing the Q matrix??]

1 Like

I don’t think the LinearAlgebra package has such a function; perhaps the following will work:

julia> function gs(H::Matrix{Td}) where {Td}
Q = copy(H)
n = size(H,2)
R = Matrix{Td}(I,n,n)
for i=1:n
for j=1:i-1
R[j,i] = H[:,i]'*Q[:,j]
Q[:,i]-= R[j,i]*Q[:,j]
end
end
return Q,R
end
gs (generic function with 1 method)

julia> N=2; a=rand(0:10,N,N);   Q,R = gs(a);   Q*R-a
2×2 Array{Int64,2}:
0  0
0  0

julia> N=2; a=rand(0:10,N,N).//rand(1:10,N,N); Q,R = gs(a); Q*R-a
2×2 Array{Rational{Int64},2}:
0//1  0//1
0//1  0//1

2 Likes

Thanks for inspiration. Your algorithm doesn’t work as I need it, though. In general, the matrix (a) will be rectangular, not square. With, e.g., a 2\times 3 matrix, your algorithm gives the wrong dimensions, and it is not clear what is Q and what is R. Consider he following example:

julia> A = [1 0 2;2 1 3]
2×3 Array{Int64,2}:
1  0  2
2  1  3

julia> Q,R = gs(A);
julia> Q
2×3 Array{Int64,2}:
1  -2  -32
2  -3  -52

julia> R
3×3 Array{Int64,2}:
1  2    8
0  1  -13
0  0    1

julia> Q*R
2×3 Array{Int64,2}:
1  0  2
2  1  3


The correct matrices when using Julia’s LinearAlgebra routine is, however:

julia> F = qr(A);
julia> qr.Q
2×2 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
-0.447214  -0.894427
-0.894427   0.447214

julia> qr.R
2×3 Array{Float64,2}:
-2.23607  -0.894427  -3.57771
0.0       0.447214  -0.447214


As seen, the gs() Q matrix has the wrong dimension, and the gs() R matrix has the wrong structure (R should be upper triangular, Q should not).

Inspired by your algorithm, I found a Wikipedia page about the QR matrix, which included MATLAB code. Comparing these, and looking into the Gram Schmidt procedure, I observe that the line:

R[j,i] = H[:,i]'*Q[:,j]


needs be changed to:

R[j,i] = H[:,i]'*Q[:,j]/(Q[:,j]'*Q[:,j])


to ensure that the columns of Q are orthogonal. Then, I need to pick out Q[:,1:2]Q, and compute R as Q\A. The R from the gs() algorithm is not really needed, I think. Thus:

julia> A1 = Rational.(A)
julia> Q1,R1 = gs(A1)
julia> Q1
2×3 Array{Rational{Int64},2}:
1//1  -2//1  -32//1
2//1  -3//1  -52//1

julia> R_new = Q1[:,1:2]\A1
2×3 Array{Rational{Int64},2}:
1//1  2//1   0//1
0//1  1//1  -1//1

julia> Q1[:,1:2]*R_new
2×3 Array{Rational{Int64},2}:
1//1  0//1  2//1
2//1  1//1  3//1


However, for my use, the A matrix will, in general, have both linearly dependent rows and columns. So it is necessary to use column pivoting. I’ll think it over.

Anyway, inspiring attempt!

Would Givens rotations but without the normalisation work? That is, multiply by [a b; -b a] to introduce zeros?

2 Likes

Possibly. I need to check whether it works while keeping the field type of the system (or rather: preserve rational types, as a minimum). Thanks for suggestion.

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


I think an improved version of function rrqr() would be useful.

1. Rank revealing QR is superior to standard QR since it gives information about the conditioning of the system. And with a good implementation, it should be faster than SVD.
2. Educationally, I think RRQR is useful in that it serves to clarify concepts such as column space (range) and null space (kernel).
3. By relaxing on the requirement that Q should be orthonormal, it is possible to do exact computations. This is useful in quite a few applications.
4. Most likely, using Rational numbers requires BigInt for larger systems than academic examples. But the idea of rank revealing QR is worthwhile also when using Floats.

If anyone is interested in my code, let me know. It is not huge in any way – a few tens of lines, so I can include it here somewhere. However, it needs improvement in many areas (efficiency, memory allocation, handling types better, better names for variables and functions, etc.). I used the Gram-Schmidt procedure; it would be interesting to check if something similar/better can be achieved with Givens rotation, etc. I’m not good at code optimization – I even struggle with finding errors (huge thanks to the people behind VS code, Julia extension with debugger – it is not perfect, but to me, it represents a milestone).

1 Like

Hey Bernt, I would be interested in your code, it seems to be exactly what I’m looking for! In particular, I’m attempting to decompose a large (but rank-2) matrix into a sum of two rank-1 matrices (outer products of vectors) where the vectors are also rational. I think that can be done with the rrqr you’re describing.

Edit PS: Have you noticed this? Rational matrix division in Julia - Stack Overflow
Perhaps an implementation of this already exists in Julia so that rational-preserving linear solves can be done.

1 Like

Why can’t you just use SVD?

An SVD uses an iterative method, and I need a direct method to ensure rational entries in my vectors

Edit:
I threw together this quickly, which does column-pivoted "Q"R (where Q has orthogonal but not normalized columns) on square and short+fat matrices. Apologies for it being slightly trashy, I’m still a Julia newbie It seems it’d be easier to integrate if there was any kind of name for this factorization.

function proj(a, u)
return u*(dot(a,u)//dot(u,u))
end

function rationalrrqr(A)
mat = copy(A)
Q = Array{Rational{BigInt}}(undef, size(mat, 1), size(mat, 1))
perm = Matrix(1I, size(mat,2), size(mat, 2))

for col in 1:size(mat, 1)
maxnorm = 0
maxcol = col
for piv in col:size(mat,2)
pivcoln = dot(mat[:, piv], mat[:, piv])
if pivcoln > maxnorm
maxnorm = pivcoln
maxcol = piv
end
end

swap = copy(mat[:, col])
mat[:,col] = copy(mat[:,maxcol])
mat[:,maxcol] = swap
swap = copy(perm[:, col])
perm[:,col] = copy(perm[:,maxcol])
perm[:,maxcol] = swap

u = mat[:, col]
u_copy = copy(u)
for j in 1:(col-1)
u -= proj(u_copy, Q[:, j])
end
if norm(u) != 0
Q[:, col] = u
else
u = rand(1:10000, size(mat, 1), 1)
u_copy = copy(u)
for j in 1:(col-1)
u -= proj(u_copy, Q[:, j])
end
Q[:, col] = u
end
end
R = Q \ mat
return Q, R, perm
end


2 Likes

There is something known as the Smith decomposition or Smith normal form (not to be confused with Schmidt decomposition, which is essentially the singular value decomposition). The Smith normal form is defined for numbers which do not form a field, it is sufficient if they have the structure of a principal ideal domain. Integers have this structure. However, what you doing is not that, so maybe the Smith normal form is too much for what you need for your application.

2 Likes