bugGNU Octave - Bugs: bug #56232, Octave crash when inverting an...

 
 

bug #56232: Octave crash when inverting an empty sparse matrix.

Submitted by:  Kai Torben Ohlhus <siko1056>
Submitted on:  Sun 28 Apr 2019 06:38:51 PM UTC  
 
Category:  Octave Function Severity:  4 - Important
Priority:  5 - Normal Item Group:  Segfault, Bus Error, etc.
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  Closed
Release:  5.1.0 Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Sat 18 May 2019 09:45:41 PM UTC, comment #62: 

Thanks to Marco for taking the lead in resolving this problem.  I made some modifications for syntax or for performance, and checked it in here (https://hg.savannah.gnu.org/hgweb/octave/rev/6e18f0ce268c).

Documenting my additions, I changed the BIST tests which were using some side effects of our implementation to work, to use more ordinary syntax.

Instead of

+%!test
+%!warning <matrix singular> A = inv (sparse ([1, 2;0 ,0]));
+%! assert (A, sparse ([Inf, Inf; 0, 0]));

I wrote

+%!test
+%! fail ("A = inv (sparse ([1, 2;0 ,0]))", "warning", "matrix singular");
+%! assert (A, sparse ([Inf, Inf; 0, 0]));

For the diagonal matrices, the implementation relied on creating a local buffer that could hold up to nnz() elements.  The algorithm also needed to cycle over the entire number of diagonal elements before arriving at a decision about whether the matrix was singular.  The initial implementation is below.

info = 0;
OCTAVE_LOCAL_BUFFER (octave_idx_type, nzelem, len);
octave_idx_type count_nz = 0;
for (octave_idx_type i = 0; i < len; i++)
  {
    if (elem (i, i) == 0.0)
      retval.elem (i, i) = octave::numeric_limits<double>::Inf ();
    else
      {
        retval.elem (i, i) = 1.0 / elem (i, i);
        nzelem[count_nz] = i;
        count_nz++;
      }
  }
if (count_nz == 0)
  {
    (*current_liboctave_error_handler)
      ("inverse of the null matrix not defined");
  }
else if (count_nz < len)
  {
    info = -1;
    for (octave_idx_type i = 0; i < count_nz; i++)
      retval.elem (nzelem[i], nzelem[i]) = octave::numeric_limits<double>::Inf ();
  }

Instead, I found that it was possible with two count variables to eliminate the OCTAVE_LOCAL_BUFFER and to short-circuit out of the loop as soon as a zero is detected.  For a very small speed-up, I also used xelem() to access elements where I could be certain that I was not going to go outside the bounds of the array.

info = 0;
octave_idx_type len = r;        // alias for readability
octave_idx_type z_count  = 0;   // zeros
octave_idx_type nz_count = 0;   // non-zeros
for (octave_idx_type i = 0; i < len; i++)
  {
    if (xelem (i, i) == 0.0)
      {
        z_count++;
        if (nz_count > 0)
          break;
      }
    else
      {
        nz_count++;
        if (z_count > 0)
          break;
        retval.elem (i, i) = 1.0 / xelem (i, i);
      }
  }
if (nz_count == 0)
  {
    (*current_liboctave_error_handler)
      ("inverse of the null matrix not defined");
  }
else if (z_count > 0)
  {
    info = -1;
    element_type *data = retval.fortran_vec ();
    std::fill (data, data + len, octave::numeric_limits<double>::Inf ());
  }

For Sparse and CompleSparse arrays, the original additional code was

+          if (rcond == 0.)
+            {
+              for (octave_idx_type i = 0; i < cols () + 1; i++)
+                ret.cidx (i) = cidx (i);
+
+              for (octave_idx_type i = 0; i < nnz (); i++)
+                {
+                  ret.data (i) = octave::numeric_limits<double>::Inf ();
+                  ret.ridx (i) = ridx (i);
+                }
+
+              return ret;
+            }

Instead of using hand-rolled for loops, I decided on functions from the C++ standard library.

+          if (rcond == 0.0)
+            {
+              // Return all Inf matrix with sparsity pattern of input.
+              octave_idx_type nz = nnz ();
+              ret = SparseComplexMatrix (rows (), cols (), nz);
+              std::fill (ret.xdata (), ret.xdata () + nz,
+                         octave::numeric_limits<double>::Inf ());
+              std::copy_n (ridx (), nz, ret.xridx ());
+              std::copy_n (cidx (), cols () + 1, ret.xcidx ());
+
+              return ret;
+            }

Rik <rik5>
Project Administrator
Fri 17 May 2019 04:03:14 PM UTC, comment #61: 

Thanks, I can confirm, that the fourth patch of Marco (file #46913) has no failing tests and all examples of this report result in a useful error or warning message.

Kai Torben Ohlhus <siko1056>
Project Member
Thu 16 May 2019 07:26:29 AM UTC, comment #60: 

@Rik, comment #58: I agree to keep the behavior of inv(zeros(2)). Then, diagonal matrices are the most prominent example of sparse matrices (although they are displayed in octave as full), so I would treat the singular case as for sparse matrices. I fixed the test, hence.

@Kai, comment #59: so mattype.ishermitian() really means "is the matrix hermitian AND nonsingular", good to know. Then, I removed the check of rcond in that branch.

I attach my fourth patch, thanks for your feedback.

(file #46913)

Marco Caliari <caliari>
Project Member
Wed 15 May 2019 08:59:33 PM UTC, comment #59: 

@Marco: Thanks for the patch.

Same observation with "make check" and opinion regarding Matlab compatibility as Rik in comment #58.

Regarding comment #57 1): Do you mean, when this branch is used?

https://hg.savannah.gnu.org/hgweb/octave/file/9326c2258e60/liboctave/array/CSparse.cc#l1005

By inserting a std::cout "YEAH sparse!!!", I could verify, that this branch is taken for the following matrix A:

N = 100;                                     # dimension
A = sprand(N, N, 4/N) + sprand(N, N, 4/N)*i; # complex randomness
A(1:(N+1):end) = 1:N;                        # set diagonal real > 0
A = sparse (triu (A) + triu (A, 1)');        # Hermitian

# Properties of A
matrix_type (A), ishermitian (A), iscomplex (A), issparse (A)
ans = Positive Definite
ans = 1
ans = 1
ans = 1

inv(A);

YEAH sparse!!!

Does this help?  I think you cannot enter this branch by using a singular matrix, when stepping through MatrixType

https://hg.savannah.gnu.org/hgweb/octave/file/9326c2258e60/liboctave/array/MatrixType.cc#l235

the necessary value "maybe_hermitian = true;" is only set in the "if (! singular)" branch.

Kai Torben Ohlhus <siko1056>
Project Member
Wed 15 May 2019 05:46:24 PM UTC, comment #58: 

@Marco: This is looking good.

On the question of inv(zeros(2)), I think we should retain current behavior and Matlab compatibility and return a matrix of all Inf.

I ran 'make check' and got this one failing test

>>>>> processing /home/rik/wip/Projects_Mine/octave-dbg/test/diag-perm.tst
***** test
x = diag (1:3);
assert (inv (x), diag ([1 1/2 1/3]));
x = diag (0:2);
assert (inv (x), diag ([Inf 1 1/2]));
!!!!! test failed
ASSERT errors for:  assert (inv (x),diag ([Inf, 1, 1 / 2]))

  Location  |  Observed  |  Expected  |  Reason
   (2,2)         Inf           1         'Inf' mismatch
   (3,3)         Inf          0.5        'Inf' mismatch

In Matlab and Octave,

x = full (diag ([Inf, 1, 1/2]))
inv (x)
x =

       Inf   0.00000   0.00000
   0.00000   1.00000   0.00000
   0.00000   0.00000   0.50000

inv (x)
warning: matrix singular to machine precision
ans =

   0  -0  -0
   0   1  -0
   0   0   2

 

Should the implementation of diagonal matrices be more like sparse (Inf whenever the element is non-zero) or should it be more like full matrices (the result above which is 1./diagonal_matrix)?

This is simple to fix

Rik <rik5>
Project Administrator
Wed 15 May 2019 07:40:56 AM UTC, comment #57: 

With this third patch, the warnings are triggered for the singular diagonal cases. It remains:

1) how to trigger the sparse hermitian case (to test it and to understand what mattype.ishermitian() is)
2) decide if we want that inv(zeros(2)) returns an error. I was initially in favor, but maybe in this case it is better to preserve the compatibility with matlab.

(file #46911)

Marco Caliari <caliari>
Project Member
Tue 14 May 2019 06:14:31 PM UTC, comment #56: 

Please consider sparse-singular2.cset, instead. NOW the diag case behaves like the sparse one, although I do not know how to give the singular matrix warning.

(file #46908)

Marco Caliari <caliari>
Project Member
Tue 14 May 2019 03:09:57 PM UTC, comment #55: 

I started the implementation (attached). It seems to work, but I cannot trigger the case mattype.ishermitian() in SparseMatrix::inverse. For instance, the following

> A = sparse ([1, 0, 0;0 ,0, 0; 0, 0, 1];
> ishermitian (A)
ans = 1

falls in ! mattype.ishermitian ().

I copied the result for the diag case and fixed a misalignment between the double and the complex version. Should I continue with inv (zeros (2)) returning an error? In Matlab and in Octave now the result is Inf (2).

(file #46906)

Marco Caliari <caliari>
Project Member
Mon 13 May 2019 04:40:26 PM UTC, comment #54: 

@Marco: No real comments from Maintainer's list so you are safe to implement the solution you suggested.

Rik <rik5>
Project Administrator
Fri 10 May 2019 11:33:10 AM UTC, comment #53: 

Thanks Marco for starting the thread.

https://lists.gnu.org/archive/html/octave-maintainers/2019-05/msg00053.html

Everything done, closing report.

Kai Torben Ohlhus <siko1056>
Project Member
Thu 09 May 2019 08:32:10 PM UTC, comment #52: 

I agree that the original bug report, a segfault, is resolved.  Could someone on this report summarize the issues on the Octave Maintainer's list so that we can get additional perspectives?

Rik <rik5>
Project Administrator
Thu 09 May 2019 07:46:52 PM UTC, comment #51: 

Thank you all for taking care.

@Rik: Your patches solve my original issue:

>> C = sparse (2, 2, 0); inv (C)
error: inv: division by zero
>> C = sparse (ones (2)); C = C - C; inv (C)
error: inv: division by zero


>> ver
----------------------------------------------------------------------
GNU Octave Version: 5.1.1 (hg id: c0d8ce61c1c9)
GNU Octave License: GNU General Public License
Operating System: Linux 4.12.14-lp150.12.58-default #1 SMP Mon Apr 115:20:46 UTC 2019 (58fcc15) x86_64
---------------------------------------------------------------------

Personally, I like the current error behavior of Octave.  I am surprised, that Matlab actually creates Matrices with a huge fill-up when inverting a sparse or empty matrix.  But as Marco says this might be a discussion for the developer list.

The segfault is fixed, so this item for me :-)

Kai Torben Ohlhus <siko1056>
Project Member
Thu 09 May 2019 07:32:44 PM UTC, comment #50: 

@Rik:

1) I see that in the full case, octave already returns a matrix full of Inf values. Preserving the sparsity pattern in the sparse case is the most similar thing I can imagine (although it does not work in the null case, I know). Anyway, no real high value.

2) If I had to check whether the inversion of a matrix was successful, as first option I would test the result with isinf, not isnan. But if we decide for NaN, we should change the behavior in the full case, too.

For things like these, shouldn't we use the developers mailing list?

Marco Caliari <caliari>
Project Member
Thu 09 May 2019 06:33:45 PM UTC, comment #49: 

I agree with Marco.  Although the inverse could, in general, be a fully dense matrix.  It is only worth paying the computation time, and the memory cost, if you are deriving information and getting a usable result.  If the result is meaningless then we might as well pick a return value that we like.

I tried some sample singular matrices, including the empty matrix, in Matlab and I don't see a pattern to the results.  The matrices themselves are fully populated with entries consisting only of -Inf, Inf, NaN.  So no real guidance there.

Two philosophical questions:

1) How high is the value of preserving the sparsity pattern in the output?

2) Do we want to poison downstream calculations, so that even if a user missed the warning about a singular matrix, the end result is clearly incorrect?  That would favor using NaN rather than Inf as a fill value.

Rik <rik5>
Project Administrator
Thu 09 May 2019 03:00:32 AM UTC, comment #48: 

I took care of reversing the hack that was put in to stop segfaults witch cholmod.  See changeset https://hg.savannah.gnu.org/hgweb/octave/rev/df9e3cb79379.

Final question to resolve is what to return when inverting a singular matrix and an empty matrix.

Rik <rik5>
Project Administrator
Thu 09 May 2019 01:29:11 AM UTC, comment #47: 

I think I took care of the MEX api with this changeset (https://hg.savannah.gnu.org/hgweb/octave/rev/23761e83756f).

Rik <rik5>
Project Administrator
Wed 08 May 2019 10:39:17 PM UTC, comment #46: 

Remaining things to review

1) Use of &dummy in dSparse.cc, CSparse.cc, sparse-chol.cc

I think these were put in when the data pointer for a sparse matrix could be NULL.  This can't happen anymore and so can probably be removed.

2) Use of nzmax in mex.cc

I didn't look at this, but I think mxSetNzmax (or somewhere in the chain) needs to be modified so that it can't set a value of 0.  Also, the constructors for class mxArray_sparse probably need to be modified like the Octave ones in Sparse.h so that they always create valid pointers.

Rik <rik5>
Project Administrator
Wed 08 May 2019 10:32:19 PM UTC, comment #45: 

It took several hours to verify all uses of nzmax in Octave.  I think changing Octave to universally reserve memory for a single value is okay.  I checked in the change on the stable branche here (https://hg.savannah.gnu.org/hgweb/octave/rev/c0d8ce61c1c9).

One nice thing is that it fixed a bad interaction with the amd libraries.  Previously,

octave:1> symamd ([])

symamd version 2.9, May 4, 2016: ERROR.  Array A (row indices of matrix) not present.
error: symamd: internal error!

But now,

octave:1> symamd ([])
ans = [](1x0)

which is the same as Matlab.

Rik <rik5>
Project Administrator
Wed 08 May 2019 06:37:09 PM UTC, comment #44: 

Sorry if I insist, but what about if we return a sparse matrix, with exactly the same number of element different from zeros of the input matrix, with those elements set to Inf?

invA = spalloc (rows (A), columns (A), nnz (A));
invA (A != 0) = Inf;

The only problem would be the sparse null matrix, again...

Marco Caliari <caliari>
Project Member
Wed 08 May 2019 06:00:00 PM UTC, comment #43: 

It does seem bad to return a full array, but the same problems can happen for non-singular sparse matrices since inversion can create a completely full matrix.

In any case, we do need to avoid the segfault.  It seems fine to me to always allocate the index and data pointers and ensure nzmax is at least 1 on both stable and default.

John W. Eaton <jwe>
Project Administrator
Wed 08 May 2019 05:17:46 PM UTC, comment #42: 

I think we can begin to proceed on resolving some parts of this bug.  I believe all of these changes should go on the dev branch.  An open question is whether we want to backport something minimal just to stop the segfault in 5.1.0.  That seems useful as the interpreter shouldn't segfault (losing user variables) for innocuous, valid inputs like C - C.

The first part of the fix on dev is a move to maintaining memory at all times for at least one element.  This will ensure valid data, ridx, and cidx pointers at all times.  Any routine which sets nzmax needs to ensure that a value of 0 is replaced with a value of 1.

Rik <rik5>
Project Administrator
Wed 08 May 2019 05:06:16 PM UTC, comment #41: 

In one sense, the current behavior is actually nicer.  When you attempt to invert a sparse singular matrix it prints an error message and then refuses to waste memory creating a full array filled with not particularly useful values (Inf or NaN).

I would be comfortable with this behavior (but, make the message explicit that this was a singular matrix rather than division-by-zero).  There are many other instances where Octave does more extensive input validation and doesn't perform a garbage in / garbage out calculation.

But, if we really want to return something, then returning an effectively full array seems like a bad idea.  There will be fat-fingered people who create an empty large spare array and then mistakenly try to invert it.  If the size is large enough the request for memory will fail immediately.  That's the best outcome.  If the size is larger than main memory but less than main + swap then Octave will proceed and I find that it can be a minute or more before a Ctrl+C is recognized.

Rik <rik5>
Project Administrator
Wed 08 May 2019 03:35:50 PM UTC, comment #40: 

It would also be fine with me if we detect an exactly singular matrix and just return an appropriately sized sparse matrix filled with Inf.  You could easily run out of memory.  But you should be prepared for that anyway if you are inverting a sparse matrix since that operation doesn't preserve sparsity, correct?

John W. Eaton <jwe>
Project Administrator
Wed 08 May 2019 02:59:53 PM UTC, comment #39: 

What happens if, instead of the divide by zero warnings, you create a +/-Inf or NaN element in the output at the location in the output array where the divide by zero would have occurred?  +/-Inf if the numerator is non-zero and NaN if it is 0 (not present).  I think this change will require some care, because the zero in the denominator could (or would?) be a value not stored in the sparse array.  I've looked at tinverse, but I don't know immediately see how to modify it to do what I'm suggesting.

John W. Eaton <jwe>
Project Administrator
Wed 08 May 2019 02:06:53 PM UTC, comment #38: 

Maybe inv(A) should return a matrix with the same sparsity pattern of A and Inf as non-zero elements?

Marco Caliari <caliari>
Project Member
Wed 08 May 2019 01:55:19 PM UTC, comment #37: 

Marco: Well, that is weirdly inconsistent with the other singular matrices you showed.

John W. Eaton <jwe>
Project Administrator
Wed 08 May 2019 01:31:27 PM UTC, comment #36: 

@jwe: you are right, probably we do not want it. But this is matlab

>> full(inv(sparse(10,10,0)))
Warning: Matrix is singular to working precision.

ans =

   Inf   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   Inf   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   Inf   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   Inf   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   Inf   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   Inf   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   Inf   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN   Inf   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   Inf   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   Inf

Marco Caliari <caliari>
Project Member
Wed 08 May 2019 01:23:18 PM UTC, comment #35: 

Marco, are you saying that any time we detect a division by zero in the sparse inverse function that we should return a sparse matrix of the appropriate dimension that is all Inf?  If so, that will probably often lead to out of memory errors.  Is that really what we want?

John W. Eaton <jwe>
Project Administrator
Wed 08 May 2019 12:19:09 PM UTC, comment #34: 

For a sparse singular matrix, it is not clear to me what matlab does:

>> inv(sparse([1,0;1,0]))
Warning: Matrix is singular to working precision.

ans =

   (1,1)        1
   (2,1)     -Inf
   (2,2)      Inf
>> inv(sparse([1,2;0,0]))
Warning: Matrix is singular to working precision.

ans =

   (1,1)      NaN
   (2,1)      NaN
   (1,2)     -Inf
   (2,2)      Inf

For the same examples, in the full form, matlab gives a matrix of Inf. The same does octave. I would suggest that octave gives a sparse matrix of Inf. There are five cases in which tinverse gives the division by zero error. It should be enough to modify them in such a way that a sparse matrix of Inf is returned.

Marco Caliari <caliari>
Project Member
Wed 08 May 2019 12:05:11 PM UTC, comment #33: 

@jwe: yes, tinverse has to be changed. Please notice that you get the division by zero error also with standard singular matrices

octave:1> s=sparse([1,2;0,0])
s =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 2 [50%])

  (1, 1) ->  1
  (1, 2) ->  2

octave:2> inv(s)
error: inv: division by zero

The triangular U factor has zeros on the diagonal and tinverse requires to invert it.

Marco Caliari <caliari>
Project Member
Wed 08 May 2019 11:50:07 AM UTC, comment #32: 

I'm not worried about the allocations and it seems like always having valid arrays will make the implementation simpler for the mex interface and for calling suitesparse functions.

Making the minimum value for nzmax 1 also seems fine, and if that's what Matlab does, there may be code that relies on it so we might as well do the same.

For the problem with inverse, I guess we just need to fix SparseMatrix::tinverse to generate +/-Inf and NaN values instead of reporting division by zero.  I looked at it and am hopeful that someone else understands that code and the sparse data structure better than I do.

John W. Eaton <jwe>
Project Administrator
Wed 08 May 2019 07:28:49 AM UTC, comment #31: 

@Rik, comment #30:

>> ver
-----------------------------------------------------------------------------------------------------
MATLAB Version: 9.3.0.713579 (R2017b)
-----------------------------------------------------------------------------------------------------
>> mex_sparse
Reported nzmax size when created with 0 elements: 1
mxGetIr returns valid pointer when nzmax == 0
mxGetPr returns valid pointer when nzmax == 0
Reported nzmax size when created with 1 elements: 1
mxGetIr returns valid pointer when nzmax == 1
mxGetPr returns valid pointer when nzmax == 1

I will answer to the final part of comment #29 later.

Marco Caliari <caliari>
Project Member
Tue 07 May 2019 10:04:53 PM UTC, comment #30: 

I'm attaching a file mex_sparse.c that we need someone to compile and run in Matlab to verify what it does.  On looking through the documentation, it seems that nzmax must always be a minimum of 1.  I.e., they always allocate space for at least one element.

Within Matlab, I tried

S = sparse (zeros (2));
nzmax (S)
=> 1

so this minimum is preserved in the interpreter as well, not just the mex interface.

After compiling and running the mex_sparse.c in Octave I get

octave:4> mex_sparse
Reported nzmax size when created with 0 elements: 0
mxGetIr returns valid pointer when nzmax == 0
mxGetPr returns valid pointer when nzmax == 0
Reported nzmax size when created with 1 elements: 1
mxGetIr returns valid pointer when nzmax == 1
mxGetPr returns valid pointer when nzmax == 1

So Octave allows the creation of fully empty (nzmax = 0, nnz = 0) sparse matrices.  The pointer returned from mxGetIr and mxGetPr is not NULL even in this instance.

(file #46880)

Rik <rik5>
Project Administrator
Tue 07 May 2019 09:56:52 PM UTC, comment #29: 

@jwe: This was the approach I was trying to avoid since it requires three changes instead of one, and the overhead of managing small blocks on the heap (https://stackoverflow.com/questions/1087042/c-new-int0-will-it-allocate-memory).

But, if it will make the implementation of the mex interface simpler we could do it.  Since Matlab's documentation is often only loosely coupled to its behavior, we should verify with a sample mex program, that mxGetJc really does return a valid pointer for an empty sparse array.

Attached is a patch that modifies the sparse constructors, but uses C++11 aggregate value initialization '{}' in the mem-init rather than writing a constructor to which uses std::fill_n.  It makes for shorter code without, I think, any loss of clarity.

Regardless, with either patch applied, I get

error: inv: division by zero
error: called from
    tst_inv at line 3 column 1

when trying to take the inverse of an empty matrix.

(file #46879)

Rik <rik5>
Project Administrator
Tue 07 May 2019 08:19:08 PM UTC, comment #28: 

After seeing comment #16, I built umfpackmex from UMFPACK for Octave.  If I call it with

[l, u, p, q] = umfpackmex (sparse (2, 2, 0))

then it errors with the symbolic factorization failed message.  As David says in comment #18, this failure can happen if Ai or Ap are not allocated.  So just before the call to umfpack_dl_qsymbolic, I added

   if (!Ai)
      Ai = (Long *) malloc (0);

and then it works with both Ap and Ai valid pointers and it returns the following result:

octave:1> [l, u, p, q] = umfpackmex (sparse (2, 2, 0))
l =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 2 [50%])

  (1, 1) ->  1
  (2, 2) ->  1

u =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 0 [0%])


p =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 2 [50%])

  (1, 1) ->  1
  (2, 2) ->  1

q =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 2 [50%])

  (2, 1) ->  1
  (1, 2) ->  1

I would also prefer if we didn't have to allocate arrays that will never be used (because NNZ is zero, for example, or the size or the sparse array is 0x0) but I also found that mxGetJc and mxGetIr are documented to return NULL on error but it looks like they should return valid pointers otherwise, even if the lengths of the arrays the point to are zero.  We could arrange for Octave's mxGetJc and mxGetIr functions to fake that behavior, but it might be easiest if we always allocated them in Octave's internals as well, even if there are no elements in the sparse array.

It seems that UMFPACK also will fail if the data pointer is NULL.

The attached change seems to work for me.

(file #46877)

John W. Eaton <jwe>
Project Administrator
Tue 07 May 2019 08:05:38 PM UTC, comment #27: 

I think that is a step in the right direction.  I had made the second change to line 686 and then I was getting a new segfault that was related to the free'ing of memory within UMFPACK library itself.

Rik <rik5>
Project Administrator
Tue 07 May 2019 07:12:34 PM UTC, comment #26: 

@Rik: the change to Ai should be applied, in case, also around line 686 of sparse-lu. But then octave segfaults. Could it be that Ax has to be different from nullptr, too? I found in umfpack documentation

However, nz can be zero, since all singular matrices are
handled correctly. If you attempt to malloc an array of size nz = 0, however, malloc will return
a null pointer which UMFPACK will report as a missing argument. If you malloc an array of size
nz to pass to UMFPACK, make sure that you handle the nz = 0 case correctly (use a size equal to
the maximum of nz and 1, or use a size of nz+1).

Marco Caliari <caliari>
Project Member
Tue 07 May 2019 05:04:34 PM UTC, comment #25: 

Tried the suggestion from comment #21 and UMFPACK still fails at the symbolic factorization for

s = sparse (zeros (2));
inv (s)

and fails with

error: inv: division by zero

for

s = sparse (ones (2));
s = s - s;
inv (s)

I can check later this evening with gdb to see what is happening.

Rik <rik5>
Project Administrator
Tue 07 May 2019 04:53:16 PM UTC, comment #24: 

I want to avoid "the tail wagging the dog".  We want to use UMFPACK, but someone else's decisions on internal code organization should not dictate how Octave works.  For interoperation we should do whatever is necessary at the API join between the two code bases.  I'll test David's suggestion to replace the ridx with the cidx when interfacing to UMFPACK.

Rik <rik5>
Project Administrator
Tue 07 May 2019 04:47:43 PM UTC, comment #23: 

I think this is going to require several steps and several patches to properly fix.  I believe the patch for lu-sparse.cc is still a part of the overall solution.

For starters, what should the internal representation of a sparse, empty matrix be in Octave?  Right now there are two models used simultaneously which is not good coding practice.  The first model is from the Sparse constructor which assigns a nullptr to the data and row indices.  The second model is used by an existing Sparse matrix which becomes empty.  In that case, a pointer is assigned to the row and data indices, but it is expected that they not be used because there is no valid data there.

Either will work, but I think some practical concerns lean towards the first model.  First, it is more intuitive to me that if the matrix is empty, there are no row or data entries.  Second, although there is not a real performance hit on modern computers with lots of memory, why allocate two pointers worth of memory with new() and have the runtime maintain that information when it isn't necessary.  Third, changing over to the second convention will, I think, involve more changes to the code base.  The second representation is, I believe, only used in the change_length routine so a single modification converts all of Octave to the first representation.

Rik <rik5>
Project Administrator
Tue 07 May 2019 03:02:32 PM UTC, comment #22: 

Opps that should read a.cidx(nc. Sorry off by one error. I don't have an octave source tree to be able to check my work before commenting

D.

Anonymous
Tue 07 May 2019 02:52:53 PM UTC, comment #21: 

Proposition for a simple fix.

1/ Rather than use Rik's patch in sparse-lu.cc around line 459, male a change like

// If nz is zero then ridx is probably a nullptr. In this case
// the value of ridx is not important but UMFPACK needs a valid
// pointer in any case, so just set it to cidx that is always
// valid
const octave_idx_type *Ai = (a.cidx(nc+1) == 0 ? a.cidx() : a.ridx ());

This should fix the LU factorization problem

2/ In tinverse (line 718 in dsparse.cc) make the change

      octave_idx_type nz2 = nz == 0 ? 1 ; nz;

and the same in CSparse.cc. This should fix the segfault

I can't test it but this will probably work, even if inelegant

David

Anonymous
Tue 07 May 2019 02:28:40 PM UTC, comment #20: 

Here is another interesting data point

>> [L,U,P,Q] = lu(sparse([],[],[],2,2,1))
L =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 2 [50%])

  (1, 1) ->  1
  (2, 2) ->  1

U =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 0 [0%])


P =

Permutation Matrix

   1   0
   0   1

Q =

Permutation Matrix

   0   1
   1   0

I use the value of NZMAX (set to 1) to create a sparse matrix that although it is empty has valid row, col et data indexes and LU returns a valid factorization. So again setting the sparse column et data indices to a null pointer is not the right thing to do.

David

Anonymous
Tue 07 May 2019 01:04:12 PM UTC, comment #19: 

octave 5.1.1

warning: lu: function may fail when called with less than 4 output arguments and a sparse input
UMFPACK V5.7.6 (May 4, 2016), Control:
    Matrix entry defined as: double
    Int (generic integer) defined as: SuiteSparse_long

    0: print level: 6
    1: dense row parameter:    0.2
        "dense" rows have    > max (16, (0.2)*16*sqrt(n_col) entries)
    2: dense column parameter: 0.2
        "dense" columns have > max (16, (0.2)*16*sqrt(n_row) entries)
    3: pivot tolerance: 0.1
    4: block size for dense matrix kernels: 32
    5: strategy: 0 (auto)
    10: ordering: 1 AMD/COLAMD
    11: singleton filter: enabled
    6: initial allocation ratio: 0.7
    7: max iterative refinement steps: 2
    13: Q fixed during numerical factorization: 1 (yes)
    14: AMD dense row/col parameter:    10
       "dense" rows/columns have > max (16, (10)*sqrt(n)) entries
        Only used if the AMD ordering is used.
    15: diagonal pivot tolerance: 0.001
        Only used if diagonal pivoting is attempted.
    16: scaling: 0 (no)
    17: frontal matrix allocation ratio: 0.5
    18: drop tolerance: 0
    19: AMD and COLAMD aggressive absorption: 1 (yes)

    The following options can only be changed at compile-time:
    8: BLAS library used:  Fortran BLAS.  size of BLAS integer: 4
    compiled for ANSI C
    POSIX C clock_getttime.
    computer/operating system: Linux
    size of int: 4 SuiteSparse_long: 8 Int: 8 pointer: 8 double: 8 Entry: 8 (in bytes)

column-form matrix, n_row 2 n_col 2, nz = 0. ERROR: Ai missing


UMFPACK:  Copyright (c) 2005-2013 by Timothy A. Davis.  All Rights Reserved.


UMFPACK License:  refer to UMFPACK/Doc/License.txt for the license.
   UMFPACK is available under alternate licenses,
   contact T. Davis for details.


Availability: http://www.suitesparse.com
UMFPACK V5.7.6 (May 4, 2016): ERROR: required argument(s) missing

UMFPACK V5.7.6 (May 4, 2016), Info:
    matrix entry defined as:          double
    Int (generic integer) defined as: SuiteSparse_long
    BLAS library used: Fortran BLAS.  size of BLAS integer: 4
    MATLAB:                           no.
    CPU timer:                        POSIX C clock_getttime ( ) routine.
    number of rows in matrix A:       2
    number of columns in matrix A:    2
    memory usage reported in:         16-byte Units
    size of int:                      4 bytes
    size of SuiteSparse_long:         8 bytes
    size of pointer:                  8 bytes
    size of numerical entry:          8 bytes

    strategy used:                    unsymmetric
    ordering used:                    not computed
    symbolic factorization defragmentations:       0

    symbolic/numeric factorization:      upper bound               actual      %
    variable-sized part of Numeric object:

error: lu: sparse_lu: symbolic factorization failed

matlab 9.6.0.1099231 (R2019a) Update 1

pivot threshold: 1
sprealloc in lu(L): 2 1 1 3
Pruned nzs: 2
L nonzeros: 2
U nonzeros: 0

ans =

   All zero sparse: 2×2

Marco Caliari <caliari>
Project Member
Tue 07 May 2019 01:04:01 PM UTC, comment #18: 

Looking at https://github.com/PetterS/SuiteSparse/blob/master/UMFPACK/Source/umfpack_qsymbolic.c the code for qsymbolic contains the lines

    if (!Ai || !Ap || !SymbolicHandle)
    {
        Info [UMFPACK_STATUS] = UMFPACK_ERROR_argument_missing ;
        return (UMFPACK_ERROR_argument_missing) ;
    }

right at the top. So the null pointers in the column index causes UMFPACK to fail. Sorry Rik, the at least the column index needs to be initialized, though the data can be a null pointer

David

Anonymous
Tue 07 May 2019 12:48:58 PM UTC, comment #17: 

@Marco

Can you try

sparms("spumoni", 6); lu(sparse(2,2,0))

under both matlab and octave. UMFPACK is probably configured differently. Doing this under Octave, UMFPACK returns the error "ERROR: Ai missing". So UMFPACK is testing the column indexes, finding this its a null pointer and crashing

David

Anonymous
Tue 07 May 2019 09:59:00 AM UTC, comment #16: 

@Rik: your patch fixes the problem with inv, but I think that the source of the problem is the LU factorization of a sparse null matrix. The result should be

>> [L, U, P, Q] = umfpack (sparse (2, 2, 0))

L =

   (1,1)        1
   (2,2)        1


U =

   All zero sparse: 2×2


P =

   (1,1)        1
   (2,2)        1


Q =

   (2,1)        1
   (1,2)        1

Here umfpack () is the mex interface of UMFPACK, run in matlab. So, it seems that UMFPACK is able to compute the factorization of a sparse null matrix, but not in octave. I fear that your patch may hide the true source of troubles.

Marco Caliari <caliari>
Project Member
Tue 07 May 2019 12:11:03 AM UTC, comment #15: 

Attached is a patch which does what I suggested, using nullptr for r and d, if the sparse matrix is empty.  This actually doesn't have any performance impact because the re-allocation code was being executed anyways.

This stops the crash, although the error is still

error: inv: sparse_lu: symbolic factorization failed

The second part of resolving this is to return a warning and some sort of matrix as Matlab does.

For the full case, I get

octave:1> inv (zeros (2))
warning: matrix singular to machine precision
ans =

   Inf   Inf
   Inf   Inf

So that would be one possibility, to return a warning and the matrix

sparse (Inf (2))

Matlab returns

sparse ([Inf, NaN;
         NaN, Inf])

but that is probably just an artifact of the implementation and needn't be copied exactly unless we want to.

(file #46868)

Rik <rik5>
Project Administrator
Mon 06 May 2019 11:10:31 PM UTC, comment #14: 

Actually, my suggestion to zero out r and d is inefficient.  Octave doesn't automatically cull memory when removing zero elements because doing so every single time would be a burden.  It only re-optimizes memory when a significant number (1/5) of elements have been removed.

Rik <rik5>
Project Administrator
Mon 06 May 2019 10:48:54 PM UTC, comment #13: 

Incidentally, I checked with gdb that setting the pointer r to 0x0 did actually fix the situation and cause umfpack_qsymbolic to fail.

Rik <rik5>
Project Administrator
Mon 06 May 2019 10:47:02 PM UTC, comment #12: 

After much time with gdb, I've found a significant difference.

For a newly created empty sparse matrix the underlying representation is

rep = {d = 0x0, r = 0x0, c = 0x7ff1a415c8f0, nzmx = 0, nrows = 4, ncols = 4,

For the case of C - C it is

rep = {d = 0x7ff1a41817d0, r = 0x7ff1a4315df0, c = 0x7ff1a415c8f0, nzmx = 0, nrows = 4, ncols = 4,

What you'll notice is that the list of row indices is not a null pointer.  When the matrix is empty, because there are no row indices, the call to umfpack_qsymbolic correctly fails.  However, when using C - C, there is apparently leftover garbage in the memory locations pointed to by r and d which cause the symbolic factorization to succeed.

This could be fixed in sparse-lu.cc by checking for an empty matrix at the very start before calling any routines from UMFPACK.

However, I think I might prefer to zero out the data (d) and row (r) pointers when there are no entries.  This seems safer because I don't know where else in the Octave code base we might be calling other libraries with empty matrices.

Thoughts?

Rik <rik5>
Project Administrator
Mon 06 May 2019 07:57:32 PM UTC, comment #11: 

Surprise: [L, U] = lu (C - C) does the right thing (L is the identity, U is null, as in the full case, as in Matlab) but gives the "symbolic factorization failed" message with sparse (2, 2, 0).
So, I see now two problems:

1) why C - C is different from sparse (2, 2, 0)
2) tinverse and dinverse should be adapted for zeros in the diagonal, if we want to invert a matrix by inverting its factors, and the matrix is singular

Marco Caliari <caliari>
Project Member
Mon 06 May 2019 07:21:50 PM UTC, comment #10: 

For what it's worth, I've been single stepping through the code with gdb and the matrix type is marked as Full.

Rik <rik5>
Project Administrator
Mon 06 May 2019 01:24:13 PM UTC, comment #9: 

It is not C-C to be marked as triangular. SparseMatrix::inverse is called in dSparse.cc. The branch if (! mattype.ishermitian ()) is entered, octave::math::sparse_lu<SparseMatrix> fact is computed and, instead of failing, InvL and InvU are tried to be computed. Notice that if you replace C with

C = complex (sparse (ones (2)));

fact is computed, but then you get the nicer message

error: inv: sparse_lu: symbolic factorization failed

from sparse_lu. I don't know if it is related with bug #53390 (it is currently not safe to use UMFPACK to factorize into LU a sparse matrix without column reordering).

@David: nice to hear you again.

Marco Caliari <caliari>
Project Member
Mon 06 May 2019 10:13:55 AM UTC, comment #8: 

David Bateman here (can't be bothered trying to recover my password)

Several things I note here:

The function "tinverse" was written to handle triangular or permuted triangular matrices. Though

C = sparse (ones (2)); C = C - C; matrix_type(C)
>> Full
C = sparse (2,2,0); Matrix_type(C)
>> Full

So something is screwing up in the detection of the matrix type of C-C and its being detcted in the subtraction operation as triangular

I see to way of fixing this. Fix the matrix type detection code to detect "C - C" as Full. In fact an empty sparse matrix could be considered as triangular, so this seems wrong to me. Otherwise the function tinverse should be fixed to allow for empty matrices. In that case dinverse might alos need to be adapted for zeros en the diagonal

D.

Anonymous
Wed 01 May 2019 04:35:37 PM UTC, comment #7: 

Maybe the problem is not in the computation of C-C, but shouldn't be the macro SPARSE_SMSM_BIN_OP_1(R, F, OP, M1, M2) in liboctave/operators/Sparse-op-defs.h performing that operation? I would think so, but it is not. Where is the code which computes C-C? And what is the purpose of SPARSE_SMSM_BIN_OP_1?

Marco Caliari <caliari>
Project Member
Tue 30 Apr 2019 10:00:01 PM UTC, comment #6: 

It is actually the double difference:

octave:1> C = sparse (ones (2))
C =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 4 [100%])

  (1, 1) ->  1
  (2, 1) ->  1
  (1, 2) ->  1
  (2, 2) ->  1

octave:2> a = C - C
a =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 0 [0%])


octave:3> b = a - a
b =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 0 [0%])


octave:4> inv(b)
error: inv: sparse_lu: symbolic factorization failed
octave:4> inv(a)

Thread 10 "QThread" received signal SIGSEGV, Segmentation fault.
...

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Tue 30 Apr 2019 09:34:31 PM UTC, comment #5: 

The observation in comment #4 is probably because (C - C) creates a temporary and the function works on the temporary object.  That might be a clue to the problem.

Rik <rik5>
Project Administrator
Tue 30 Apr 2019 05:43:10 PM UTC, comment #4: 

strange:

>> C = sparse (ones (2)); C = C - C;
>> inv(C-C)
error: inv: sparse_lu: symbolic factorization failed

but inv(C) crashes Octave

Avinoam Kalma <avinoam>
Project Member
Mon 29 Apr 2019 03:26:34 PM UTC, comment #3: 

I think I may have been wrong about (C - C) corrupting things.  It seems that it was only the gdb macro display-sparse-array which has problems, not necessarily the code itself.

Rik <rik5>
Project Administrator
Mon 29 Apr 2019 03:13:50 PM UTC, comment #2: 

I browsed the sources for half an hour... where is the code for sparse_matrix - sparse_matrix?

Marco Caliari <caliari>
Project Member
Sun 28 Apr 2019 09:39:24 PM UTC, comment #1: 

Confirmed.  That is exceedingly odd.  I can also confirm that the segfault goes back through version 4.2.1 so this behavior has probably always been present, but never quite hit upon.

Slightly longer backtrace to show the full code path

#0  0x00007fc47cabdbfb in SparseMatrix::tinverse (this=0x7fc44fff6110, mattype=...,
    info=@0x7fc44fff62f0: 0, rcond=@0x7fc44fff6038: 0, calccond=false)
    at liboctave/array/dSparse.cc:734
#1  0x00007fc47cabf271 in SparseMatrix::inverse (this=0x7fc44fff6420, mattype=...,
    info=@0x7fc44fff62f0: 0, rcond=@0x7fc44fff62f8: 0, calc_cond=true)
    at liboctave/array/dSparse.cc:983
#2  0x00007fc47e4b3500 in Finv (args=..., nargout=0) at libinterp/corefcn/inv.cc:146

Inspecting things in SparseMatrix::inverse, I think (C -C) may be corrupting the dim_vector within the Sparse representation.

Rik <rik5>
Project Administrator
Sun 28 Apr 2019 06:38:51 PM UTC, original submission:  

I was pointed to an Octave crash, when entering the following line into Octave 4.4.1, 5.1.0, or the 6.0 (hg id: 726d945f23de):

C = sparse (ones (2)); C = C - C; inv (C)

The point is, that by computation the 2x2 sparse matrix becomes empty and calling `inv` crashes Octave.  On the other hand, when creating an empty matrix, only an error is thrown:

C = sparse (2, 2, 0); inv (C)
error: inv: sparse_lu: symbolic factorization failed

Matlab R2019a warns in both cases only:

Warning: Matrix is singular to working precision.

Using GDB debugging, I got the following:

Thread 7 "QThread" received signal SIGSEGV, Segmentation fault.
0x00007ffff58f344e in SparseMatrix::tinverse (this=this@entry=0x7fffc9895030, mattype=..., info=@0x7fffc9895240: 0,
    rcond=@0x7fffc9894f68: 0, calccond=false) at liboctave/array/dSparse.cc:734
734               retval.xridx (cx) = i;

Somehow octave::math::sparse_lu<SparseMatrix>

https://hg.savannah.gnu.org/hgweb/octave/file/ad77f3204eda/liboctave/array/dSparse.cc#l976

seems not to catch this case of an empty matrix and does not prevent the following computation.  Maybe I can add a check for this case?

Kai Torben Ohlhus <siko1056>
Project Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #46913:  sparse-singular4.cset added by caliari (9KiB - application/octet-stream)
file #46911:  sparse-singular3.cset added by caliari (10KiB - application/octet-stream)
file #46908:  sparse-singular2.cset added by caliari (8KiB - application/octet-stream)
file #46906:  sparse-singular.cset added by caliari (8KiB - application/octet-stream)
file #46880:  mex_sparse.c added by rik5 (2KiB - text/x-csrc)
file #46879:  sparse_constructor.diffs added by rik5 (1KiB - application/octet-stream)
file #46877:  sparse-diffs.txt added by jwe (1KiB - text/plain)
file #46868:  bug56232.cset added by rik5 (2KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by avinoam (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by siko1056 (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 20 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-05-18 rik5 StatusPatch Reviewed => Fixed
    2019-05-17 siko1056 StatusPatch Submitted => Patch Reviewed
    2019-05-16 caliari Attached File- => Added sparse-singular4.cset, #46913
        StatusFixed => Patch Submitted
    2019-05-15 caliari Attached File- => Added sparse-singular3.cset, #46911
    2019-05-14 caliari Attached File- => Added sparse-singular2.cset, #46908
    2019-05-14 caliari Attached File- => Added sparse-singular.cset, #46906
    2019-05-13 rik5 StatusPatch Submitted => Fixed
        Open/ClosedOpen => Closed
    2019-05-13 rik5 StatusFixed => Patch Submitted
        Open/ClosedClosed => Open
    2019-05-10 siko1056 StatusReady For Test => Fixed
        Open/ClosedOpen => Closed
    2019-05-09 siko1056 StatusPatch Submitted => Ready For Test
    2019-05-07 rik5 Attached File- => Added mex_sparse.c, #46880
    2019-05-07 rik5 Attached File- => Added sparse_constructor.diffs, #46879
    2019-05-07 jwe Attached File- => Added sparse-diffs.txt, #46877
    2019-05-07 rik5 Attached File- => Added bug56232.cset, #46868
        StatusConfirmed => Patch Submitted
    2019-04-28 rik5 StatusNone => Confirmed

    Back to the top


    Powered by Savane 3.5