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

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.

Sign in to comment.

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.

Jason
on 5 Jun 2013

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?

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.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.