MATLAB Answers

Jason
0

Sparse and nonspare QR

Asked by Jason
on 5 Jun 2013
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

  0 Comments

Sign in to comment.

3 Answers

Answer by Richard Brown on 14 Jun 2013
 Accepted Answer

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.

  0 Comments

Sign in to comment.


James Tursa
Answer by James Tursa
on 5 Jun 2013
Edited by James Tursa
on 5 Jun 2013

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 Comment

Jason
on 5 Jun 2013
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?

Sign in to comment.


Answer by Richard Brown on 5 Jun 2013
Edited by Richard Brown on 5 Jun 2013

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 Comments

Jason
on 6 Jun 2013
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 on 7 Jun 2013
thanks - I didn't get time to get back to this question this week, but try back on Monday sometime!

Sign in to comment.