Sparse and nonspare QR

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 1

[Richard Brown]

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.

Answer 2

[James Tursa]

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]

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?

Answer 3

[Richard Brown]

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

thanks - I didn't get time to get back to this question this week, but try back on Monday sometime!