QR factorization of a low rank matrix
4 次查看(过去 30 天)
显示 更早的评论
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 个评论
Bruno Luong
2023-9-20
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.
回答(3 个)
Bruno Luong
2023-9-20
编辑:Bruno Luong
2023-9-20
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
2023-9-20
编辑:Bruno Luong
2023-9-20
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 个评论
Bruno Luong
2023-9-20
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
2023-9-20
编辑:Bruno Luong
2023-9-20
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
2023-9-20
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 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!