Complete QR-decomposition

[atouminet]

Consider the following code

let a = arr2(&[[1.], [0.]]);
let (q, r) = a.qr().unwrap();
println!("q = {}", q);
println!("r = {}", r);

The output is

q = [[1],
 [0]]
r = [[1]]

however I was expecting q to be of size 2 x 2 and r of size 2 x 1. For instance, in Matlab, [q, r] = qr([1; 0]) prints

q =

   1   0
  -0   1

r =

   1
   0

[atouminet]

I am using the version 0.10.0 with the intel-mkl backend.

[termoshtt]

NumPy returns

In [5]: a
Out[5]: 
array([[1],
       [0]])

In [6]: np.linalg.qr(a)
Out[6]: 
(array([[1.],
        [0.]]), array([[1.]]))

Do you have reason that the Matlab answer is "correct"?

[atouminet]

By definition, the QR factorization of a matrix M of size m x n is M = QR where

• Q is an orthonormal matrix of size m x m
• R is an upper triangular matrix of size m x n

See https://en.wikipedia.org/wiki/QR_decomposition By definition an orthogonal matrix is always square. If you look at the numpy documentation (https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.linalg.qr.html) you can read that their qr function actually defaults to a reduced version of the real QR factorization.

After double checking the results what ndarray-linalg returns is actually that reduced decomposition. However, the reduced QR factorization contains less information than the full one, and should not be the default behavior of the function.

It is not possible with the current API to get the real QR factorization of a non-square matrix. I think it should be fixed by either:

• Modifying the current qr function to return the full factorization and modify the API to optionally compute the reduced factorization instead
• Keeping the current API and adding a new function to compute the full factorization

I think that defaulting to the full factorization is a better API design choice (as Matlab does).

[termoshtt]

I agree that we should have both option.

the reduced QR factorization contains less information than the full one,

This needs some discussion. QR factorization of a MxN matrix A=QR (M > N) has N vectors in R^M which cannot span full space. We need (M-N) more vectors to determine full-rank Q matrix, i.e. we have to orthogonalize the tangent space of Im(A), and this choice is additional information than A originally has, while reduced QR has same "amount" of information as A.

Keeping the current API and adding a new function to compute the full factorization

I prefer this option. The non-square matrices usually appear with large M (or N), and we should not allocate MxM matrix in these cases. Full-QR is natural for small matrices, but ndarray structure (i.e. unsized vector + strides) is suitable for large array basically.

[jedbrown]

Julia uses an elegant solution that retains the elementary reflector representation. This avoids the caller needing to ask up-front for complete or reduced representation. The resulting Q matrix would implement its operations using dormqr. It's cheaper, takes less storage, and sometimes more numerically stable than reducing to naive storage using dorgqr. The difference is extreme when you want a complete QR of a tall-skinny matrix.

[vlad17]

I was considering using ndarray-linalg for a least squares problem, and @jedbrown raises a good point! Householder info would be very convenient if that's an acceptable sidestep.

For least squares, the full QR representation may not even be representable. However, storing a sequence of householder transforms h_1,...,h_k is. Assuming pivoting is captured by a permutation pi, the least squares solution to ||Ax-y|| is then just pi(R^{-1} h_k(...(h_1(y))...) ) (where h_k...h_1 = Q and A pi = Q R)