bugGNU Octave - Bugs: bug #59488, eigs error in dsaupd for zero...

 
 

bug #59488: eigs error in dsaupd for zero matrix

Submitter:  JD <jdbancal>
Submitted:  Thu 19 Nov 2020 09:36:17 PM UTC
   
 
Category:  Octave Function Severity:  2 - Minor
Priority:  3 - Low Item Group:  Regression
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 26 Nov 2020 06:26:38 AM UTC, comment #16: 

Thanks for the overhaul of that function Rik 🙂  The input validation seems indeed a little "grown" over time and might need more attention 😇

I pushed an additional changeset

https://hg.savannah.gnu.org/hgweb/octave/rev/64528a919674

as I no longer saw any output for the OPs example (nargout == 0).

Kai Torben Ohlhus <siko1056>
Group Member
Wed 25 Nov 2020 10:20:33 PM UTC, comment #15: 

I checked in an overhaul of the eigs.m function which also addresses this bug report.  The changeset was made on the development branch here: http://hg.savannah.gnu.org/hgweb/octave/rev/00baf82bf56e

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Wed 25 Nov 2020 09:38:48 PM UTC, comment #14: 

Before a programmer can make use of eigenvalues and eigenvectors, they have to have some understanding of the linear algebra involved.  An all-zeros matrix doesn't represent any useful linear transformation so it won't be encountered in a practical case where a programmer is modeling a real-world phenomenon.  Maybe somebody doing something in pure math would have a zero-transform that changes all inputs in a vector space in to the zero vector, but the pure math folks don't need a numerical computation language for that, they can work this out symbolically with pen and paper.

The all-zeros matrix is a corner case for which it would be nice to return a value, rather than an error, but isn't critical to most Octave users.

But, I have noticed some issues with the documentation, input validation, and naming conventions within eigs.m.  I've got most of a changeset prepared which fixes these items as well as handling the case of an all-zeros matrix.  I plan to check the code in on the development branch later today.

Rik <rik5>
Group administrator
Wed 25 Nov 2020 08:59:26 PM UTC, comment #13: 

Hmm, how will the programmer know that all eigenvalues are zero, is the error message from octave going to say so?

JD <jdbancal>
Fri 20 Nov 2020 05:17:15 PM UTC, comment #12: 

Thanks for testing.  4.4.1 and 5.2.0 are the same thing, just with a nicer decoding of the message.  I checked the code in libinterp and liboctave and it doesn't appear that Octave changes the vector in any way.  This means that the error message from 4.4.1 ("error -9 in dsaupd") is correct, but attributing it to an empty starting vector is incorrect in this case.

There is a dependence on the version of the external ARPACK library.  In my case, I am using 3.5.0+real-2 and it just works.

It seems like ARPACK may be getting more brittle in the types of input matrices it can accept.  That is annoying, but means we at Octave probably have to program defensively around it.  I would suggest the nnz() test in _eigs_.cc.

I would probably put this change on the development branch rather than the about-to-be-released stable branch.  It is a corner case, and there is an easy workaround.  If Octave throws this error then the programmer already knows that all eigenvalues are zero and there is no need for any computation.

Rik <rik5>
Group administrator
Fri 20 Nov 2020 04:51:41 PM UTC, comment #11: 

With 4.2.2:


octave:1> opts.v0 = [1;0;1];
octave:2> eigs (zeros(3), 1, 'lm', opts)
ans = 0


With 4.4.1:


octave:1> opts.v0 = [1;0;1];
octave:2> eigs (zeros(3), 1, 'lm', opts)
error: __eigs__: eigs: error -9 in dsaupd
error: called from
    eigs at line 285 column 18


With 5.2.0:


octave:1> opts.v0 = [1;0;1];
octave:2> eigs (zeros(3), 1, 'lm', opts)
error: __eigs__: eigs: error in dsaupd: Starting vector is zero
This should not happen.  Please, see https://www.gnu.org/software/octave/bugs.html, and file a bug report
error: called from
    eigs at line 285 column 20


All on the same system with the same version of ARPACK installed from the Debian package:


devnull:810> dpkg -l '*arpack*'
Desired=Unknown/Install/Remove/Purge/Hold
| Status=Not/Inst/Conf-files/Unpacked/halF-conf/Half-inst/trig-aWait/Trig-pend
|/ Err?=(none)/Reinst-required (Status,Err: uppercase=bad)
||/ Name                 Version      Architecture Description
+++-====================-============-============-============================================================================
ii  libarpack2:amd64     3.7.0-3      amd64        Fortran77 subroutines to solve large scale eigenvalue problems
ii  libarpack2-dev:amd64 3.7.0-3      amd64        Fortran77 subroutines to solve large scale eigenvalue problems (development)


Also, for the one that succeeds, eig is called instead of _eigs_, so ARPACK isn't even used in that case.

John W. Eaton <jwe>
Group administrator
Fri 20 Nov 2020 04:48:55 PM UTC, comment #10: 

Here is what I get for the code you mentioned:

octave:1> opts.v0 = [1;0;1];
octave:2> eigs (zeros(3), 1, 'lm', opts)
error: eigs: error in dsaupd: Starting vector is zero
This should not happen.  Please, see https://www.gnu.org/software/octave/bugs.html, and file a bug report
error: called from
    eigs at line 298 column 20


JD <jdbancal>
Fri 20 Nov 2020 04:31:17 PM UTC, comment #9: 

Just to try and disambiguate whether there is a problem with Octave or a problem upstream in the way ARPACK now generates its starting vector could someone try


opts.v0 = [1;0;1];
eigs (zeros(3), 1, 'lm', opts)


Rik <rik5>
Group administrator
Fri 20 Nov 2020 04:21:44 PM UTC, comment #8: 

The issue of interacting with Fortran and BLAS/LAPACK libraries when the matrices contain non-finite values is a large one.  If we want to tackle that, we should develop a strategy first.  Right now we have patched things up in various ways as problems have arisen (such as in inverse routines).  The danger here is that the code base grows substantially (more maintenance burden) to contain a lot of if/elseif checking for specific conditions.  And we may forget to port each case over to a newly discovered example of a bad C++/Fortran interaction.

I'm not even particularly in love with the "solution" we have for the inverse of a matrix.  Our current approach is to calculate the condition number and check it for Inf/NaN before calling the main inversion routines.  This is a performance hit if all the user wanted was the inverse.  If try/catch blocks are cheap in C++ maybe we should just be wrapping the external library calls in them and using the catch block to translate the undecipherable message "parameter number 4 is invalid" to something nice like "non-finite value encountered in matrix".

For empty matrices, it probably makes sense to have a shortcut path in the code.  Besides generating the correct result, it will be a performance win as well.  eigs() is designed to be called on sparse matrices and the nnz() function is cheap on sparse matrices.  So


if (nnz () == 0)
 ...
endif




Rik <rik5>
Group administrator
Fri 20 Nov 2020 03:13:43 PM UTC, comment #7: 

In addition to matrices of zeros, do we also need to have special cases for matrices that contain Inf or NaN?  For example, I see the following:


octave> x = sprand (20, 20, 0.1); x(10) = NaN; eigs (x, 1)
DLASCL: parameter number 4 is invalid
error: Fortran procedure terminated by call to XERBLA
error: called from
    eigs at line 285 column 20

octave> x = sprand (20, 20, 0.1); x(10) = Inf; eigs (x, 1)
error: __eigs__: eigs: error in dneupd: DNAUPD did not find any eigenvalues to sufficient accuracy.
This should not happen.  Please, see https://www.gnu.org/software/octave/bugs.html, and file a bug report
error: called from
    eigs at line 285 column 20


John W. Eaton <jwe>
Group administrator
Fri 20 Nov 2020 03:10:11 PM UTC, comment #6: 

It works for me with 4.2.2 but fails with 4.4.1, so the bug is not new.

Bisecting for old versions of Octave on a current system doesn't work very well because of changes in tools and things, so I tried to just run the debugger on my 4.2.2 and 4.4.1 builds and I found that dsaupd is not called for this problem with version 4.2.2.  So I stepped through eigs.m in both versions and found that in 4.2.2, _eigs_ was not called.

Next, I extracted all versions of eigs.m from the repo since the time around version 4.2.2 and tried the example with each one until I found the first that failed:


## list of revisions with changes to eigs.m
hg log scripts/*/eigs.m --no-merges --template "{rev}\n" > /tmp/revs

## extract eigs.m from each of those revisions
mkdir /tmp/eigs
for rev in $(cat /tmp/revs) ; do hg update -r $rev ; cp scripts/sparse/eigs.m /tmp/eigs/$rev-eigs.m ; done

## test each until we find a failure
## in this case, *-eigs.m is sorts correctly
## in chronological order so the first that fails
## should be the first bad revision
cd /tmp/eigs
for f in *-eigs.m ; do echo $f ; cp $f eigs.m ; octave --eval "eigs (zeros (3), 1)"; rm -f eigs.m ; done


The first revision that showed the error for me was


# HG changeset patch
# User Marco Caliari <marco.caliari@univr.it>
# Date 1524560299 -7200
#      Tue Apr 24 10:58:19 2018 +0200
# Branch stable
# Node ID 1a632692a58e3204d8597404cdee8d7200deac88
# Parent  4c98a9e5ce2557e64256a8f02d857157187a80a3
Use eig in eigs when p is equal to matrix dimension n (bug #53719)


However, even the previous version that doesn't fail for zeros(3) will fail for matrices of zeros that are larger than the limit eigs may also use to determine when to call eig instead of _eigs_.  So we may just need to have a special case for matrices that are all zeros.

John W. Eaton <jwe>
Group administrator
Fri 20 Nov 2020 12:48:53 PM UTC, comment #5: 

My version of libarpack2 and libarpack2-dev (on ubuntu 20.04.1) is 3.7.0-3

JD <jdbancal>
Fri 20 Nov 2020 12:09:24 PM UTC, comment #4: 

Maybe a different version of ARPACK?

Markus Mützel <mmuetzel>
Group administrator
Fri 20 Nov 2020 11:38:28 AM UTC, comment #3: 

I get the same error with dev Octave on Mageia 7.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 20 Nov 2020 10:46:47 AM UTC, comment #2: 

Thanks for checking!
I just tried on the latest hg, namely

changeset:   29079:c43e748de388
bookmark:    @
tag:         tip
parent:      29077:b85f72ec844a
parent:      29078:90a8e6fd85c3
user:        John W. Eaton <jwe@octave.org>
date:        Thu Nov 19 18:59:10 2020 -0500
summary:     maint: merge stable to default.

and I obtain the same error. The only configuration options I used are "--without-qt --without-opengl". Any idea why we see obtain results?

JD <jdbancal>
Thu 19 Nov 2020 10:49:19 PM UTC, comment #1: 

Could you try again with a recent tip from the stable branch of Octave?  I am at cset


parent: 29075:dc3ee9616267 tip
 eigs.m: Fix error message when second input is empty matrix (bug #59486)
branch: stable


and this code works fine.

Rik <rik5>
Group administrator
Thu 19 Nov 2020 09:36:17 PM UTC, original submission:  


octave:1> eigs(zeros(3), 1)
error: eigs: error in dsaupd: Starting vector is zero
This should not happen.  Please, see https://www.gnu.org/software/octave/bugs.html, and file a bug report
error: called from
    eigs at line 298 column 20


Used to work on Octave 4.2

JD <jdbancal>

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by siko1056 (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by jdbancal (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-11-25 rik5 Open/ClosedOpen Closed
    2020-11-25 rik5 StatusWorks For Me Fixed
    2020-11-25 rik5 Severity3 - Normal 2 - Minor
        Priority5 - Normal 3 - Low
        Release6.0.92 dev
    2020-11-19 rik5 StatusNone Works For Me

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code