[Octave-bug-tracker] [bug #51950] sparse qr may return a malformed sparse matrix R with out of bounds entries

[Joseph Young]

Follow-up Comment #14, bug #51950 (project octave):

>From what I can tell, Octave doesn't call cs_qr_mex.  However, it appears to
parrot how that file is constructed.  For example, the comments in
sparse-qr.cc about whether or not the double transpose to sort is necessary. 
This double transpose can also be found in cs_qr_mex.c.  That said, I don't
know the history of the files, so I don't know for sure.

Anyway, it looks like a number of things are going on that are simultaneously
breaking.  At its core, though, it's not clear to me that cs_qr can correctly
compute what Octave currently asks for it.  As I mentioned before, p78 of Tim
Davis' book Direct Methods for Sparse Linear Systems notes that extra rows can
be added when the matrix is structurally rank-deficient.  For how he calls the
routine, this isn't too much of a problem since he generally doesn't form Q. 
Rather, he keeps around the Householder reflectors and then uses them.  In
that way, it doesn't really matter if the extra rows are there because if we
solve when the factorization we can always remove the extra zeros from the
solution later.  In our case, this doesn't work since we want Q.

The easiest way to see why is to use cs_qr_mex.c from CXSparse along with
cs_qright.  Also, comment out the row/column check in cs_qr_mex.c to make this
simpler:

      //if (m < n) mexErrMsgTxt ("A must have # rows >= #
columns") ;

Anyway, we can ask CXSparse for the QR factorization with the following code:

% Create the matrix
A=sparse(5,6);A(3,1)=0.8;A(2,2)=1.4;A(1,6)=-0.5;

% Find the factorization
[V beta p R q]= cs_qr(A);

% Find Q
Q = cs_qright(V,beta,p,speye(size(V,1)));

Here, we get:

octave:6> full(Q)
ans =

   0.00000   0.00000   0.00000   0.00000   0.00000  -1.00000   0.00000  
0.00000
   0.00000  -1.00000   0.00000   0.00000   0.00000   0.00000   0.00000  
0.00000
  -1.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000  
0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   1.00000  
0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000  
1.00000
   0.00000   0.00000   1.00000   0.00000   0.00000   0.00000   0.00000  
0.00000
   0.00000   0.00000   0.00000   1.00000   0.00000   0.00000   0.00000  
0.00000
   0.00000   0.00000   0.00000   0.00000   1.00000   0.00000   0.00000  
0.00000

octave:7> full(R)
ans =

  -0.80000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000  -1.40000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.50000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000

Or really:

R =

Compressed Column Sparse (rows = 8, cols = 6, nnz = 3 [6.2%])

  (1, 1) -> -0.80000
  (2, 2) -> -1.4000
  (6, 6) ->  0.50000

This is where we see the element in the (6,6) position.  Notice, that these
matrices are much bigger than what we want.  We want Q to be 5x5 and R to be
5x6, but we added three rows.  Now, I believe the current QR factorization
code in Octave just truncates these matrices.  Note:
octave:107> [qqq rrr]=qr(A)
qqq =

   0.00000   0.00000   0.00000  -0.50000   0.00000
   0.00000  -1.00000   0.00000  -0.50000   0.00000
  -1.00000   0.00000   0.00000  -0.50000   0.00000
   0.00000   0.00000   0.00000  -0.50000   0.00000
   0.00000   0.00000   0.00000  -0.50000   0.00000

rrr =

Compressed Column Sparse (rows = 5, cols = 6, nnz = 3 [10%])

  (1, 1) -> -0.80000
  (2, 2) -> -1.4000
  (6, 6) ->  0.50000

This is bad for two reasons.  One, the (6,6) element can't exist in a 5x6
matrix.  Two, we have an all zero column in Q, so Q is not unitary (Q'*Q=I).

Now, can this be fixed?  I'm really not sure.  We can't safely truncate R
since adding additional rows may very well mean that nonzero elements are
pushed into the extra rows in the factorization.  Without changing the
algorithm, I'm not sure this can be prevented.  Again, if we don't need Q,
this is not necessarily a problem since we can still solve linear systems
using the Householder reflectors directly and then truncate the result. 
However, there are situations where we want Q.

One possible, bad, fix is to play a game with a Dulmage-Mendelsohn permutation
of R.  Basically, the DM permutation permutes things to try and obtain a block
diagonal structure.  We know that we should have this in R, so it should be
safe to use.  For example, the following code kind of fixes cs_qr_mex:

% Computes the QR factorization in a slightly odd way that sort of handles
% structurally rank deficient matrices.
function [Q R]=myqr(A)

    % Determine the size of the matrix
    [m n]=size(A);

    % Find the QR factorization
    [V beta p R q]= cs_qr(A);

    % Find Q
    Q = cs_qright(V,beta,p,speye(size(V,1)));

    % Find the Dulmage-Mendelsohn permutation of R
    [p q]=dmperm(R);

    % Reorder the rows of R and the columns of Q.  The goal here is to move
the
    % elements of R up in the case where it's structurually rank deficient
    % while keeping the block diagonal structure.
    Q=Q(:,p);
    R=R(p,:);

    % Truncate the size of Q and R
    Q = Q(1:m,1:m);
    R = R(1:m,1:n);
end

For example, if we try the code:

% Create the matrix
A=sparse(5,6);A(3,1)=0.8;A(2,2)=1.4;A(1,6)=-0.5;

% Find the factorization
[Q R]=myqr(A);

We get:

octave:2> Q
Q =

Compressed Column Sparse (rows = 5, cols = 5, nnz = 3 [12%])

  (3, 1) -> -1
  (2, 2) -> -1.0000
  (1, 3) -> -1

octave:3> R
R =

Compressed Column Sparse (rows = 5, cols = 6, nnz = 3 [10%])

  (1, 1) -> -0.80000
  (2, 2) -> -1.4000
  (3, 6) ->  0.50000

which is almost what we want.  Basically, the R matrix now looks good, but Q
is messed up.  Note, Q is not unitary and has zero columns.  However, here,
A=QR.  Anyway, it works because it essentially replaces the old factorization
A=QR with A=QP'PR where P is a permutation designed to move everything above
the added rows.  I think this would be better than what Octave currently does,
but it's not a great solution.

Likely a better solution would be to migrate the code from CXSparse to SPQR,
which uses a different algorithm where this isn't an issue.  I believe MATLAB
already uses SPQR; it's faster; and Octave is already using codes from
SuiteSparse.

Anyway, if someone has an easier solution, I'd like to hear.  Thanks for
taking a look at things.

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?51950>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/