QR factorization of a low rank matrix

6 visualizzazioni (ultimi 30 giorni)
Sara
Sara il 20 Set 2023
Modificato: Bruno Luong il 20 Set 2023
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 Commenti
Bruno Luong
Bruno Luong il 20 Set 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.
Sara
Sara il 20 Set 2023
Yes, sorry, I also need the permutation matrix P

Accedi per commentare.

Risposte (3)

Bruno Luong
Bruno Luong il 20 Set 2023
Modificato: Bruno Luong il 20 Set 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
  1 Commento
Sara
Sara il 20 Set 2023
My hope was to somehow exploit the limitation on rank

Accedi per commentare.


Bruno Luong
Bruno Luong il 20 Set 2023
Modificato: Bruno Luong il 20 Set 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
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 Commenti
Bruno Luong
Bruno Luong il 20 Set 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
Bruno Luong il 20 Set 2023
Modificato: Bruno Luong il 20 Set 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
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

Accedi per commentare.


Christine Tobler
Christine Tobler il 20 Set 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)
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.

Categorie

Scopri di più su Linear Algebra in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by