QR factorization of a low rank matrix

8 次查看（过去 30 天）
Sara
Sara 2023-9-20
编辑： Bruno Luong 2023-9-20
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 ?
  2 个评论
Bruno Luong
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.
Sara
Sara 2023-9-20
Yes, sorry, I also need the permutation matrix P

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Bruno Luong
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
Bruno Luong 2023-9-20
编辑：Bruno Luong 2023-9-20
Ran in:
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
  2 个评论
Bruno Luong
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
Bruno Luong 2023-9-20
编辑：Bruno Luong 2023-9-20
Ran in:
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

请先登录，再进行评论。

Christine Tobler
Christine Tobler 2023-9-20
Ran in:
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.

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by