# QR factorization of a low rank matrix

11 ビュー (過去 30 日間)
Sara 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 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 2023 年 9 月 20 日
Yes, sorry, I also need the permutation matrix P

サインインしてコメントする。

### 回答 (3 件)

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
##### 1 件のコメント-1 件の古いコメントを表示-1 件の古いコメントを非表示
Sara 2023 年 9 月 20 日
My hope was to somehow exploit the limitation on rank

サインインしてコメントする。

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
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 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 日

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 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)
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 ExchangeLinear Algebra についてさらに検索

### Community Treasure Hunt

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

Start Hunting!

Translated by