QR factorization of a low rank matrix

Suppose A is an matrix and suppose we know that the rank of A is bounded above by 2.

What is the most efficient method to calculate a QR factorisation of A with Q matrix and R matrix ?

[Bruno Luong]

You also need a permutation matrix P of (n x n) such that

A*P = Q*R

Otherwise the decomposition might not be possible.

Note a permutation matrix is a matrix P with entries 0/1 such that sum(P,1) == 1 and sum(P,2) == 1.

[Sara]

Yes, sorry, I also need the permutation matrix P

Answer 1

[Bruno Luong]

Just use MATLAB qr. I don't think you have a chance to beat Blas function with standard MATLAB program, even if you know rank <= 2

[Sara]

My hope was to somehow exploit the limitation on rank

Answer 2

[Bruno Luong]

Here is the Gram-Schmidt orthoganalisation algo that stops at the given ranks (arbitrary, not necessary 2).

Again I don't think it can beats MATLAB qr (I didn't test but I'm pretty sure that is the fact)

% Generat test low-rank matrix
m = 10;
n = 5;
rkmax = 2;
A = rand(m,rkmax)*rand(rkmax,n);

p = (1:n)';
Q = zeros(m,rkmax);
R = zeros(rkmax,n);
B = A;
for r = 1:rkmax
    [maxs,jmax] = max(sum(B.*conj(B),1)); % pivot position
    Rrr = sqrt(maxs);
    Qr = B(:,jmax)/Rrr;
    % swap columns to bring the pivot to column r
    j = [r, (r-1+jmax)];
    p(j) = p(j([2 1]));
    B(:,jmax) = B(:,1);
    B(:,1) = []; % As it move to jmax  we no longer need it
    Rr = Qr'*B;
    Q(:,r) = Qr;
    R(r,p(r:end)) = [Rrr, Rr];
    if r < rkmax % no need to update for the last iteration (where B should be close to O)
        B = B-Qr*Rr;
    end
end
P = accumarray([p, (1:n)'],1);
R = R(:,p);

% Check correctness
Q*R

ans = 10×5

0.6654 0.1944 0.4332 0.4734 0.2505 0.1606 0.4941 0.3577 0.1971 0.1451 0.7676 0.8249 0.8398 0.6573 0.4027 0.1540 0.1466 0.1578 0.1284 0.0772 0.0362 0.0306 0.0349 0.0295 0.0174 0.5014 0.8344 0.7159 0.4841 0.3190 0.6931 0.4340 0.5823 0.5360 0.3048 0.7426 0.4058 0.5904 0.5633 0.3153 0.4160 0.6712 0.5820 0.3977 0.2607 0.8260 0.2636 0.5504 0.5918 0.3152

Ap = A*P, % = A(:,p)

Ap = 10×5

0.6654 0.1944 0.4332 0.4734 0.2505 0.1606 0.4941 0.3577 0.1971 0.1451 0.7676 0.8249 0.8398 0.6573 0.4027 0.1540 0.1466 0.1578 0.1284 0.0772 0.0362 0.0306 0.0349 0.0295 0.0174 0.5014 0.8344 0.7159 0.4841 0.3190 0.6931 0.4340 0.5823 0.5360 0.3048 0.7426 0.4058 0.5904 0.5633 0.3153 0.4160 0.6712 0.5820 0.3977 0.2607 0.8260 0.2636 0.5504 0.5918 0.3152

norm(Q*R-Ap,'fro')/norm(A,'fro')

ans = 1.0286e-16

[Bruno Luong]

Some of the indexing in my code reduce the flops but might hurts the runtime (MATLAB is not fast when indexing). I did not bother to opimize the code, but if you want to play with it, you should gain some performance.

[Bruno Luong]

tic/toc with a big matrix to compare with MATLAB qr.

I build this case where MATLAB qr is very defavorable (matrix is size 1000 og rank 2)

m = 1000;
n = 1000;
rkmax = 2;
A = rand(m,rkmax)*(rand(rkmax,n)+1i*rand(rkmax,n));
tic
[Q,R,P] = qr(A);
toc

Elapsed time is 0.294079 seconds.

tic
p = (1:n)';
Q = zeros(m,rkmax);
R = zeros(rkmax,n);
B = A;
for r = 1:rkmax
    [maxs,jmax] = max(sum(B.*conj(B),1)); % pivot position
    Rrr = sqrt(maxs);
    Qr = B(:,jmax)/Rrr;
    % swap columns to bring the pivot to column r
    j = [r, (r-1+jmax)];
    p(j) = p(j([2 1]));
    B(:,jmax) = B(:,1);
    B(:,1) = []; % As it move to jmax  we no longer need it
    Rr = Qr'*B;
    Q(:,r) = Qr;
    R(r,p(r:end)) = [Rrr, Rr];
    if r < rkmax % no need to update for the last iteration (where B should be close to O)
        B = B-Qr*Rr;
    end
end
R = R(:,p);
toc

Elapsed time is 0.036170 seconds.

P = accumarray([p, (1:n)'],1);
Ap = A*P; % A(:,p)
norm(Q*R-Ap,'fro')/norm(A,'fro')

ans = 2.7891e-16

Answer 3

[Christine Tobler]

Let's start with computing a rank-2 approximation of a large matrix. In general, this would be done using the SVD, but if you know your matrix to be rank-2 or lower, a cheaper algorithm can work:

rng default;
A = randn(1e3, 2) * randn(2, 1e3); % Random rank-2 matrix

% Find column of maximum element by absolute value:
[~, linearIndex] = max(abs(A(:)));
[~, j] = ind2sub(size(A), linearIndex);
U = A(:, j);
U = U / norm(U);

V = U'*A;
Aremainder = A - U*V;

% Find column of maximum element by absolute value in remainder:
[~, linearIndex] = max(abs(Aremainder(:)));
[~, j2] = ind2sub(size(A), linearIndex);

U = orth([A(:, j) A(:, j2)]);
V = U'*A;

norm(A - U*V)

ans = 2.4235e-13

Now that we have this, let's make it fit as a Q and an R factor. U is already orthogonal, so we only need to make V be triangular:

[Qsmall, R] = qr(V);
Q = U*Qsmall;

norm(A - Q*R)

ans = 2.5439e-13

For sufficiently large matrices, this should be faster than calling QR. Otherwise, calling permuted QR and cutting off all but the first two columns of Q will be faster.