QR factorization of a low rank matrix
4 views (last 30 days)
Show older comments
Suppose A is an
matrix and suppose we know that the rank of A is bounded above by 2.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1487917/image.png)
What is the most efficient method to calculate a QR factorisation of A with Q matrix
and R matrix
?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1487922/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1487927/image.png)
2 Comments
Bruno Luong
on 20 Sep 2023
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.
Answers (3)
Bruno Luong
on 20 Sep 2023
Edited: Bruno Luong
on 20 Sep 2023
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
Bruno Luong
on 20 Sep 2023
Edited: Bruno Luong
on 20 Sep 2023
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
Ap = A*P, % = A(:,p)
norm(Q*R-Ap,'fro')/norm(A,'fro')
2 Comments
Bruno Luong
on 20 Sep 2023
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
on 20 Sep 2023
Edited: Bruno Luong
on 20 Sep 2023
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
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
P = accumarray([p, (1:n)'],1);
Ap = A*P; % A(:,p)
norm(Q*R-Ap,'fro')/norm(A,'fro')
Christine Tobler
on 20 Sep 2023
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)
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)
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.
0 Comments
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!