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??]
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
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?
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.
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).
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.
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
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.
Hi, thanks for your interest.
I like the idea of rrqr – there is a code for MATLAB, but I haven’t looked at that. I’m not 100% sure that my code actually gives rrqr, although I suspect so (I would need to prove it). Properly implemented (and for Floats), rrqr is more efficient than svd, and gives similar information related to condition number (not the same information, but similar). Also, qr algorithms should work for sparse matrices. The Q matrix will probably not be sparse, so it would be useful to only save the rank A first columns of Q + the R matrix.
I have had to look at some other things, the last couple of weeks. I will get back to my code, but it needs some fixing…
If you have good command of Julia programming, I’m sure you can improve on my code. Give me some time into next week, and I’ll see if I can clean it up somewhat.
I would be thrilled to also see a Smith normal form (for rectangular matrix, not restricted to square matrices!). I did some Maple programming some 27 years ago using the Smith form – I looked into Linear Systems Theory (a la Kailath’s book) and the Smith-McMillan form. I don’t know how large systems can be studied via the Smith form, but it gives very useful info for linear dynamic systems. Incidentally, the reason why I looked into rrqr was to do some computations on coprime polynomial matrices.