octave-bug-tracker
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Octave-bug-tracker] [bug #51950] sparse qr may return a malformed spars


From: Joseph Young
Subject: [Octave-bug-tracker] [bug #51950] sparse qr may return a malformed sparse matrix R with out of bounds entries
Date: Thu, 7 Sep 2017 02:42:59 -0400 (EDT)
User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:55.0) Gecko/20100101 Firefox/55.0

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

This is likely an issue inherited from the original cs_qr_mex file.  This call
is a little wonky:

octave:2> A=sparse(5,6);A(3,1)=0.8;A(2,2)=1.4;A(1,6)=-0.5;
octave:3> [V beta p R q]= cs_qr(A)
error: cs_qr: A must have # rows >= # columns

If we add the additional row to force it to work:

octave:4> A=[A;sparse(1,6)]
A =

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

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

octave:5> [V beta p R]= cs_qr(A)
V =

Compressed Column Sparse (rows = 9, cols = 6, nnz = 6 [11%])

  (1, 1) ->  1.6000
  (2, 2) ->  2.8000
  (3, 3) ->  1
  (4, 4) ->  1
  (5, 5) ->  1
  (6, 6) -> -1

beta =

   0.78125
   0.25510
   0.00000
   0.00000
   0.00000
   2.00000

p =

   3   2   7   8   9   1   4   5   6

R =

Compressed Column Sparse (rows = 9, cols = 6, nnz = 3 [5.6%])

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

Notice that we have the magical (6,6) element.  However:

octave:6> [q r]=qr(full(A))
q =

   0   0   1   0   0   0
   0   1   0   0   0   0
   1  -0   0   0   0   0
   0  -0  -0   1   0   0
   0  -0  -0  -0   1   0
   0  -0  -0  -0  -0   1

r =

   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.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
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000

Candidly, it shouldn't be there.  What's happening?  Other than cs_qr_mex.c
likely having a bug, there is documentation in cs_qr.m:

If A is structurally rank deficient, additional empty rows may have been added
to V and R.

In other words, CXSparse can and will add additional rows to R when it deems
fit.  This is discussed on page 78 in Tim Davis' book Direct Methods for
Sparse Linear Systems:

The cs_qr function uses the symbolic analysis computed by cs_sqr: the column
elimination tree S->parent, column preordering S->q, row permutation
S->pinv, the S->leftmost array, the number of nonzeros in R and V
(S->unz and S->lnz, respectively), and the number of rows S->m2 after
adding fictitious rows if A is structurally rank deficient.

Here's an example of the call in C:

#include "cs.h"

int main() {
    // Create a 5x6 matrix with a max of 10 values.  Allocate memory for the
    // values and denote that we're using triplet form. 
    cs * A = cs_spalloc(5,6,10,1,1);

    // Set the entries of A
    cs_entry(A,2,0,0.8);
    cs_entry(A,1,1,1.4);
    cs_entry(A,0,5,-0.5);
    #if 1
    cs_entry(A,0,2,0.0);
    cs_entry(A,0,3,0.0);
    cs_entry(A,0,4,0.0);
    #endif

    // Convert A to a compressed matrix
    cs * Ac = cs_compress(A);

    // Perform a symbolic QR factorization of A.  The 0 means natural
ordering
    // and the 1 means we're doing a sparse qr factorization.
    css * As = cs_sqr(0,Ac,1);

    // Perform a QR factorization of A
    csn * QR = cs_qr(Ac,As);

    // Free memory from A and Ac
    cs_spfree(A);
    cs_spfree(Ac);

    // Free the symbolic factorization
    cs_sfree(As);

    // Free the QR factorization
    cs_nfree(QR);
}

In gdb we can see:

(gdb) print QR[0].U[0]
$4 = {nzmax = 3, m = 8, n = 6, p = 0x100203550, i = 0x100203580, 
  x = 0x100203320, nz = -1}
(gdb) print QR[0].U[0].p[0] @ 7
$5 = {0, 1, 2, 2, 2, 2, 3}
(gdb) print QR[0].U[0].i[0] @ 3
$6 = {0, 1, 2}
(gdb) print QR[0].U[0].x[0] @ 3
$7 = {-0.80000000000000004, -1.3999999999999999, -0.5}

In other words, we have

R(1,1) = -0.8
R(2,2) = -1.4
R(6,3) = -0.5

which is what we want except for that R has size 8 x 6 since CXSparse added 3
fictitious rows since we were structurally rank deficient.

How do we fix this?  Unless I'm missing something, I can't seem to find a
routine in CXSparse that drops the fictitious rows automatically.  We can
probably just calculate the size by hand and ignore what comes out in the cs
matrix since the number of rows may be greater than what we want.

    _______________________________________________________

Reply to this item at:

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

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




reply via email to

[Prev in Thread] Current Thread [Next in Thread]