Sparse and nonspare QR

6 次查看（过去 30 天）
Jason
Jason 2013-6-5
I have question regarding the difference between the sparse and non-sparse QR as implimented R2012b. I have a matrix, X, that is mostly zeros which I wish to upper-triangularize. Because the matrix X is formed piecemeal from several sparse multiplications, the X is sparse. The results of:
a1 = qr(X)
and
a2 = triu(qr(full(X)))
are completely different. The few nonzero values end up in very different places. The code is for a square root DWY backward (Fisher type) filter as in Park and Kailath 96. Based on what I know of the problem and what the answer should be, a2 is correct (or at least I know a1 is incorrect because what should be the multiplication of 2 nonzero submatrices is zero in a1 but not a2).
Does anyone know why the two are different? Further is there a way to force qr(X) to produce the same result as triu(qr(full(X))) without resorting to the full() call?
Thank you

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

采纳的回答

Richard Brown
Richard Brown 2013-6-14
Sorry to take a while, busy week. Neither are wrong. You'll notice that the only difference in the outputs is that the sparse version puts all the rows with nonzero diagonals at the top of the matrix, while the full version leaves them in place. Thus, to get the nonzero pattern you want from the sparse matrix, do something like this
R = qr(A); % A is sparse
[I, J, S] = find(R); % Dismantle R
pivots = accumarray(I, J, [], @min); % Find min col in each nz row
I = pivots(I); % Translate the rows
R = sparse(I, J, S, size(Rs, 1), size(Rs, 2));
The matrices should be now be the same up the the signs of the rows. You can probably do this slightly more efficiently, but it should do the trick, avoiding full matrices.

更多回答（2 个）

James Tursa
James Tursa 2013-6-5
编辑：James Tursa 2013-6-5
What happens when you compare the complete result, e.g.,
[a1 r1] = qr(X);
[a2 r2] = qr(full(X));
Then compare a1*r1 with a2*r2. There may be some sign ambiguity in the factored results that you cannot control.
  1 个评论
Jason
Jason 2013-6-5
Both qr's are "correct" in that sense (a1*r1 == a2*r2). a2 is an identity (except for the last 3 columns which are 0). r2 matches the "q-less" qr (from dgeqrf.f), r1 doesn't. In a perfect world, I would like the r's from the q-less sparse and non-sparse qr to match. I need to take submatrices out of r. If the r is permuted this is obviously a problem. I guess this is not possible?

请先登录，再进行评论。

Richard Brown
Richard Brown 2013-6-5
编辑：Richard Brown 2013-6-5
On my system they match (up to the signs of the rows, as you'd expect).
m = 20;
n = 10;
X = sprand(m,n,0.3);
Rs = qr(X);
Rs = spdiags(sign(spdiags(Rs, 0)),0,n,n) * Rs(1:n,1:n);
Rf = triu(qr(full(X)));
Rf = diag(sign(diag(Rf))) * Rf(1:n,1:n);
disp(norm(Rf - Rs));
Can you possibly provide a MWE? Is your matrix full rank?
  2 个评论
Jason
Jason 2013-6-6
On my matrix the above results in norm(rf-rs) = 1.31999599941703
No, the matrix is not full rank.
Below is a similar problem. It is the start of an information filter so there are more zeros then there will be later. After the qr, rf(58:114,58:114) and rf(58:114,115) should have nonzero values. Once this happens, block4, block5, and block6 will be non-zero.
tmp1 = cov(randn(500,3));
tmp2 = cov(randn(500,3));
block1 = blkdiag(chol(inv(tmp1)),zeros(54,54));
block2 = zeros(57,57);
block3 = blkdiag(chol(tmp1)',zeros(54,54))*randn(57,1);
block4 = zeros(57,57);
block5 = zeros(57,57);
block6 = zeros(57,1);
block7 = zeros(3,57);
block8 = chol(inv(tmp2))*randn(3,57);
block9 = chol(inv(tmp2))*randn(3,1);
a = [block1,block2,block3;block4,block5,block6;block7,block8,block9];
b = sparse(a);
rs = qr(b);
rf = triu(qr(a));
figure(1);
imagesc(rf);
figure(2);
imagesc(full(rs));
Richard Brown
Richard Brown 2013-6-7
thanks - I didn't get time to get back to this question this week, but try back on Monday sometime!

请先登录，再进行评论。

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by