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? https://stackoverflow.com/questions/20985783/rational-matrix-division-in-julia
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 :stuck_out_tongue:

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

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…

  • The implementation is probably not very efficient. I’m not a super programmer, and I’m sure a good programmer could make the code faster. I’m aware of a few tricks that would make it better.
  • Currently, my code assumes that the matrix elements are rational numbers. It would be good to allow for any types, and then return rational numbers (or integers with a scaling diagonal) when the input type is rational (integer), or floating point numbers when the matrix elements are floats. Should probably be done by type dispatching.
  • I have used the Gram-Schmidt procedure. I’m pretty sure that for matrices larger than toy examples, it would be necessary to use BigInt. Perhaps using Givens rotation without scaling tackles this better? (Someone suggested that.) However, when I compute with rational numbers, there is, of course, no rounding error.

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.